reformatted code and doc page
This commit is contained in:
Steve Plimpton
@ -1,106 +1,119 @@
|
||||
.. index:: fix hmc
|
||||
| .. index:: fix hmc
|
||||
|
||||
fix hmc command
|
||||
===============
|
||||
| fix hmc command
|
||||
| ===============
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
.. code-block:: LAMMPS
|
||||
| Syntax
|
||||
| """"""
|
||||
| .. code-block:: LAMMPS
|
||||
|
||||
fix ID group-ID hmc N seed temp integrator keyword values ...
|
||||
| fix ID group-ID hmc N seed temp integrator keyword values ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* hmc = style name of this fix command
|
||||
* N = invoke a Monto Carlo step every N steps
|
||||
* seed = random # seed (positive integer)
|
||||
* temp = temperature for assigning velocities
|
||||
* integrator = integrator fix: flexible (for nve) or rigid (for rigid/small)
|
||||
* keyword = *mom* or *resample*
|
||||
| * ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
| * hmc = style name of this fix command
|
||||
| * N = invoke a Monto Carlo step every N steps
|
||||
| * seed = random # seed (positive integer)
|
||||
| * temp = temperature for assigning velocities
|
||||
| * integrator = *flexible* (for nve) or *rigid* (for rigid/small)
|
||||
| * one or more keyword/value pairs may be appended
|
||||
|
||||
.. parsed-literal::
|
||||
| .. parsed-literal::
|
||||
|
||||
*mom* value = *yes* or *no*
|
||||
*resample* value = *yes* or *no*
|
||||
| keyword = *rigid* or *mom* or *resample*
|
||||
| *mom* value = *yes* or *no*
|
||||
| *resample* value = *yes* or *no*
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
| Examples
|
||||
| """"""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
| .. code-block:: LAMMPS
|
||||
|
||||
fix 1 all hmc 10 123 500 flexible
|
||||
fix hmc_water all hmc 100 123 298.15 rigid
|
||||
fix 2 all hmc 10 12345 300 flexible mom no resample yes
|
||||
| fix 1 all hmc 10 123 500 flexible
|
||||
| fix hmc_water all hmc 100 123 298.15 rigid
|
||||
| fix 2 all hmc 10 12345 300 flexible mom no resample yes
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
| Description
|
||||
| """""""""""
|
||||
|
||||
.. versionadded:: TBD
|
||||
| .. versionadded:: TBD
|
||||
|
||||
This fix performs the the Hybrid/Hamiltonian Monte Carlo (HMC) algorithm
|
||||
in line with the following order of steps:
|
||||
| This fix implements the so-called hybrid or Hamiltonian Monte Carlo
|
||||
| (HMC) algorithm. The basic idea is
|
||||
|
||||
The new particle configuration (positions and velocities) is calculated
|
||||
by invoking the velocity-Verlet time integration algorithm.
|
||||
Before these configuration changes are performed, the proposed change
|
||||
in the Hamiltonian, :math:`\Delta{H}` is calculated following the equation:
|
||||
|
||||
.. math::
|
||||
| in line with the following order of steps:
|
||||
|
||||
\Delta{H} = H(q',p') - H(q,p)
|
||||
| The new particle configuration (positions and velocities) is calculated
|
||||
| by invoking the velocity-Verlet time integration algorithm.
|
||||
| Before these configuration changes are performed, the proposed change
|
||||
| in the Hamiltonian, :math:`\Delta{H}` is calculated following the equation:
|
||||
|
||||
This new proposed configuration is then accepted/rejected according to
|
||||
the Metropolis criterion with probability:
|
||||
| .. math::
|
||||
|
||||
.. math::
|
||||
| \Delta{H} = H(q',p') - H(q,p)
|
||||
|
||||
p^{acc} = min(1,e^{\frac{-\Delta{H}}{k_B T}})
|
||||
| This new proposed configuration is then accepted/rejected according to
|
||||
| the Metropolis criterion with probability:
|
||||
|
||||
Upon acceptance, the new proposed particle configuration positions and
|
||||
velocities are updated. Upon rejection, the old particle configuration
|
||||
is kept, and particle momenta (and therefore velocities) are randomly
|
||||
resampled from a normal distribution:
|
||||
| .. math::
|
||||
|
||||
.. math::
|
||||
| p^{acc} = min(1,e^{\frac{-\Delta{H}}{k_B T}})
|
||||
|
||||
p_{x,y,z} = \textbf{N}(0,1) \sqrt{\frac{k_B T}{2 m^2}}
|
||||
| Upon acceptance, the new proposed particle configuration positions and
|
||||
| velocities are updated. Upon rejection, the old particle configuration
|
||||
| is kept, and particle momenta (and therefore velocities) are randomly
|
||||
| resampled from a normal distribution:
|
||||
|
||||
The algorithm then continues, proposing a new configuration of particles
|
||||
and velocities N integration steps later.
|
||||
| .. math::
|
||||
|
||||
Typically, HMC is run with a larger timestep than molecular dynamics.
|
||||
The timestep provides control over the acceptance ratio. A larger
|
||||
timestep will lead to larger and more extreme MC moves which are
|
||||
less likely to be accepted.
|
||||
| p_{x,y,z} = \textbf{N}(0,1) \sqrt{\frac{k_B T}{2 m^2}}
|
||||
|
||||
This fix is not designed to be used with anything but an NVE simulation.
|
||||
Only atom-specific data are restored on MC move rejection, so anything which
|
||||
adds or remove atoms, changes the box size, or has some external state not
|
||||
dependent on atomic data will have undefined behavior.
|
||||
| The algorithm then continues, proposing a new configuration of particles
|
||||
| and velocities N integration steps later.
|
||||
|
||||
The group assigned to this fix has no meaning. All
|
||||
| Typically, HMC is run with a larger timestep than would be used for
|
||||
| traditional molecular dynamics (MD). The timestep provides control
|
||||
| over the acceptance ratio. A larger timestep will lead to larger and
|
||||
| more extreme MC moves which are less likely to be accepted.
|
||||
|
||||
----------
|
||||
| This fix is not designed to be used with anything but an NVE
|
||||
| simulation. Only atom-specific data are restored on MC move
|
||||
| rejection, so anything which adds or remove atoms, changes the box
|
||||
| size, or has some external state not dependent on atomic data will
|
||||
| have undefined behavior.
|
||||
|
||||
The keyword/value options are used in the following ways.
|
||||
| The group assigned to this fix has no meaning is ignored.
|
||||
|
||||
The *mom* keyword sets the linear momentum of the ensemble of particles.
|
||||
If mom = yes, the linear momentum of the ensemble of velocities is
|
||||
zeroed. If mom = no, the linear momentum of the ensemble of velocities
|
||||
is not zeroed.
|
||||
| NOTE: mention adds several new computes
|
||||
| when how do those computes get triggered ?
|
||||
| why does global virial need to be stored with state ?
|
||||
| is it simply for output at end of step
|
||||
| will thermo output at end-of-step be correct ?
|
||||
| no references in the text
|
||||
|
||||
The *resample* keyword decides whether velocities are resampled upon acceptance.
|
||||
If resample = yes, velocities are resampled upon acceptance. If resample = no,
|
||||
velocities are not resampled upon acceptance.
|
||||
| ----------
|
||||
|
||||
| The keyword/value options are used in the following ways.
|
||||
|
||||
| The *mom* keyword sets the linear momentum of the ensemble of
|
||||
particles. If mom = yes, the linear momentum of the ensemble of
|
||||
velocities is zeroed. If mom = no, the linear momentum of the ensemble
|
||||
of velocities is not zeroed.
|
||||
|
||||
The *resample* keyword decides whether velocities are resampled upon
|
||||
acceptance. If resample = yes, velocities are resampled upon
|
||||
acceptance. If resample = no, velocities are not resampled upon
|
||||
acceptance.
|
||||
|
||||
----------
|
||||
|
||||
Ouput info
|
||||
""""""""""
|
||||
|
||||
This fix computes a global scalar and global vector of length 5,
|
||||
which can be accessed by various :doc:`output commands
|
||||
<Howto_output>`. The scalar is the fraction of attempted MC moves which have been accepted. The vector stores the
|
||||
following quantities:
|
||||
This fix computes a global scalar and global vector of length 5, which
|
||||
can be accessed by various :doc:`output commands <Howto_output>`. The
|
||||
scalar is the fraction of attempted MC moves which have been accepted.
|
||||
The vector stores the following quantities:
|
||||
|
||||
* 1 = number of accepted moves
|
||||
* 2 = number of rejected moves
|
||||
@ -108,31 +121,37 @@ following quantities:
|
||||
* 4 = change in kinetic energy
|
||||
* 5 = change in total energy (kinetic + potential energy)
|
||||
|
||||
These values are calculated every N timesteps
|
||||
These values are calculated once every N timesteps
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This fix is part of the MC package and requires the RIGID package to
|
||||
be installed. It is only enabled if LAMMPS was built with both packages.
|
||||
See the :doc:`Build package <Build_package>` doc page for more info.
|
||||
be installed. It is only enabled if LAMMPS was built with both
|
||||
packages. See the :doc:`Build package <Build_package>` doc page for
|
||||
more info.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix nvt <fix_nh>`, :doc:`fix gcmc <fix_gcmc>`, :doc:`fix tfmc <fix_tfmc>`
|
||||
:doc:`fix nvt <fix_nh>`, :doc:`fix gcmc <fix_gcmc>`, :doc:`fix tfmc
|
||||
<fix_tfmc>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
The option default is mom = yes, resample = no.
|
||||
The option defaults are mom = yes, resample = no.
|
||||
|
||||
----------
|
||||
|
||||
**(Watkins)** Watkins and Jorgensen, J Phys Chem A, 105, 4118-4125 (2001).
|
||||
**(Watkins)** Watkins and Jorgensen, J Phys Chem A, 105, 4118-4125
|
||||
(2001).
|
||||
|
||||
**(Betancourt)** Betancourt, M. A Conceptual Introduction to Hamiltonian Monte Carlo, 2018.
|
||||
**(Betancourt)** Betancourt, Conceptual Introduction to Hamiltonian
|
||||
Monte Carlo, 2018.
|
||||
|
||||
**(Duane)** Duane, S.; Kennedy, A. D.; Pendleton, B. J.; Roweth, D. Hybrid Monte Carlo. Physics Letters B 1987, 195 (2), 216-222. https://doi.org/10.1016/0370-2693(87)91197-X.
|
||||
**(Duane)** Duane, Kennedy, Pendleton, and Roweth, Physics Letters B,
|
||||
195 (2), 216-222 (1987). https://doi.org/10.1016/0370-2693(87)91197-X
|
||||
|
||||
**(Metropolis)** Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. The journal of chemical physics 1953, 21, 1087-1092.
|
||||
**(Metropolis)** Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller,
|
||||
and E. Teller, J Chemical Physics, 21, 1087-1092 (1953).
|
||||
|
||||
@ -4,7 +4,7 @@ units real
|
||||
boundary p p p
|
||||
atom_style full
|
||||
|
||||
read_data data.hmc_spce
|
||||
read_data data.hmc.spce
|
||||
replicate 4 4 4
|
||||
|
||||
mass 1 15.9994
|
||||
@ -65,27 +65,27 @@ FixHMC::FixHMC(LAMMPS *lmp, int narg, char **arg) :
|
||||
eglobal(nullptr), eglobalptr(nullptr), vglobal(nullptr), vglobalptr(nullptr), pe(nullptr),
|
||||
ke(nullptr), peatom(nullptr), press(nullptr), pressatom(nullptr), buf_store(nullptr)
|
||||
{
|
||||
// set some defaults
|
||||
// defaults
|
||||
|
||||
mom_flag = 1;
|
||||
resample_on_accept_flag = 0;
|
||||
rigid_flag = false;
|
||||
if (narg < 5) utils::missing_cmd_args(FLERR, "fix hmc", error);
|
||||
|
||||
// Retrieve user-defined options:
|
||||
// required arguments
|
||||
|
||||
nevery = utils::numeric(FLERR, arg[3], false, lmp); // Number of MD steps per MC step
|
||||
int seed = utils::numeric(FLERR, arg[4], false, lmp); // Seed for random number generation
|
||||
double temp = utils::numeric(FLERR, arg[5], false, lmp); // System temperature
|
||||
if (seed <= 0)
|
||||
error->all(FLERR, "Illegal fix hmc seed argument: {}. Seed must be greater than 0.0", seed);
|
||||
if (temp <= 0)
|
||||
error->all(FLERR, "Illegal fix hmc temp argument: {}. Temp must be greater than 0.0", temp);
|
||||
nevery = utils::numeric(FLERR, arg[3], false, lmp);
|
||||
int seed = utils::numeric(FLERR, arg[4], false, lmp);
|
||||
double temp = utils::numeric(FLERR, arg[5], false, lmp);
|
||||
|
||||
if (seed <= 0) error->all(FLERR, "Fix hmc seed must be > 0");
|
||||
if (temp <= 0) error->all(FLERR, "Fix hmc temperature must be > 0.0");
|
||||
|
||||
KT = force->boltz * temp / force->mvv2e; // K*T in mvv units
|
||||
mbeta = -1.0 / (force->boltz * temp); // -1/(K*T) in energy units
|
||||
|
||||
// Check optional keywords:
|
||||
// optional keywords
|
||||
|
||||
mom_flag = 1;
|
||||
resample_on_accept_flag = 0;
|
||||
rigid_flag = false;
|
||||
|
||||
int iarg = 6;
|
||||
while (iarg < narg) {
|
||||
@ -112,15 +112,17 @@ FixHMC::FixHMC(LAMMPS *lmp, int narg, char **arg) :
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize RNG with a different seed for each process:
|
||||
// initialize RNG with a different seed for each process
|
||||
|
||||
random = new RanPark(lmp, seed + comm->me);
|
||||
for (int i = 0; i < 100; i++) random->uniform();
|
||||
random_equal = new RanPark(lmp, seed);
|
||||
|
||||
// Register callback:
|
||||
// register callback
|
||||
|
||||
atom->add_callback(0);
|
||||
|
||||
// Add new computes for global and per-atom properties:
|
||||
// add new computes for global and per-atom properties
|
||||
|
||||
ke = modify->add_compute(fmt::format("hmc_ke_{} all ke", id));
|
||||
pe = modify->add_compute(fmt::format("hmc_pe_{} all pe", id));
|
||||
@ -128,7 +130,7 @@ FixHMC::FixHMC(LAMMPS *lmp, int narg, char **arg) :
|
||||
press = modify->add_compute(fmt::format("hmc_press_{} all pressure NULL virial", id));
|
||||
pressatom = modify->add_compute(fmt::format("hmc_pressatom_{} all stress/atom NULL virial", id));
|
||||
|
||||
// Define non-default fix attributes:
|
||||
// fix attributes
|
||||
|
||||
global_freq = 1;
|
||||
scalar_flag = 1;
|
||||
@ -138,10 +140,12 @@ FixHMC::FixHMC(LAMMPS *lmp, int narg, char **arg) :
|
||||
size_vector = 5;
|
||||
force_reneighbor = 1;
|
||||
next_reneighbor = -1;
|
||||
|
||||
first_init_complete = false;
|
||||
first_setup_complete = false;
|
||||
|
||||
// initialize quantities
|
||||
// initializations
|
||||
|
||||
nattempts = 0;
|
||||
naccepts = 0;
|
||||
DeltaPE = 0;
|
||||
@ -156,12 +160,15 @@ FixHMC::~FixHMC()
|
||||
|
||||
memory->destroy(eglobal);
|
||||
memory->destroy(vglobal);
|
||||
|
||||
delete[] eglobalptr;
|
||||
if (vglobalptr)
|
||||
for (int m = 0; m < nv; m++) delete[] vglobalptr[m];
|
||||
delete[] vglobalptr;
|
||||
|
||||
delete random;
|
||||
delete random_equal;
|
||||
|
||||
modify->delete_compute(std::string("hmc_ke_") + id);
|
||||
modify->delete_compute(std::string("hmc_pe_") + id);
|
||||
modify->delete_compute(std::string("hmc_peatom_") + id);
|
||||
@ -181,7 +188,8 @@ void FixHMC::setup_arrays_and_pointers()
|
||||
int improper_flag;
|
||||
int kspace_flag;
|
||||
|
||||
// Determine which energy contributions must be computed:
|
||||
// determine which energy contributions must be computed
|
||||
|
||||
ne = 0;
|
||||
if (force->pair) {
|
||||
pair_flag = 1;
|
||||
@ -214,23 +222,27 @@ void FixHMC::setup_arrays_and_pointers()
|
||||
} else
|
||||
kspace_flag = 0;
|
||||
|
||||
// Initialize arrays for managing global energy terms:
|
||||
// initialize arrays for managing global energy terms
|
||||
|
||||
neg = pair_flag ? ne + 1 : ne;
|
||||
memory->create(eglobal, neg, "fix_hmc:eglobal");
|
||||
delete[] eglobalptr;
|
||||
eglobalptr = new double *[neg];
|
||||
|
||||
m = 0;
|
||||
if (pair_flag) {
|
||||
eglobalptr[m++] = &force->pair->eng_vdwl;
|
||||
eglobalptr[m++] = &force->pair->eng_coul;
|
||||
}
|
||||
|
||||
if (bond_flag) eglobalptr[m++] = &force->bond->energy;
|
||||
if (angle_flag) eglobalptr[m++] = &force->angle->energy;
|
||||
if (dihedral_flag) eglobalptr[m++] = &force->dihedral->energy;
|
||||
if (improper_flag) eglobalptr[m++] = &force->improper->energy;
|
||||
if (kspace_flag) eglobalptr[m++] = &force->kspace->energy;
|
||||
|
||||
// Initialize arrays for managing global virial terms:
|
||||
// initialize arrays for managing global virial terms
|
||||
|
||||
nv = ne;
|
||||
for (const auto &ifix : modify->get_fix_list())
|
||||
if (ifix->virial_global_flag) nv++;
|
||||
@ -239,7 +251,9 @@ void FixHMC::setup_arrays_and_pointers()
|
||||
for (m = 0; m < nv; m++) delete[] vglobalptr[m];
|
||||
delete[] vglobalptr;
|
||||
vglobalptr = new double **[nv];
|
||||
|
||||
for (m = 0; m < nv; m++) vglobalptr[m] = new double *[6];
|
||||
|
||||
for (i = 0; i < 6; i++) {
|
||||
m = 0;
|
||||
if (pair_flag) vglobalptr[m++][i] = &force->pair->virial[i];
|
||||
@ -252,11 +266,10 @@ void FixHMC::setup_arrays_and_pointers()
|
||||
if (ifix->virial_global_flag) vglobalptr[m++][i] = &ifix->virial[i];
|
||||
}
|
||||
|
||||
// Determine maximum number of per-atom variables in forward
|
||||
// communications when dealing with rigid bodies:
|
||||
if (rigid_flag) {
|
||||
comm_forward = 12;
|
||||
}
|
||||
// set maximum number of per-atom variables in forward comm
|
||||
// when dealing with rigid bodies
|
||||
|
||||
if (rigid_flag) comm_forward = 12;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -274,43 +287,50 @@ void FixHMC::init()
|
||||
{
|
||||
int ntimestep = update->ntimestep;
|
||||
|
||||
if (atom->tag_enable == 0) error->all(FLERR, "Cannot use fix hmc unless atoms have IDs");
|
||||
if (atom->tag_enable == 0)
|
||||
error->all(FLERR, "Cannot use fix hmc unless atoms have IDs");
|
||||
|
||||
// Check whether there is any fixes that change box size and/or shape:
|
||||
// check whether there is any fixes that change box size and/or shape
|
||||
|
||||
for (const auto &ifix : modify->get_fix_list())
|
||||
if (ifix->box_change)
|
||||
error->all(FLERR, "fix hmc is incompatible with fixes that change box size or shape");
|
||||
error->all(FLERR, "Fix hmc is incompatible with fixes that change box size or shape");
|
||||
|
||||
// Check whether there are subsequent fixes with active virial_flag:
|
||||
// check whether there are subsequent fixes with active virial_flag
|
||||
|
||||
int past_this_fix = false;
|
||||
int past_rigid = !rigid_flag;
|
||||
for (const auto &ifix : modify->get_fix_list()) {
|
||||
if (past_this_fix && past_rigid && ifix->virial_global_flag) {
|
||||
if (comm->me == 0) utils::logmesg(lmp, "Fix {} defined after fix hmc.\n", ifix->style);
|
||||
error->all(FLERR, "fix hmc cannot precede fixes that modify the system pressure");
|
||||
error->all(FLERR, "Fix hmc cannot precede fixes that contribute to the system pressure");
|
||||
}
|
||||
if (!strcmp(ifix->id, id)) past_this_fix = true;
|
||||
if (rigid_flag && !strcmp(ifix->id, fix_rigid->id)) past_rigid = true;
|
||||
}
|
||||
|
||||
if (!first_init_complete)
|
||||
{
|
||||
// Look for computes with active press_flag:
|
||||
if (!first_init_complete) {
|
||||
|
||||
// look for computes with active press_flag
|
||||
|
||||
press_flag = 0;
|
||||
for (const auto &icompute : modify->get_compute_list())
|
||||
if (strncmp(icompute->id, "hmc_", 4)) {
|
||||
press_flag = press_flag | icompute->pressflag;
|
||||
}
|
||||
|
||||
// Initialize arrays and pointers for saving/restoration of states:
|
||||
// initialize arrays and pointers for saving/restoring state
|
||||
|
||||
setup_arrays_and_pointers();
|
||||
|
||||
// Count per-atom properties to be exchanged:
|
||||
// count per-atom properties to be exchanged
|
||||
|
||||
nvalues = 0;
|
||||
if (peatom_flag) nvalues += ne;
|
||||
if (pressatom_flag) nvalues += 6 * nv;
|
||||
|
||||
// Activate potential energy and other necessary calculations at setup:
|
||||
// activate potential energy and other necessary calculations at setup
|
||||
|
||||
pe->addstep(ntimestep);
|
||||
if (peatom_flag) peatom->addstep(ntimestep);
|
||||
if (press_flag) press->addstep(ntimestep);
|
||||
@ -327,10 +347,12 @@ void FixHMC::init()
|
||||
void FixHMC::setup(int vflag)
|
||||
{
|
||||
if (!first_setup_complete) {
|
||||
// initialize rigid first to avoid saving uninitialized state
|
||||
|
||||
// initialize fix rigid first to avoid saving uninitialized state
|
||||
|
||||
if (rigid_flag) fix_rigid->setup(vflag);
|
||||
|
||||
// Compute properties of the initial state:
|
||||
// compute properties of initial state
|
||||
|
||||
if (rigid_flag) {
|
||||
rigid_body_random_velocities();
|
||||
@ -348,7 +370,8 @@ void FixHMC::setup(int vflag)
|
||||
memory->create(buf_store, maxstore + bufextra, "fix_hmc:buf_store");
|
||||
save_current_state();
|
||||
|
||||
// Activate potential energy and other necessary calculations:
|
||||
// activate potential energy and other necessary calculations
|
||||
|
||||
int nextstep = update->ntimestep + nevery;
|
||||
pe->addstep(nextstep);
|
||||
if (peatom_flag) peatom->addstep(nextstep);
|
||||
@ -356,6 +379,7 @@ void FixHMC::setup(int vflag)
|
||||
if (pressatom_flag) pressatom->addstep(nextstep);
|
||||
|
||||
// create buffer, store fixes for warnings
|
||||
|
||||
int maxexchange_fix = 0;
|
||||
int maxexchange_atom = atom->avec->maxexchange;
|
||||
|
||||
@ -366,41 +390,47 @@ void FixHMC::setup(int vflag)
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Apply the Metropolis acceptance criterion
|
||||
Restore saved system state if move is rejected
|
||||
Activate computes for the next MC step
|
||||
apply the Metropolis acceptance criterion
|
||||
restore saved system state if move is rejected
|
||||
activate computes for the next MC step
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixHMC::end_of_step()
|
||||
{
|
||||
nattempts++;
|
||||
|
||||
// Compute potential and kinetic energy variations:
|
||||
// compute potential and kinetic energy variations
|
||||
|
||||
update->eflag_global = update->ntimestep;
|
||||
double newPE = pe->compute_scalar();
|
||||
double newKE = ke->compute_scalar();
|
||||
DeltaPE = newPE - PE;
|
||||
DeltaKE = newKE - KE;
|
||||
|
||||
// Apply the Metropolis criterion:
|
||||
// apply the Metropolis criterion
|
||||
|
||||
double DeltaE = DeltaPE + DeltaKE;
|
||||
int accept = DeltaE < 0.0;
|
||||
if (!accept) {
|
||||
accept = random_equal->uniform() <= exp(mbeta * DeltaE);
|
||||
MPI_Bcast(&accept, 1, MPI_INT, 0, world);
|
||||
}
|
||||
|
||||
// for accept: update potential energy and save current state
|
||||
// for reject: restore saved state, exchange, and
|
||||
// enforce check_distance/reneighboring on next step
|
||||
|
||||
if (accept) {
|
||||
// Update potential energy and save the current state:
|
||||
naccepts++;
|
||||
PE = newPE;
|
||||
save_current_state();
|
||||
} else {
|
||||
// Restore saved state, exchange, and enforce check_distance/reneighboring in the next step:
|
||||
restore_saved_state();
|
||||
next_reneighbor = update->ntimestep + 1;
|
||||
}
|
||||
|
||||
// Choose new velocities and compute kinetic energy:
|
||||
// choose new velocities and compute kinetic energy
|
||||
|
||||
if (!accept || resample_on_accept_flag) {
|
||||
if (rigid_flag)
|
||||
rigid_body_random_velocities();
|
||||
@ -409,7 +439,8 @@ void FixHMC::end_of_step()
|
||||
KE = ke->compute_scalar();
|
||||
}
|
||||
|
||||
// Activate potential energy and other necessary calculations:
|
||||
// activate potential energy and other necessary calculations
|
||||
|
||||
int nextstep = update->ntimestep + nevery;
|
||||
if (nextstep <= update->laststep) {
|
||||
pe->addstep(nextstep);
|
||||
@ -420,7 +451,7 @@ void FixHMC::end_of_step()
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Return the acceptance fraction of proposed MC moves
|
||||
return acceptance fraction of proposed MC moves
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixHMC::compute_scalar()
|
||||
@ -431,7 +462,7 @@ double FixHMC::compute_scalar()
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Return the acceptance fraction of proposed MC moves, or
|
||||
return the acceptance fraction of proposed MC moves, or
|
||||
the total energy variation of the last proposed MC move, or
|
||||
the mean-square atom displacement in the last proposed MC move
|
||||
------------------------------------------------------------------------- */
|
||||
@ -476,8 +507,9 @@ void FixHMC::grow_store(int n, int flag)
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Save the system state for eventual restoration if move is rejected
|
||||
save current system state for eventual restoration if move is rejected
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixHMC::save_current_state()
|
||||
{
|
||||
int m;
|
||||
@ -488,21 +520,25 @@ void FixHMC::save_current_state()
|
||||
AtomVec *avec = atom->avec;
|
||||
nstore = 0;
|
||||
|
||||
// store all needed info about owned atoms via pack_exchange()
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (nstore > maxstore) grow_store(nstore, 1);
|
||||
nstore += avec->pack_exchange(i, &buf_store[nstore]);
|
||||
if (nstore > maxstore) grow_store(nstore, 1);
|
||||
nstore += avec->pack_exchange(i, &buf_store[nstore]);
|
||||
}
|
||||
|
||||
// Save global energy terms:
|
||||
// save global energy terms
|
||||
|
||||
for (m = 0; m < neg; m++) eglobal[m] = *eglobalptr[m];
|
||||
|
||||
// Save global virial terms:
|
||||
// save global virial terms
|
||||
|
||||
if (press_flag)
|
||||
for (m = 0; m < nv; m++) memcpy(vglobal[m], *vglobalptr[m], six);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Restore the system state saved at the beginning of the MC step
|
||||
restore system state saved at the beginning of the MC step
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixHMC::restore_saved_state()
|
||||
@ -510,34 +546,45 @@ void FixHMC::restore_saved_state()
|
||||
int i;
|
||||
AtomVec *avec = atom->avec;
|
||||
|
||||
// clear atom map
|
||||
|
||||
if (atom->map_style != Atom::MAP_NONE) atom->map_clear();
|
||||
|
||||
// delete all owned + ghost atoms
|
||||
|
||||
atom->nghost = 0;
|
||||
atom->avec->clear_bonus();
|
||||
// delete all atoms
|
||||
atom->nlocal = 0;
|
||||
|
||||
// unpack exchange
|
||||
|
||||
int m = 0;
|
||||
while (m < nstore) {
|
||||
m += avec->unpack_exchange(&buf_store[m]);
|
||||
m += avec->unpack_exchange(&buf_store[m]);
|
||||
}
|
||||
|
||||
// NOTE: still no ghost atoms ?
|
||||
|
||||
// reinit atom_map
|
||||
|
||||
if (atom->map_style != Atom::MAP_NONE) {
|
||||
atom->map_clear();
|
||||
atom->map_init();
|
||||
atom->map_set();
|
||||
}
|
||||
|
||||
// Restore global energy terms:
|
||||
// restore global energy terms
|
||||
|
||||
for (i = 0; i < neg; i++) *eglobalptr[i] = eglobal[i];
|
||||
|
||||
// Restore global virial terms:
|
||||
// restore global virial terms
|
||||
|
||||
if (press_flag)
|
||||
for (i = 0; i < nv; i++) memcpy(*vglobalptr[i], vglobal[i], six);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Randomly choose velocities from a Maxwell-Boltzmann distribution
|
||||
randomly choose velocities from a Maxwell-Boltzmann distribution
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixHMC::random_velocities()
|
||||
@ -548,12 +595,12 @@ void FixHMC::random_velocities()
|
||||
|
||||
double stdev;
|
||||
int nlocal, dimension;
|
||||
//
|
||||
double *rmass = atom->rmass;
|
||||
double *mass = atom->mass;
|
||||
nlocal = atom->nlocal;
|
||||
dimension = domain->dimension;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass)
|
||||
@ -562,6 +609,7 @@ void FixHMC::random_velocities()
|
||||
stdev = sqrt(KT / mass[type[i]]);
|
||||
for (int j = 0; j < dimension; j++) v[i][j] = stdev * random->gaussian();
|
||||
}
|
||||
|
||||
if (mom_flag) {
|
||||
double vcm[3];
|
||||
group->vcm(igroup, group->mass(igroup), vcm);
|
||||
|
||||
@ -24,8 +24,6 @@
|
||||
|
||||
#include <cstring>
|
||||
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
static constexpr int DELTA = 10000;
|
||||
@ -91,7 +89,7 @@ ComputeRigidLocal::ComputeRigidLocal(LAMMPS *lmp, int narg, char **arg) :
|
||||
if (nvalues == 1) size_local_cols = 0;
|
||||
else size_local_cols = nvalues;
|
||||
|
||||
ncount = nmax = 0;
|
||||
ncount = nmax = 0;
|
||||
vlocal = nullptr;
|
||||
alocal = nullptr;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user