Merge branch 'lammps:develop' into master

This commit is contained in:
Evangelos Voyiatzis
2023-02-24 08:46:27 +02:00
committed by GitHub
18 changed files with 203 additions and 191 deletions

View File

@ -168,7 +168,7 @@ OPT.
* :doc:`pafi <fix_pafi>`
* :doc:`pair <fix_pair>`
* :doc:`phonon <fix_phonon>`
* :doc:`pimd <fix_pimd>`
* :doc:`pimd/nvt <fix_pimd>`
* :doc:`planeforce <fix_planeforce>`
* :doc:`plumed <fix_plumed>`
* :doc:`poems <fix_poems>`

View File

@ -451,7 +451,10 @@ piece of python code:
for j in range(j1-j2,min(twojmax,j1+j2)+1,2):
if (j>=j1): print j1/2.,j2/2.,j/2.
For even twojmax = 2(*m*\ -1), :math:`K = m(m+1)(2m+1)/6`, the *m*\ -th pyramidal number. For odd twojmax = 2 *m*\ -1, :math:`K = m(m+1)(m+2)/3`, twice the *m*\ -th tetrahedral number.
There are :math:`m(m+1)/2` descriptors with last index *j*,
where *m* = :math:`\lfloor j \rfloor + 1`.
Hence, for even *twojmax* = 2(*m*\ -1), :math:`K = m(m+1)(2m+1)/6`, the *m*\ -th pyramidal number,
and for odd *twojmax* = 2 *m*\ -1, :math:`K = m(m+1)(m+2)/3`, twice the *m*\ -th tetrahedral number.
.. note::

View File

@ -320,7 +320,7 @@ accelerated styles exist.
* :doc:`pafi <fix_pafi>` - constrained force averages on hyper-planes to compute free energies (PAFI)
* :doc:`pair <fix_pair>` - access per-atom info from pair styles
* :doc:`phonon <fix_phonon>` - calculate dynamical matrix from MD simulations
* :doc:`pimd <fix_pimd>` - Feynman path integral molecular dynamics
* :doc:`pimd/nvt <fix_pimd>` - Feynman path integral molecular dynamics with Nose-Hoover thermostat
* :doc:`planeforce <fix_planeforce>` - constrain atoms to move in a plane
* :doc:`plumed <fix_plumed>` - wrapper on PLUMED free energy library
* :doc:`poems <fix_poems>` - constrain clusters of atoms to move as coupled rigid bodies

View File

