rename fix pimd to fix pimd/nvt (with backward compatibility)

This commit is contained in:
Axel Kohlmeyer
2023-02-23 12:26:42 -05:00
parent e422ae9477
commit 6ae177f37e
10 changed files with 129 additions and 117 deletions

View File

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

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:`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:`pair <fix_pair>` - access per-atom info from pair styles
* :doc:`phonon <fix_phonon>` - calculate dynamical matrix from MD simulations * :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:`planeforce <fix_planeforce>` - constrain atoms to move in a plane
* :doc:`plumed <fix_plumed>` - wrapper on PLUMED free energy library * :doc:`plumed <fix_plumed>` - wrapper on PLUMED free energy library
* :doc:`poems <fix_poems>` - constrain clusters of atoms to move as coupled rigid bodies * :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 Syntax
@ -8,10 +8,10 @@ Syntax
.. parsed-literal:: .. 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 * 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 * zero or more keyword/value pairs may be appended
* keyword = *method* or *fmass* or *sp* or *temp* or *nhc* * keyword = *method* or *fmass* or *sp* or *temp* or *nhc*
@ -28,11 +28,15 @@ Examples
.. code-block:: LAMMPS .. 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 Description
""""""""""" """""""""""
.. versionchanged:: TBD
Fix pimd was renamed to fix pimd/nvt.
This command performs quantum molecular dynamics simulations based on This command performs quantum molecular dynamics simulations based on
the Feynman path integral to include effects of tunneling and the Feynman path integral to include effects of tunneling and
zero-point motion. In this formalism, the isomorphism of a quantum 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:: .. 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 combined with NH massive chain thermostat, so no other time
integration fix should be used. 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 normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics
(CMD). The difference between the styles is as follows. (CMD). The difference between the styles is as follows.
In standard PIMD, the value used for a bead's fictitious mass is 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 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 mass of the entire ring-polymer being equal to the real quantum
particle. But it can be difficult to efficiently integrate the particle. But it can be difficult to efficiently integrate the
equations of motion for the stiff harmonic interactions in the ring equations of motion for the stiff harmonic interactions in the ring
polymers. polymers.
A useful way to resolve this issue is to integrate the equations of A useful way to resolve this issue is to integrate the equations of
motion in a normal mode representation, using Normal Mode motion in a normal mode representation, using Normal Mode
Path-Integral Molecular Dynamics (NMPIMD) :ref:`(Cao1) <Cao1>`. In NMPIMD, Path-Integral Molecular Dynamics (NMPIMD) :ref:`(Cao1) <Cao1>`. In
the NH chains are attached to each normal mode of the ring-polymer and NMPIMD, the NH chains are attached to each normal mode of the
the fictitious mass of each mode is chosen as Mk = the eigenvalue of ring-polymer and the fictitious mass of each mode is chosen as Mk =
the Kth normal mode for k > 0. The k = 0 mode, referred to as the the eigenvalue of the Kth normal mode for k > 0. The k = 0 mode,
zero-frequency mode or centroid, corresponds to overall translation of referred to as the zero-frequency mode or centroid, corresponds to
the ring-polymer and is assigned the mass of the real particle. 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 Motion of the centroid can be effectively uncoupled from the other
normal modes by scaling the fictitious masses to achieve a partial normal modes by scaling the fictitious masses to achieve a partial
adiabatic separation. This is called a Centroid Molecular Dynamics adiabatic separation. This is called a Centroid Molecular Dynamics
(CMD) approximation :ref:`(Cao2) <Cao2>`. The time-evolution (and resulting (CMD) approximation :ref:`(Cao2) <Cao2>`. The time-evolution (and
dynamics) of the quantum particles can be used to obtain centroid time resulting dynamics) of the quantum particles can be used to obtain
correlation functions, which can be further used to obtain the true centroid time correlation functions, which can be further used to
quantum correlation function for the original system. The CMD method obtain the true quantum correlation function for the original system.
also uses normal modes to evolve the system, except only the k > 0 The CMD method also uses normal modes to evolve the system, except
modes are thermostatted, not the centroid degrees of freedom. only the k > 0 modes are thermostatted, not the centroid degrees of
freedom.
The keyword *fmass* sets a further scaling factor for the fictitious The keyword *fmass* sets a further scaling factor for the fictitious
masses of beads, which can be used for the Partial Adiabatic CMD 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 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 quasi-beads to :doc:`binary restart files <restart>`. See the
:doc:`read_restart <read_restart>` command for info on how to re-specify :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 a fix in an input script that reads a restart file, so that the
operation of the fix continues in an uninterrupted fashion. operation of the fix continues in an uninterrupted fashion.
None of the :doc:`fix_modify <fix_modify>` options 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 Fix pimd/nvt computes a global 3-vector, which can be accessed by
:doc:`output commands <Howto_output>`. The three quantities in the various :doc:`output commands <Howto_output>`. The three quantities in
global vector are the global vector are:
#. the total spring energy of the quasi-beads, #. the total spring energy of the quasi-beads,
#. the current temperature of the classical system of ring polymers, #. the current temperature of the classical system of ring polymers,
#. the current value of the scalar virial estimator for the kinetic #. the current value of the scalar virial estimator for the kinetic
energy of the quantum system :ref:`(Herman) <Herman>`. 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". temperature, which is "intensive".
No parameter of this fix can be used with the *start/stop* keywords of No parameter of fix pimd/nvt can be used with the *start/stop* keywords
the :doc:`run <run>` command. This fix is not invoked during of the :doc:`run <run>` command. Fix pimd/nvt is not invoked during
:doc:`energy minimization <minimize>`. :doc:`energy minimization <minimize>`.
Restrictions Restrictions
"""""""""""" """"""""""""
This fix is part of the REPLICA package. It is only enabled if This fix is part of the REPLICA package. It is only enabled if LAMMPS
LAMMPS was built with that package. See the was built with that package. See the :doc:`Build package
:doc:`Build package <Build_package>` page for more info. <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 A PIMD simulation can be initialized with a single data file read via
the :doc:`read_data <read_data>` command. However, this means all the :doc:`read_data <read_data>` command. However, this means all
quasi-beads in a ring polymer will have identical positions and quasi-beads in a ring polymer will have identical positions and
velocities, resulting in identical trajectories for all quasi-beads. velocities, resulting in identical trajectories for all quasi-beads. To
To avoid this, users can simply initialize velocities with different avoid this, users can simply initialize velocities with different random
random number seeds assigned to each partition, as defined by the number seeds assigned to each partition, as defined by the uloop
uloop variable, e.g. variable, e.g.
.. code-block:: LAMMPS .. code-block:: LAMMPS
@ -201,8 +207,8 @@ uloop variable, e.g.
Default Default
""""""" """""""
The keyword defaults are method = pimd, fmass = 1.0, sp = 1.0, temp = 300.0, The keyword defaults for fix pimd/nvt are method = pimd, fmass = 1.0, sp
and nhc = 2. = 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) print("Total number of style index entries:", total_index)
skip_angle = ('sdk') 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_pair = ('meam/c','lj/sf','reax/c','lj/sdk','lj/sdk/coul/long','lj/sdk/coul/msm')
skip_compute = ('pressure/cylinder') skip_compute = ('pressure/cylinder')
@ -286,7 +286,7 @@ if counter:
counter = 0 counter = 0
counter += check_style_index("compute", compute, index["compute"], skip=['pressure/cylinder']) 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("angle_style", angle, index["angle_style"], skip=['sdk'])
counter += check_style_index("bond_style", bond, index["bond_style"]) counter += check_style_index("bond_style", bond, index["bond_style"])
counter += check_style_index("dihedral_style", dihedral, index["dihedral_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 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_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 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 #read_restart system_${ibead}.rest1
special_bonds charmm 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 10
thermo_style custom step temp pe etotal f_1[1] f_1[2] f_1[3] 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_mol_swap.h
/fix_pimd.cpp /fix_pimd.cpp
/fix_pimd.h /fix_pimd.h
/fix_pimd_nvt.cpp
/fix_pimd_nvt.h
/fix_qbmsst.cpp /fix_qbmsst.cpp
/fix_qbmsst.h /fix_qbmsst.h
/fix_qtb.cpp /fix_qtb.cpp

View File

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

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Package FixPIMD Package FixPIMDNVT
Purpose Quantum Path Integral Algorithm for Quantum Chemistry Purpose Quantum Path Integral Algorithm for Quantum Chemistry
Copyright Voth Group @ University of Chicago Copyright Voth Group @ University of Chicago
Authors Chris Knight & Yuxing Peng (yuxing at uchicago.edu) Authors Chris Knight & Yuxing Peng (yuxing at uchicago.edu)
@ -21,7 +21,7 @@
Version 1.0 Version 1.0
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "fix_pimd.h" #include "fix_pimd_nvt.h"
#include "atom.h" #include "atom.h"
#include "comm.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; max_nsend = 0;
tag_send = nullptr; 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) else if (strcmp(arg[i + 1], "cmd") == 0)
method = CMD; method = CMD;
else else
error->universe_all(FLERR, fmt::format("Unknown method parameter {} for fix pimd", error->universe_all(
arg[i + 1])); FLERR, fmt::format("Unknown method parameter {} for fix pimd/nvt", arg[i + 1]));
} else if (strcmp(arg[i], "fmass") == 0) { } else if (strcmp(arg[i], "fmass") == 0) {
fmass = utils::numeric(FLERR, arg[i + 1], false, lmp); fmass = utils::numeric(FLERR, arg[i + 1], false, lmp);
if ((fmass < 0.0) || (fmass > np)) 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) { } else if (strcmp(arg[i], "sp") == 0) {
sp = utils::numeric(FLERR, arg[i + 1], false, lmp); 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) { } else if (strcmp(arg[i], "temp") == 0) {
nhc_temp = utils::numeric(FLERR, arg[i + 1], false, lmp); 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) { } else if (strcmp(arg[i], "nhc") == 0) {
nhc_nchain = utils::inumeric(FLERR, arg[i + 1], false, lmp); 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 } 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) 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 */ /* Initiation */
@ -138,7 +138,7 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixPIMD::~FixPIMD() FixPIMDNVT::~FixPIMDNVT()
{ {
delete[] mass; delete[] mass;
atom->delete_callback(id, Atom::GROW); atom->delete_callback(id, Atom::GROW);
@ -170,7 +170,7 @@ FixPIMD::~FixPIMD()
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int FixPIMD::setmask() int FixPIMDNVT::setmask()
{ {
int mask = 0; int mask = 0;
mask |= POST_FORCE; mask |= POST_FORCE;
@ -181,13 +181,13 @@ int FixPIMD::setmask()
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixPIMD::init() void FixPIMDNVT::init()
{ {
if (atom->map_style == Atom::MAP_NONE) 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) 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 // prepare the constants
@ -222,7 +222,7 @@ void FixPIMD::init()
fbond = -_fbond * force->mvv2e; fbond = -_fbond * force->mvv2e;
if (universe->me == 0) 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; dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v; 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) if (universe->me == 0 && universe->uscreen)
fprintf(universe->uscreen, "Setting up Path-Integral ...\n"); 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_v();
nhc_update_x(); nhc_update_x();
@ -259,14 +259,14 @@ void FixPIMD::initial_integrate(int /*vflag*/)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixPIMD::final_integrate() void FixPIMDNVT::final_integrate()
{ {
nhc_update_v(); 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 i = 0; i < atom->nlocal; i++)
for (int j = 0; j < 3; j++) atom->f[i][j] /= np; 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 Nose-Hoover Chains
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixPIMD::nhc_init() void FixPIMDNVT::nhc_init()
{ {
double tau = 1.0 / omega_np; double tau = 1.0 / omega_np;
double KT = force->boltz * nhc_temp; double KT = force->boltz * nhc_temp;
@ -309,7 +309,7 @@ void FixPIMD::nhc_init()
nhc_eta_dotdot[i][ichain] = 0.0; nhc_eta_dotdot[i][ichain] = 0.0;
nhc_eta_mass[i][ichain] = mass0; nhc_eta_mass[i][ichain] = mass0;
if ((method == CMD || method == NMPIMD) && universe->iworld == 0) if ((method == CMD || method == NMPIMD) && universe->iworld == 0)
; // do nothing ; // do nothing
else else
nhc_eta_mass[i][ichain] *= fmass; 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; int n = atom->nlocal;
double **x = atom->x; 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 n = atom->nlocal;
int *type = atom->type; int *type = atom->type;
@ -447,14 +447,14 @@ void FixPIMD::nhc_update_v()
Normal Mode PIMD Normal Mode PIMD
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixPIMD::nmpimd_init() void FixPIMDNVT::nmpimd_init()
{ {
memory->create(M_x2xp, np, np, "fix_feynman:M_x2xp"); memory->create(M_x2xp, np, np, "fix_feynman:M_x2xp");
memory->create(M_xp2x, np, np, "fix_feynman:M_xp2x"); memory->create(M_xp2x, np, np, "fix_feynman:M_xp2x");
memory->create(M_f2fp, np, np, "fix_feynman:M_f2fp"); memory->create(M_f2fp, np, np, "fix_feynman:M_f2fp");
memory->create(M_fp2f, np, np, "fix_feynman:M_fp2f"); 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 // 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_ptr = ptr;
comm->forward_comm(this); 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 n = atom->nlocal;
int m = 0; 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; spring_energy = 0.0;
@ -576,7 +576,7 @@ void FixPIMD::spring_force()
Comm operations Comm operations
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixPIMD::comm_init() void FixPIMDNVT::comm_init()
{ {
if (size_plan) { if (size_plan) {
delete[] plan_send; 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; int nlocal = atom->nlocal;
if (nlocal > max_nlocal) { if (nlocal > max_nlocal) {
max_nlocal = nlocal + 200; max_nlocal = nlocal + 200;
int size = sizeof(double) * max_nlocal * 3; 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++) 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 // copy local positions
@ -666,9 +666,9 @@ void FixPIMD::comm_exec(double **ptr)
if (nsend > max_nsend) { if (nsend > max_nsend) {
max_nsend = nsend + 200; max_nsend = nsend + 200;
tag_send = tag_send =
(tagint *) memory->srealloc(tag_send, sizeof(tagint) * max_nsend, "FixPIMD:tag_send"); (tagint *) memory->srealloc(tag_send, sizeof(tagint) * max_nsend, "FixPIMDNVT:tag_send");
buf_send = buf_send = (double *) memory->srealloc(buf_send, sizeof(double) * max_nsend * 3,
(double *) memory->srealloc(buf_send, sizeof(double) * max_nsend * 3, "FixPIMD:x_send"); "FixPIMDNVT:x_send");
} }
// send tags // 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; 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; int i, m, last;
@ -744,28 +744,28 @@ void FixPIMD::unpack_forward_comm(int n, int first, double *buf)
Memory operations Memory operations
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double FixPIMD::memory_usage() double FixPIMDNVT::memory_usage()
{ {
return (double) atom->nmax * size_peratom_cols * sizeof(double); 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; if (nmax == 0) return;
int count = nmax * 3; int count = nmax * 3;
memory->grow(array_atom, nmax, size_peratom_cols, "FixPIMD::array_atom"); memory->grow(array_atom, nmax, size_peratom_cols, "pimd_nvt:array_atom");
memory->grow(nhc_eta, count, nhc_nchain, "FixPIMD::nh_eta"); memory->grow(nhc_eta, count, nhc_nchain, "pimd_nvt:nh_eta");
memory->grow(nhc_eta_dot, count, nhc_nchain + 1, "FixPIMD::nh_eta_dot"); memory->grow(nhc_eta_dot, count, nhc_nchain + 1, "pimd_nvt:nh_eta_dot");
memory->grow(nhc_eta_dotdot, count, nhc_nchain, "FixPIMD::nh_eta_dotdot"); memory->grow(nhc_eta_dotdot, count, nhc_nchain, "pimd_nvt:nh_eta_dotdot");
memory->grow(nhc_eta_mass, count, nhc_nchain, "FixPIMD::nh_eta_mass"); 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 i_pos = i * 3;
int j_pos = j * 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 offset = 0;
int pos = i * 3; 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 offset = 0;
int pos = nlocal * 3; 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 offset = 0;
int pos = i * 3; 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; 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; return size_peratom_cols + 1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int FixPIMD::size_restart(int /*nlocal*/) int FixPIMDNVT::size_restart(int /*nlocal*/)
{ {
return size_peratom_cols + 1; 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 == 0) { return spring_energy; }
if (n == 1) { return t_sys; } if (n == 1) { return t_sys; }

View File

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