@ -1,6 +1,6 @@
.. index:: fix pimd
.. index:: fix pimd/nvt
fix pimd command
fix pimd/nvt command
================
Syntax
@ -8,10 +8,10 @@ Syntax
.. parsed-literal::
fix ID group-ID pimd keyword value ...
fix ID group-ID pimd/nvt keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* pimd = style name of this fix command
* pimd/nvt = style name of this fix command
* zero or more keyword/value pairs may be appended
* keyword = *method* or *fmass* or *sp* or *temp* or *nhc*
@ -28,11 +28,15 @@ Examples
.. code-block:: LAMMPS
fix 1 all pimd method nmpimd fmass 1.0 sp 2.0 temp 300.0 nhc 4
fix 1 all pimd/nvt method nmpimd fmass 1.0 sp 2.0 temp 300.0 nhc 4
Description
"""""""""""
.. versionchanged:: TBD
Fix pimd was renamed to fix pimd/nvt.
This command performs quantum molecular dynamics simulations based on
the Feynman path integral to include effects of tunneling and
zero-point motion. In this formalism, the isomorphism of a quantum
@ -67,7 +71,7 @@ number of NH thermostats would be 3 x N x P x Nc.
.. note::
This fix implements a complete velocity-verlet integrator
Fix pimd/nvt implements a complete velocity-verlet integrator
combined with NH massive chain thermostat, so no other time
integration fix should be used.
@ -76,31 +80,33 @@ value of *pimd* is standard PIMD. A value of *nmpimd* is for
normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics
(CMD). The difference between the styles is as follows.
In standard PIMD, the value used for a bead's fictitious mass is
arbitrary. A common choice is to use Mi = m/P, which results in the
mass of the entire ring-polymer being equal to the real quantum
particle. But it can be difficult to efficiently integrate the
equations of motion for the stiff harmonic interactions in the ring
polymers.
In standard PIMD, the value used for a bead's fictitious mass is
arbitrary. A common choice is to use Mi = m/P, which results in the
mass of the entire ring-polymer being equal to the real quantum
particle. But it can be difficult to efficiently integrate the
equations of motion for the stiff harmonic interactions in the ring
polymers.
A useful way to resolve this issue is to integrate the equations of
motion in a normal mode representation, using Normal Mode
Path-Integral Molecular Dynamics (NMPIMD) :ref:`(Cao1) <Cao1>`. In NMPIMD,
the NH chains are attached to each normal mode of the ring-polymer and
the fictitious mass of each mode is chosen as Mk = the eigenvalue of
the Kth normal mode for k > 0. The k = 0 mode, referred to as the
zero-frequency mode or centroid, corresponds to overall translation of
the ring-polymer and is assigned the mass of the real particle.
A useful way to resolve this issue is to integrate the equations of
motion in a normal mode representation, using Normal Mode
Path-Integral Molecular Dynamics (NMPIMD) :ref:`(Cao1) <Cao1>`. In
NMPIMD, the NH chains are attached to each normal mode of the
ring-polymer and the fictitious mass of each mode is chosen as Mk =
the eigenvalue of the Kth normal mode for k > 0. The k = 0 mode,
referred to as the zero-frequency mode or centroid, corresponds to
overall translation of the ring-polymer and is assigned the mass of
the real particle.
Motion of the centroid can be effectively uncoupled from the other
normal modes by scaling the fictitious masses to achieve a partial
adiabatic separation. This is called a Centroid Molecular Dynamics
(CMD) approximation :ref:`(Cao2) <Cao2>`. The time-evolution (and resulting
dynamics) of the quantum particles can be used to obtain centroid time
correlation functions, which can be further used to obtain the true
quantum correlation function for the original system. The CMD method
also uses normal modes to evolve the system, except only the k > 0
modes are thermostatted, not the centroid degrees of freedom.
Motion of the centroid can be effectively uncoupled from the other
normal modes by scaling the fictitious masses to achieve a partial
adiabatic separation. This is called a Centroid Molecular Dynamics
(CMD) approximation :ref:`(Cao2) <Cao2>`. The time-evolution (and
resulting dynamics) of the quantum particles can be used to obtain
centroid time correlation functions, which can be further used to
obtain the true quantum correlation function for the original system.
The CMD method also uses normal modes to evolve the system, except
only the k > 0 modes are thermostatted, not the centroid degrees of
freedom.
The keyword *fmass* sets a further scaling factor for the fictitious
masses of beads, which can be used for the Partial Adiabatic CMD
@ -152,47 +158,47 @@ related tasks for each of the partitions, e.g.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This fix writes the state of the Nose/Hoover thermostat over all
Fix pimd/nvt writes the state of the Nose/Hoover thermostat over all
quasi-beads to :doc:`binary restart files <restart>`. See the
:doc:`read_restart <read_restart>` command for info on how to re-specify
a fix in an input script that reads a restart file, so that the
operation of the fix continues in an uninterrupted fashion.
None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix.
are relevant to fix pimd/nvt.
This fix computes a global 3-vector, which can be accessed by various
:doc:`output commands <Howto_output>`. The three quantities in the
global vector are
Fix pimd/nvt computes a global 3-vector, which can be accessed by
various :doc:`output commands <Howto_output>`. The three quantities in
the global vector are:
#. the total spring energy of the quasi-beads,
#. the current temperature of the classical system of ring polymers,
#. the current value of the scalar virial estimator for the kinetic
energy of the quantum system :ref:`(Herman) <Herman>`.
#. the total spring energy of the quasi-beads,
#. the current temperature of the classical system of ring polymers,
#. the current value of the scalar virial estimator for the kinetic
energy of the quantum system :ref:`(Herman) <Herman>`.
The vector values calculated by this fix are "extensive", except for the
The vector values calculated by fix pimd/nvt are "extensive", except for the
temperature, which is "intensive".
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during
No parameter of fix pimd/nvt can be used with the *start/stop* keywords
of the :doc:`run <run>` command. Fix pimd/nvt is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix is part of the REPLICA package. It is only enabled if
LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.
This fix is part of the REPLICA package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Fix pid cannot be used with :doc:`lj units <units>`.
Fix pimd/nvt cannot be used with :doc:`lj units <units>`.
A PIMD simulation can be initialized with a single data file read via
the :doc:`read_data <read_data>` command. However, this means all
quasi-beads in a ring polymer will have identical positions and
velocities, resulting in identical trajectories for all quasi-beads.
To avoid this, users can simply initialize velocities with different
random number seeds assigned to each partition, as defined by the
uloop variable, e.g.
velocities, resulting in identical trajectories for all quasi-beads. To
avoid this, users can simply initialize velocities with different random
number seeds assigned to each partition, as defined by the uloop
variable, e.g.
.. code-block:: LAMMPS
@ -201,8 +207,8 @@ uloop variable, e.g.
Default
"""""""
The keyword defaults are method = pimd, fmass = 1.0, sp = 1.0, temp = 300.0,
and nhc = 2.
The keyword defaults for fix pimd/nvt are method = pimd, fmass = 1.0, sp
= 1.0, temp = 300.0, and nhc = 2.
----------

View File

@ -256,7 +256,7 @@ for command_type, entries in index.items():
print("Total number of style index entries:", total_index)
skip_angle = ('sdk')
skip_fix = ('python', 'NEIGH_HISTORY/omp','acks2/reax','qeq/reax','reax/c/bonds','reax/c/species')
skip_fix = ('python', 'NEIGH_HISTORY/omp','acks2/reax','qeq/reax','reax/c/bonds','reax/c/species', 'pimd')
skip_pair = ('meam/c','lj/sf','reax/c','lj/sdk','lj/sdk/coul/long','lj/sdk/coul/msm')
skip_compute = ('pressure/cylinder')
@ -286,7 +286,7 @@ if counter:
counter = 0
counter += check_style_index("compute", compute, index["compute"], skip=['pressure/cylinder'])
counter += check_style_index("fix", fix, index["fix"], skip=['python','acks2/reax','qeq/reax','reax/c/bonds','reax/c/species'])
counter += check_style_index("fix", fix, index["fix"], skip=['python','acks2/reax','qeq/reax','reax/c/bonds','reax/c/species','pimd'])
counter += check_style_index("angle_style", angle, index["angle_style"], skip=['sdk'])
counter += check_style_index("bond_style", bond, index["bond_style"])
counter += check_style_index("dihedral_style", dihedral, index["dihedral_style"])

View File

@ -17,7 +17,7 @@ timestep 0.001
velocity all create 1.0 1985 rot yes dist gaussian
fix 1 all pimd method nmpimd fmass 1.0 temp 25.0 nhc 4
fix 1 all pimd/nvt method nmpimd fmass 1.0 temp 25.0 nhc 4
thermo_style custom step temp pe etotal pzz f_1[1] f_1[2] f_1[3]
thermo_modify colname f_1[1] espring colname f_1[2] T_ring colname f_1[3] virial

View File

@ -17,7 +17,7 @@ read_data system.data
#read_restart system_${ibead}.rest1
special_bonds charmm
fix 1 all pimd method nmpimd fmass 1.0 temp 300.0 nhc 4 sp 2.0
fix 1 all pimd/nvt method nmpimd fmass 1.0 temp 300.0 nhc 4 sp 2.0
thermo 10
thermo_style custom step temp pe etotal f_1[1] f_1[2] f_1[3]

2
src/.gitignore vendored
View File

@ -1533,6 +1533,8 @@
/fix_mol_swap.h
/fix_pimd.cpp
/fix_pimd.h
/fix_pimd_nvt.cpp
/fix_pimd_nvt.h
/fix_qbmsst.cpp
/fix_qbmsst.h
/fix_qtb.cpp

View File

@ -185,17 +185,17 @@ void AngleTable::allocate()
void AngleTable::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR, "Illegal angle_style command");
if (narg != 2) error->all(FLERR, "Illegal angle_style command: must have 2 arguments");
if (strcmp(arg[0], "linear") == 0)
tabstyle = LINEAR;
else if (strcmp(arg[0], "spline") == 0)
tabstyle = SPLINE;
else
error->all(FLERR, "Unknown table style in angle style table");
error->all(FLERR, "Unknown table style {} in angle style table", arg[0]);
tablength = utils::inumeric(FLERR, arg[1], false, lmp);
if (tablength < 2) error->all(FLERR, "Illegal number of angle table entries");
if (tablength < 2) error->all(FLERR, "Illegal number of angle table entries: {}", arg[1]);
// delete old tables, since cannot just change settings
@ -218,7 +218,7 @@ void AngleTable::settings(int narg, char **arg)
void AngleTable::coeff(int narg, char **arg)
{
if (narg != 3) error->all(FLERR, "Illegal angle_coeff command");
if (narg != 3) error->all(FLERR, "Illegal angle_coeff command: must have 3 arguments");
if (!allocated) allocate();
int ilo, ihi;
@ -234,7 +234,7 @@ void AngleTable::coeff(int narg, char **arg)
// error check on table parameters
if (tb->ninput <= 1) error->one(FLERR, "Invalid angle table length");
if (tb->ninput <= 1) error->one(FLERR, "Invalid angle table length: {}", tb->ninput);
double alo, ahi;
alo = tb->afile[0];
@ -523,7 +523,7 @@ void AngleTable::param_extract(Table *tb, char *line)
} else if (word == "EQ") {
tb->theta0 = DEG2RAD * values.next_double();
} else {
error->one(FLERR, "Invalid keyword in angle table parameters");
error->one(FLERR, "Unknown keyword {} in angle table parameters", word);
}
}
} catch (TokenizerException &e) {

View File

@ -132,7 +132,7 @@ void BondTable::allocate()
void BondTable::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR, "Illegal bond_style command");
if (narg != 2) error->all(FLERR, "Illegal bond_style command: must have 2 arguments");
tabstyle = NONE;
if (strcmp(arg[0], "linear") == 0)
@ -140,10 +140,10 @@ void BondTable::settings(int narg, char **arg)
else if (strcmp(arg[0], "spline") == 0)
tabstyle = SPLINE;
else
error->all(FLERR, "Unknown table style in bond style table");
error->all(FLERR, "Unknown table style {} in bond style table", arg[0]);
tablength = utils::inumeric(FLERR, arg[1], false, lmp);
if (tablength < 2) error->all(FLERR, "Illegal number of bond table entries");
if (tablength < 2) error->all(FLERR, "Illegal number of bond table entries: {}", arg[1]);
// delete old tables, since cannot just change settings
@ -166,23 +166,21 @@ void BondTable::settings(int narg, char **arg)
void BondTable::coeff(int narg, char **arg)
{
if (narg != 3) error->all(FLERR, "Illegal bond_coeff command");
if (narg != 3) error->all(FLERR, "Illegal bond_coeff command: must have 3 arguments");
if (!allocated) allocate();
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error);
int me;
MPI_Comm_rank(world, &me);
tables = (Table *) memory->srealloc(tables, (ntables + 1) * sizeof(Table), "bond:tables");
Table *tb = &tables[ntables];
null_table(tb);
if (me == 0) read_table(tb, arg[1], arg[2]);
if (comm->me == 0) read_table(tb, arg[1], arg[2]);
bcast_table(tb);
// error check on table parameters
if (tb->ninput <= 1) error->one(FLERR, "Invalid bond table length");
if (tb->ninput <= 1) error->all(FLERR, "Invalid bond table length: {}", tb->ninput);
tb->lo = tb->rfile[0];
tb->hi = tb->rfile[tb->ninput - 1];
@ -483,7 +481,7 @@ void BondTable::param_extract(Table *tb, char *line)
} else if (word == "EQ") {
tb->r0 = values.next_double();
} else {
error->one(FLERR, "Invalid keyword in bond table parameters");
error->one(FLERR, "Unknown keyword {} in bond table parameters", word);
}
}
} catch (TokenizerException &e) {
@ -507,9 +505,9 @@ void BondTable::bcast_table(Table *tb)
int me;
MPI_Comm_rank(world, &me);
if (me > 0) {
memory->create(tb->rfile, tb->ninput, "angle:rfile");
memory->create(tb->efile, tb->ninput, "angle:efile");
memory->create(tb->ffile, tb->ninput, "angle:ffile");
memory->create(tb->rfile, tb->ninput, "bond:rfile");
memory->create(tb->efile, tb->ninput, "bond:efile");
memory->create(tb->ffile, tb->ninput, "bond:ffile");
}
MPI_Bcast(tb->rfile, tb->ninput, MPI_DOUBLE, 0, world);

View File

@ -69,7 +69,6 @@ class BondTable : public Bond {
double splint(double *, double *, double *, int, double);
void uf_lookup(int, double, double &, double &);
void u_lookup(int, double, double &);
};
} // namespace LAMMPS_NS

View File

@ -697,15 +697,15 @@ void DihedralTable::allocate()
void DihedralTable::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Illegal dihedral_style command");
if (narg != 2) error->all(FLERR,"Illegal dihedral_style command: must have 2 arguments");
if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR;
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
else error->all(FLERR,"Unknown table style in dihedral style table");
else error->all(FLERR,"Unknown table style {} in dihedral style table", arg[0]);
tablength = utils::inumeric(FLERR,arg[1],false,lmp);
if (tablength < 3)
error->all(FLERR,"Illegal number of dihedral table entries");
error->all(FLERR,"Illegal number of dihedral table entries: {}", arg[1]);
// delete old tables, since cannot just change settings
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
@ -727,7 +727,7 @@ void DihedralTable::settings(int narg, char **arg)
void DihedralTable::coeff(int narg, char **arg)
{
if (narg != 3) error->all(FLERR,"Incorrect args for dihedral coefficients");
if (narg != 3) error->all(FLERR,"Illegal dihedral_coeff command: must have 3 arguments");
if (!allocated) allocate();
int ilo,ihi;
@ -746,10 +746,9 @@ void DihedralTable::coeff(int narg, char **arg)
// --- and resolve issues with periodicity ---
if (tb->ninput < 2)
error->all(FLERR,"Invalid dihedral table length: {}",arg[2]);
error->all(FLERR,"Invalid dihedral table length: {}", arg[2]);
else if ((tb->ninput == 2) && (tabstyle == SPLINE))
error->all(FLERR,"Invalid dihedral spline table length: {} "
"(Try linear)",arg[2]);
error->all(FLERR,"Invalid dihedral spline table length: {} (Try linear)",arg[2]);
// check for monotonicity
for (int i=0; i < tb->ninput-1; i++) {
@ -767,12 +766,10 @@ void DihedralTable::coeff(int narg, char **arg)
double phihi = tb->phifile[tb->ninput-1];
if (tb->use_degrees) {
if ((phihi - philo) >= 360)
error->all(FLERR,"Dihedral table angle range must be < 360 "
"degrees ({}).",arg[2]);
error->all(FLERR,"Dihedral table angle range must be < 360 degrees ({}).",arg[2]);
} else {
if ((phihi - philo) >= MY_2PI)
error->all(FLERR,"Dihedral table angle range must be < 2*PI "
"radians ({}).",arg[2]);
error->all(FLERR,"Dihedral table angle range must be < 2*PI radians ({}).",arg[2]);
}
// convert phi from degrees to radians
@ -1200,7 +1197,7 @@ void DihedralTable::param_extract(Table *tb, char *line)
checkU_fname = values.next_string();
} else if (word == "CHECKF") {
checkF_fname = values.next_string();
} else error->one(FLERR,"Invalid keyword in dihedral angle table parameters ({})", word);
} else error->one(FLERR,"Unknown keyword {} in dihedral table parameters", word);
}
} catch (TokenizerException &e) {
error->one(FLERR, e.what());

View File

@ -88,12 +88,13 @@ void AtomVecPeri::grow_pointers()
void AtomVecPeri::create_atom_post(int ilocal)
{
const auto xinit = atom->x;
vfrac[ilocal] = 1.0;
rmass[ilocal] = 1.0;
s0[ilocal] = DBL_MAX;
x0[ilocal][0] = x[ilocal][0];
x0[ilocal][1] = x[ilocal][1];
x0[ilocal][2] = x[ilocal][2];
x0[ilocal][0] = xinit[ilocal][0];
x0[ilocal][1] = xinit[ilocal][1];
x0[ilocal][2] = xinit[ilocal][2];
}
/* ----------------------------------------------------------------------
@ -103,10 +104,11 @@ void AtomVecPeri::create_atom_post(int ilocal)
void AtomVecPeri::data_atom_post(int ilocal)
{
const auto xinit = atom->x;
s0[ilocal] = DBL_MAX;
x0[ilocal][0] = x[ilocal][0];
x0[ilocal][1] = x[ilocal][1];
x0[ilocal][2] = x[ilocal][2];
x0[ilocal][0] = xinit[ilocal][0];
x0[ilocal][1] = xinit[ilocal][1];
x0[ilocal][2] = xinit[ilocal][2];
if (rmass[ilocal] <= 0.0) error->one(FLERR, "Invalid mass in Atoms section of data file");
}

View File

@ -51,6 +51,9 @@ lmpinstalledpkgs.h
lmpgitversion.h
mliap_model_python_couple.cpp
mliap_model_python_couple.h
# renamed on 23 February 2023
fix_pimd.cpp
fix_pimd.h
# removed on 20 January 2023
atom_vec_mesont.cpp
atom_vec_mesont.h

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Package FixPIMD
Package FixPIMDNVT
Purpose Quantum Path Integral Algorithm for Quantum Chemistry
Copyright Voth Group @ University of Chicago
Authors Chris Knight & Yuxing Peng (yuxing at uchicago.edu)
@ -21,7 +21,7 @@
Version 1.0
------------------------------------------------------------------------- */
#include "fix_pimd.h"
#include "fix_pimd_nvt.h"
#include "atom.h"
#include "comm.h"
@ -44,7 +44,7 @@ enum { PIMD, NMPIMD, CMD };
/* ---------------------------------------------------------------------- */
FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
FixPIMDNVT::FixPIMDNVT(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
max_nsend = 0;
tag_send = nullptr;
@ -85,27 +85,27 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
else if (strcmp(arg[i + 1], "cmd") == 0)
method = CMD;
else
error->universe_all(FLERR, fmt::format("Unknown method parameter {} for fix pimd",
arg[i + 1]));
error->universe_all(
FLERR, fmt::format("Unknown method parameter {} for fix pimd/nvt", arg[i + 1]));
} else if (strcmp(arg[i], "fmass") == 0) {
fmass = utils::numeric(FLERR, arg[i + 1], false, lmp);
if ((fmass < 0.0) || (fmass > np))
error->universe_all(FLERR, fmt::format("Invalid fmass value {} for fix pimd", fmass));
error->universe_all(FLERR, fmt::format("Invalid fmass value {} for fix pimd/nvt", fmass));
} else if (strcmp(arg[i], "sp") == 0) {
sp = utils::numeric(FLERR, arg[i + 1], false, lmp);
if (sp < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd");
if (sp < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd/nvt");
} else if (strcmp(arg[i], "temp") == 0) {
nhc_temp = utils::numeric(FLERR, arg[i + 1], false, lmp);
if (nhc_temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd");
if (nhc_temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd/nvt");
} else if (strcmp(arg[i], "nhc") == 0) {
nhc_nchain = utils::inumeric(FLERR, arg[i + 1], false, lmp);
if (nhc_nchain < 2) error->universe_all(FLERR, "Invalid nhc value for fix pimd");
if (nhc_nchain < 2) error->universe_all(FLERR, "Invalid nhc value for fix pimd/nvt");
} else
error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix pimd", arg[i]));
error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix pimd/nvt", arg[i]));
}
if (strcmp(update->unit_style, "lj") == 0)
error->all(FLERR, "Fix pimd does not support lj units");
error->all(FLERR, "Fix pimd/nvt does not support lj units");
/* Initiation */
@ -138,7 +138,7 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
/* ---------------------------------------------------------------------- */
FixPIMD::~FixPIMD()
FixPIMDNVT::~FixPIMDNVT()
{
delete[] mass;
atom->delete_callback(id, Atom::GROW);
@ -170,7 +170,7 @@ FixPIMD::~FixPIMD()
/* ---------------------------------------------------------------------- */
int FixPIMD::setmask()
int FixPIMDNVT::setmask()
{
int mask = 0;
mask |= POST_FORCE;
@ -181,13 +181,13 @@ int FixPIMD::setmask()
/* ---------------------------------------------------------------------- */
void FixPIMD::init()
void FixPIMDNVT::init()
{
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR, "Fix pimd requires an atom map, see atom_modify");
error->universe_all(FLERR, "Fix pimd/nvt requires an atom map, see atom_modify");
if (universe->me == 0 && universe->uscreen)
fprintf(universe->uscreen, "Fix pimd initializing Path-Integral ...\n");
fprintf(universe->uscreen, "Fix pimd/nvt initializing Path-Integral ...\n");
// prepare the constants
@ -222,7 +222,7 @@ void FixPIMD::init()
fbond = -_fbond * force->mvv2e;
if (universe->me == 0)
printf("Fix pimd -P/(beta^2 * hbar^2) = %20.7lE (kcal/mol/A^2)\n\n", fbond);
utils::logmesg(lmp, "Fix pimd/nvt -P/(beta^2 * hbar^2) = {:20.7e} (kcal/mol/A^2)\n\n", fbond);
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
@ -241,7 +241,7 @@ void FixPIMD::init()
/* ---------------------------------------------------------------------- */
void FixPIMD::setup(int vflag)
void FixPIMDNVT::setup(int vflag)
{
if (universe->me == 0 && universe->uscreen)
fprintf(universe->uscreen, "Setting up Path-Integral ...\n");
@ -251,7 +251,7 @@ void FixPIMD::setup(int vflag)
/* ---------------------------------------------------------------------- */
void FixPIMD::initial_integrate(int /*vflag*/)
void FixPIMDNVT::initial_integrate(int /*vflag*/)
{
nhc_update_v();
nhc_update_x();
@ -259,14 +259,14 @@ void FixPIMD::initial_integrate(int /*vflag*/)
/* ---------------------------------------------------------------------- */
void FixPIMD::final_integrate()
void FixPIMDNVT::final_integrate()
{
nhc_update_v();
}
/* ---------------------------------------------------------------------- */
void FixPIMD::post_force(int /*flag*/)
void FixPIMDNVT::post_force(int /*flag*/)
{
for (int i = 0; i < atom->nlocal; i++)
for (int j = 0; j < 3; j++) atom->f[i][j] /= np;
@ -293,7 +293,7 @@ void FixPIMD::post_force(int /*flag*/)
Nose-Hoover Chains
------------------------------------------------------------------------- */
void FixPIMD::nhc_init()
void FixPIMDNVT::nhc_init()
{
double tau = 1.0 / omega_np;
double KT = force->boltz * nhc_temp;
@ -309,7 +309,7 @@ void FixPIMD::nhc_init()
nhc_eta_dotdot[i][ichain] = 0.0;
nhc_eta_mass[i][ichain] = mass0;
if ((method == CMD || method == NMPIMD) && universe->iworld == 0)
; // do nothing
; // do nothing
else
nhc_eta_mass[i][ichain] *= fmass;
}
@ -334,7 +334,7 @@ void FixPIMD::nhc_init()
/* ---------------------------------------------------------------------- */
void FixPIMD::nhc_update_x()
void FixPIMDNVT::nhc_update_x()
{
int n = atom->nlocal;
double **x = atom->x;
@ -359,7 +359,7 @@ void FixPIMD::nhc_update_x()
/* ---------------------------------------------------------------------- */
void FixPIMD::nhc_update_v()
void FixPIMDNVT::nhc_update_v()
{
int n = atom->nlocal;
int *type = atom->type;
@ -447,14 +447,14 @@ void FixPIMD::nhc_update_v()
Normal Mode PIMD
------------------------------------------------------------------------- */
void FixPIMD::nmpimd_init()
void FixPIMDNVT::nmpimd_init()
{
memory->create(M_x2xp, np, np, "fix_feynman:M_x2xp");
memory->create(M_xp2x, np, np, "fix_feynman:M_xp2x");
memory->create(M_f2fp, np, np, "fix_feynman:M_f2fp");
memory->create(M_fp2f, np, np, "fix_feynman:M_fp2f");
lam = (double *) memory->smalloc(sizeof(double) * np, "FixPIMD::lam");
lam = (double *) memory->smalloc(sizeof(double) * np, "pimd_nvt:lam");
// Set up eigenvalues
@ -505,7 +505,7 @@ void FixPIMD::nmpimd_init()
/* ---------------------------------------------------------------------- */
void FixPIMD::nmpimd_fill(double **ptr)
void FixPIMDNVT::nmpimd_fill(double **ptr)
{
comm_ptr = ptr;
comm->forward_comm(this);
@ -513,7 +513,7 @@ void FixPIMD::nmpimd_fill(double **ptr)
/* ---------------------------------------------------------------------- */
void FixPIMD::nmpimd_transform(double **src, double **des, double *vector)
void FixPIMDNVT::nmpimd_transform(double **src, double **des, double *vector)
{
int n = atom->nlocal;
int m = 0;
@ -528,7 +528,7 @@ void FixPIMD::nmpimd_transform(double **src, double **des, double *vector)
/* ---------------------------------------------------------------------- */
void FixPIMD::spring_force()
void FixPIMDNVT::spring_force()
{
spring_energy = 0.0;
@ -576,7 +576,7 @@ void FixPIMD::spring_force()
Comm operations
------------------------------------------------------------------------- */
void FixPIMD::comm_init()
void FixPIMDNVT::comm_init()
{
if (size_plan) {
delete[] plan_send;
@ -634,17 +634,17 @@ void FixPIMD::comm_init()
/* ---------------------------------------------------------------------- */
void FixPIMD::comm_exec(double **ptr)
void FixPIMDNVT::comm_exec(double **ptr)
{
int nlocal = atom->nlocal;
if (nlocal > max_nlocal) {
max_nlocal = nlocal + 200;
int size = sizeof(double) * max_nlocal * 3;
buf_recv = (double *) memory->srealloc(buf_recv, size, "FixPIMD:x_recv");
buf_recv = (double *) memory->srealloc(buf_recv, size, "FixPIMDNVT:x_recv");
for (int i = 0; i < np; i++)
buf_beads[i] = (double *) memory->srealloc(buf_beads[i], size, "FixPIMD:x_beads[i]");
buf_beads[i] = (double *) memory->srealloc(buf_beads[i], size, "FixPIMDNVT:x_beads[i]");
}
// copy local positions
@ -666,9 +666,9 @@ void FixPIMD::comm_exec(double **ptr)
if (nsend > max_nsend) {
max_nsend = nsend + 200;
tag_send =
(tagint *) memory->srealloc(tag_send, sizeof(tagint) * max_nsend, "FixPIMD:tag_send");
buf_send =
(double *) memory->srealloc(buf_send, sizeof(double) * max_nsend * 3, "FixPIMD:x_send");
(tagint *) memory->srealloc(tag_send, sizeof(tagint) * max_nsend, "FixPIMDNVT:tag_send");
buf_send = (double *) memory->srealloc(buf_send, sizeof(double) * max_nsend * 3,
"FixPIMDNVT:x_send");
}
// send tags
@ -709,7 +709,7 @@ void FixPIMD::comm_exec(double **ptr)
/* ---------------------------------------------------------------------- */
int FixPIMD::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
int FixPIMDNVT::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int i, j, m;
@ -727,7 +727,7 @@ int FixPIMD::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
/* ---------------------------------------------------------------------- */
void FixPIMD::unpack_forward_comm(int n, int first, double *buf)
void FixPIMDNVT::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
@ -744,28 +744,28 @@ void FixPIMD::unpack_forward_comm(int n, int first, double *buf)
Memory operations
------------------------------------------------------------------------- */
double FixPIMD::memory_usage()
double FixPIMDNVT::memory_usage()
{
return (double) atom->nmax * size_peratom_cols * sizeof(double);
}
/* ---------------------------------------------------------------------- */
void FixPIMD::grow_arrays(int nmax)
void FixPIMDNVT::grow_arrays(int nmax)
{
if (nmax == 0) return;
int count = nmax * 3;
memory->grow(array_atom, nmax, size_peratom_cols, "FixPIMD::array_atom");
memory->grow(nhc_eta, count, nhc_nchain, "FixPIMD::nh_eta");
memory->grow(nhc_eta_dot, count, nhc_nchain + 1, "FixPIMD::nh_eta_dot");
memory->grow(nhc_eta_dotdot, count, nhc_nchain, "FixPIMD::nh_eta_dotdot");
memory->grow(nhc_eta_mass, count, nhc_nchain, "FixPIMD::nh_eta_mass");
memory->grow(array_atom, nmax, size_peratom_cols, "pimd_nvt:array_atom");
memory->grow(nhc_eta, count, nhc_nchain, "pimd_nvt:nh_eta");
memory->grow(nhc_eta_dot, count, nhc_nchain + 1, "pimd_nvt:nh_eta_dot");
memory->grow(nhc_eta_dotdot, count, nhc_nchain, "pimd_nvt:nh_eta_dotdot");
memory->grow(nhc_eta_mass, count, nhc_nchain, "pimd_nvt:nh_eta_mass");
}
/* ---------------------------------------------------------------------- */
void FixPIMD::copy_arrays(int i, int j, int /*delflag*/)
void FixPIMDNVT::copy_arrays(int i, int j, int /*delflag*/)
{
int i_pos = i * 3;
int j_pos = j * 3;
@ -778,7 +778,7 @@ void FixPIMD::copy_arrays(int i, int j, int /*delflag*/)
/* ---------------------------------------------------------------------- */
int FixPIMD::pack_exchange(int i, double *buf)
int FixPIMDNVT::pack_exchange(int i, double *buf)
{
int offset = 0;
int pos = i * 3;
@ -797,7 +797,7 @@ int FixPIMD::pack_exchange(int i, double *buf)
/* ---------------------------------------------------------------------- */
int FixPIMD::unpack_exchange(int nlocal, double *buf)
int FixPIMDNVT::unpack_exchange(int nlocal, double *buf)
{
int offset = 0;
int pos = nlocal * 3;
@ -816,7 +816,7 @@ int FixPIMD::unpack_exchange(int nlocal, double *buf)
/* ---------------------------------------------------------------------- */
int FixPIMD::pack_restart(int i, double *buf)
int FixPIMDNVT::pack_restart(int i, double *buf)
{
int offset = 0;
int pos = i * 3;
@ -837,7 +837,7 @@ int FixPIMD::pack_restart(int i, double *buf)
/* ---------------------------------------------------------------------- */
void FixPIMD::unpack_restart(int nlocal, int nth)
void FixPIMDNVT::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
@ -864,21 +864,21 @@ void FixPIMD::unpack_restart(int nlocal, int nth)
/* ---------------------------------------------------------------------- */
int FixPIMD::maxsize_restart()
int FixPIMDNVT::maxsize_restart()
{
return size_peratom_cols + 1;
}
/* ---------------------------------------------------------------------- */
int FixPIMD::size_restart(int /*nlocal*/)
int FixPIMDNVT::size_restart(int /*nlocal*/)
{
return size_peratom_cols + 1;
}
/* ---------------------------------------------------------------------- */
double FixPIMD::compute_vector(int n)
double FixPIMDNVT::compute_vector(int n)
{
if (n == 0) { return spring_energy; }
if (n == 1) { return t_sys; }

View File

@ -13,21 +13,22 @@
#ifdef FIX_CLASS
// clang-format off
FixStyle(pimd,FixPIMD);
FixStyle(pimd,FixPIMDNVT);
FixStyle(pimd/nvt,FixPIMDNVT);
// clang-format on
#else
#ifndef FIX_PIMD_H
#define FIX_PIMD_H
#ifndef FIX_PIMD_NVT_H
#define FIX_PIMD_NVT_H
#include "fix.h"
namespace LAMMPS_NS {
class FixPIMD : public Fix {
class FixPIMDNVT : public Fix {
public:
FixPIMD(class LAMMPS *, int, char **);
~FixPIMD() override;
FixPIMDNVT(class LAMMPS *, int, char **);
~FixPIMDNVT() override;
int setmask() override;

View File

@ -2,18 +2,17 @@ BootStrap: docker
From: ubuntu:18.04
%environment
export PATH=/usr/lib/ccache:/usr/local/cuda-11.7/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
export CUDADIR=/usr/local/cuda-11.7
export CUDA_PATH=/usr/local/cuda-11.7
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-11.7/lib64:/opt/rocm/lib:/opt/rocm-5.1.3/llvm/lib
export LIBRARY_PATH=/usr/local/cuda-11.7/lib64/stubs
export PATH=/usr/lib/ccache:/usr/local/cuda-12.0/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
export CUDADIR=/usr/local/cuda-12.0
export CUDA_PATH=/usr/local/cuda-12.0
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-12.0/lib64:/opt/rocm/lib:/opt/rocm-5.4.3/llvm/lib
export LIBRARY_PATH=/usr/local/cuda-12.0/lib64/stubs
%post
export DEBIAN_FRONTEND=noninteractive
apt-get update
apt-get upgrade --no-install-recommends -y
apt-get install -y --no-install-recommends curl wget libnuma-dev gnupg ca-certificates
apt-get install --no-install-recommends -y software-properties-common
###########################################################################
# Latest CMake
@ -25,10 +24,10 @@ From: ubuntu:18.04
apt install -y cmake
###########################################################################
# ROCm 4.5
# ROCm 5.4.3
###########################################################################
wget https://repo.radeon.com/amdgpu-install/22.10.3/ubuntu/focal/amdgpu-install_22.10.3.50103-1_all.deb
apt-get install -y ./amdgpu-install_22.10.3.50103-1_all.deb
wget https://repo.radeon.com/amdgpu-install/5.4.3/ubuntu/focal/amdgpu-install_5.4.50403-1_all.deb
apt-get install -y ./amdgpu-install_5.4.50403-1_all.deb
apt-get update
apt-get install --no-install-recommends -y \
@ -48,6 +47,7 @@ From: ubuntu:18.04
###########################################################################
# Common Software
###########################################################################
apt-get install --no-install-recommends -y software-properties-common
add-apt-repository ppa:openkim/latest -y
apt-get update
apt-get install --no-install-recommends -y \
@ -120,13 +120,13 @@ From: ubuntu:18.04
# CUDA
###########################################################################
wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/cuda-ubuntu1804.pin
mv cuda-ubuntu1804.pin /etc/apt/preferences.d/cuda-repository-pin-600
apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/3bf863cc.pub
add-apt-repository "deb http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/ /"
wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu1804/x86_64/cuda-keyring_1.0-1_all.deb
dpkg -i cuda-keyring_1.0-1_all.deb
rm cuda-keyring_1.0-1_all.deb
apt-get update
export CUDA_PKG_VERSION=11.7
export CUDA_PKG_VERSION=12.0
export CUDA_DRV_VERSION=525
apt-get install -y --no-install-recommends \
cuda-libraries-${CUDA_PKG_VERSION} \
@ -135,7 +135,8 @@ From: ubuntu:18.04
cuda-minimal-build-${CUDA_PKG_VERSION} \
cuda-compat-$CUDA_PKG_VERSION \
libcublas-${CUDA_PKG_VERSION} \
libcublas-dev-${CUDA_PKG_VERSION}
libcublas-dev-${CUDA_PKG_VERSION} \
libnvidia-compute-${CUDA_DRV_VERSION}
# add missing symlink
ln -s /usr/local/cuda-${CUDA_PKG_VERSION}/lib64/stubs/libcuda.so /usr/local/cuda-${CUDA_PKG_VERSION}/lib64/stubs/libcuda.so.1
@ -144,8 +145,8 @@ From: ubuntu:18.04
# NVIDIA OpenCL
###########################################################################
mkdir -p /etc/OpenCL/vendors
echo "libnvidia-opencl.so.1" > /etc/OpenCL/vendors/nvidia.icd
test ! -d /etc/OpenCL/vendors && mkdir -p /etc/OpenCL/vendors && \
echo "libnvidia-opencl.so.1" > /etc/OpenCL/vendors/nvidia.icd
###########################################################################
@ -182,7 +183,7 @@ From: ubuntu:18.04
tar -xzf plumed.tar.gz
cd plumed-${PLUMED_PKG_VERSION}
./configure --disable-doc --prefix=/usr
make
make -j 4
make install
cd ../../
rm -rvf plumed

View File

@ -2,11 +2,11 @@ BootStrap: docker
From: ubuntu:20.04
%environment
export PATH=/usr/lib/ccache:/usr/local/cuda-11.7/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
export CUDADIR=/usr/local/cuda-11.7
export CUDA_PATH=/usr/local/cuda-11.7
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-11.7/lib64:/opt/rocm/lib:/opt/rocm-5.1.3/llvm/lib
export LIBRARY_PATH=/usr/local/cuda-11.7/lib64/stubs
export PATH=/usr/lib/ccache:/usr/local/cuda-12.0/bin:${PATH}:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin/x86_64
export CUDADIR=/usr/local/cuda-12.0
export CUDA_PATH=/usr/local/cuda-12.0
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-12.0/lib64:/opt/rocm/lib:/opt/rocm-5.4.3/llvm/lib
export LIBRARY_PATH=/usr/local/cuda-12.0/lib64/stubs
%post
export DEBIAN_FRONTEND=noninteractive
apt-get update
@ -15,10 +15,10 @@ From: ubuntu:20.04
apt-get install -y --no-install-recommends curl wget libnuma-dev gnupg ca-certificates
###########################################################################
# ROCm 5.1.3
# ROCm 5.4.3
###########################################################################
wget https://repo.radeon.com/amdgpu-install/22.10.3/ubuntu/focal/amdgpu-install_22.10.3.50103-1_all.deb
apt-get install -y ./amdgpu-install_22.10.3.50103-1_all.deb
wget https://repo.radeon.com/amdgpu-install/5.4.3/ubuntu/focal/amdgpu-install_5.4.50403-1_all.deb
apt-get install -y ./amdgpu-install_5.4.50403-1_all.deb
apt-get update
apt-get install --no-install-recommends -y \
@ -107,13 +107,13 @@ From: ubuntu:20.04
# CUDA
###########################################################################
wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/cuda-ubuntu2004.pin
mv cuda-ubuntu2004.pin /etc/apt/preferences.d/cuda-repository-pin-600
apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/3bf863cc.pub
add-apt-repository "deb http://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/ /"
wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/cuda-keyring_1.0-1_all.deb
dpkg -i cuda-keyring_1.0-1_all.deb
rm cuda-keyring_1.0-1_all.deb
apt-get update
export CUDA_PKG_VERSION=11.5
export CUDA_PKG_VERSION=12.0
export CUDA_DRV_VERSION=525
apt-get install -y --no-install-recommends \
cuda-libraries-${CUDA_PKG_VERSION} \
@ -122,7 +122,8 @@ From: ubuntu:20.04
cuda-minimal-build-${CUDA_PKG_VERSION} \
cuda-compat-$CUDA_PKG_VERSION \
libcublas-${CUDA_PKG_VERSION} \
libcublas-dev-${CUDA_PKG_VERSION}
libcublas-dev-${CUDA_PKG_VERSION} \
libnvidia-compute-${CUDA_DRV_VERSION}
# add missing symlink
ln -s /usr/local/cuda-${CUDA_PKG_VERSION}/lib64/stubs/libcuda.so /usr/local/cuda-${CUDA_PKG_VERSION}/lib64/stubs/libcuda.so.1
@ -131,9 +132,8 @@ From: ubuntu:20.04
# NVIDIA OpenCL
###########################################################################
mkdir -p /etc/OpenCL/vendors
echo "libnvidia-opencl.so.1" > /etc/OpenCL/vendors/nvidia.icd
test ! -d /etc/OpenCL/vendors && mkdir -p /etc/OpenCL/vendors && \
echo "libnvidia-opencl.so.1" > /etc/OpenCL/vendors/nvidia.icd
###########################################################################
# KIM-API
@ -169,7 +169,7 @@ From: ubuntu:20.04
tar -xzf plumed.tar.gz
cd plumed-${PLUMED_PKG_VERSION}
./configure --disable-doc --prefix=/usr
make
make -j 4
make install
cd ../../
rm -rvf plumed