jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

more updates to hyper docs

This commit is contained in:
Steve Plimpton 2018-10-11 16:32:30 -06:00
parent 4a6f088c0b
commit 0c8ce199af
2 changed files with 245 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

77
doc/src/fix_hyper_global.txt
View File

@ -47,12 +47,12 @@ compared to running for the same number of timesteps with normal MD.
See the "hyper"_hyper.html doc page for a more general discussion of
hyperdynamics and citations that explain both GHD and LHD.
The equations and logic used by this fix to perform GHD follow the
description given in "(Voter2013)"_#Voter2013ghd. The bond-boost form
of a bias potential for HD is due to Miron and Fichthorn as described
in "(Miron)"_#Mironghd. In LAMMPS we use a simplified version of
bond-boost GHD where a single bond in the system is biased at any one
timestep.
The equations and logic used by this fix and described here to perform
GHD follow the description given in "(Voter2013)"_#Voter2013ghd. The
bond-boost form of a bias potential for HD is due to Miron and
Fichthorn as described in "(Miron)"_#Mironghd. In LAMMPS we use a
simplified version of bond-boost GHD where a single bond in the system
is biased at any one timestep.
Bonds are defined between each pair of I,J atoms whose R0ij distance
is less than {cutbond}, when the system is in a quenched state
@ -87,7 +87,7 @@ that Vij(max) is added to the normal interatomic potential that is
computed between all atoms in the system at every step.
The derivative of Vij(max) with respect to the position of each atom
in the Emax bond gives a force Fij(max) acting on the bond as
in the Emax bond gives a bias force Fij(max) acting on the bond as
Fij(max) = - dVij(max)/dEij = 2 Vmax Eij / qfactor^2 for abs(Eij) < qfactor
= 0 otherwise :pre
@ -100,7 +100,7 @@ The time boost factor for the system is given each timestep I by
Bi = exp(beta * Vij(max)) :pre
where beta = 1/kTequil, and {Tequil} is the temperature of the system
and an arguement to this fix. Note that Bi >= 1 at every step.
and an argument to this fix. Note that Bi >= 1 at every step.
NOTE: To run GHD, the input script must also use the "fix
langevin"_fix_langevin.html command to thermostat the atoms at the
@ -108,21 +108,21 @@ same {Tequil} as specified by this fix, so that the system is running
constant-temperature (NVT) dynamics. LAMMPS does not check that this
is done.
The elapsed system time t_hyper for a GHD simulation running for {N}
The elapsed time t_hyper for a GHD simulation running for {N}
timesteps is simply
t_hyper = Sum (i = 1 to N) Bi * dt :pre
where dt is the timestep size defined by the "timestep"_timestep.html
command. The effective time acceleration due to GHD is thus t_hyper /
N*dt, where the denominator is elapsed time for a normal MD run.
N*dt, where N*dt is elapsed time for a normal MD run of N timesteps.
Note that in GHD, the boost factor varies from timestep to timestep.
Likewise, which bond has Emax strain and thus which pair of atoms the
bias potential is added to, will also vary from timestep to timestep.
This in contrast to local hyperdynamics (LHD) where the boost factor
is an input parameter; see the "fix hyper/local"_fix_hyper_local.html
doc page for details.
This is in contrast to local hyperdynamics (LHD) where the boost
factor is an input parameter; see the "fix
hyper/local"_fix_hyper_local.html doc page for details.
:line
@ -143,22 +143,32 @@ equilibrium (quenched) R0ij distance and the two atoms in the bond
could still experience a non-zero bias force.
If {qfactor} is set too large, then transitions from one energy basin
to another are affected because the bias potneial is non-zero at the
to another are affected because the bias potential is non-zero at the
transition state (e.g. saddle point). If {qfactor} is set too small
than little boost is achieved because the Eij strain of some bond in
the system will (nearly) always exceed {qfactor). A value of 0.3 for
{qfactor} is typically a reasonable value.
The {Vmax} argument is the prefactor on the bias potential. It should
be set to a value less than the smallest barrier height for an event
to occur. Otherwise the applied bias potential may be large enough
(when added to the interatomic potential) to produce a local energy
basin with a maxima in the center. This can produce artificial energy
minima in the same basin that trap an atom. Or if {Vmax} is even
larger, it may induce an atom(s) to rapidly transition to another
energy basin. Both cases are "bad dynamics" which violate the
assumptions of GHD that guarantee an accelerated time-accurate
trajectory of the system.
The {Vmax} argument is the prefactor on the bias potential. Ideally,
tt should be set to a value slightly less than the smallest barrier
height for an event to occur. Otherwise the applied bias potential
may be large enough (when added to the interatomic potential) to
produce a local energy basin with a maxima in the center. This can
produce artificial energy minima in the same basin that trap an atom.
Or if {Vmax} is even larger, it may induce an atom(s) to rapidly
transition to another energy basin. Both cases are "bad dynamics"
which violate the assumptions of GHD that guarantee an accelerated
time-accurate trajectory of the system.
Note that if {Vmax} is set too small, the GHD simulation will run
correctly. There will just be fewer events because the hyper time
(t_hyper equation above) will be shorter.
NOTE: If you have no physical intuition as to the smallest barrier
height in your system, a reasonable strategy to determine the largest
{Vmax} you can use for an LHD model, is to run a sequence of
simulations with smaller and smaller {Vmax} values, until the event
rate does not change.
The {Tequil} argument is the temperature at which the system is
simulated; see the comment above about the "fix
@ -167,7 +177,7 @@ beta term in the exponential factor that determines how much boost is
achieved as a function of the bias potential.
In general, the lower the value of {Tequil} and the higher the value
of {Vmax}, the more boost will be achieved by the GHD algorithm.
of {Vmax}, the more boost will be achievable by the GHD algorithm.
:line
@ -204,18 +214,23 @@ The first 5 quantities are for the current timestep. Quantities 6-8
are for the current hyper run. Quantities 9-11 are cummulative across
multiple runs (since the fix was defined in the input script).
For value 7, drift is the distance an atom moves between timesteps
when the bond list is reset, i.e. between events. Atoms involved in
an event will typically move the greatest distance since others are
typically oscillating around their lattice site.
For value 10, events are checked for by the "hyper"_hyper.html command
once every {Nevent} timesteps. This value is the count of those
timesteps on which one (or more) events was detected. It is NOT the
number of distinct events, since more than one physical event may
occur in the same {Nevent} time window.
number of distinct events, since more than one event may occur in the
same {Nevent} time window.
For value 11, each time the "hyper"_hyper.html command checks for an
event, it invokes a compute to flag zero or more atoms as
participating in an event. E.g. atoms that have displaced more than
some distance from the previous quench state. Value 11 is the
cummulative count of the number of atoms participating in any of the
events that were found.
participating in one or more events. E.g. atoms that have displaced
more than some distance from the previous quench state. Value 11 is
the cummulative count of the number of atoms participating in any of
the events that were found.
The scalar and vector values calculated by this fix are all
"intensive".

268
doc/src/fix_hyper_local.txt
View File

@ -10,7 +10,7 @@ fix hyper/local command :h3
[Syntax:]
fix ID group-ID hyper/local cutbond qfactor Vmax Tequil Dcut alpha boost :pre
fix ID group-ID hyper/local cutbond qfactor Vmax Tequil Dcut alpha Btarget :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
hyper/local = style name of this fix command :l
@ -20,14 +20,9 @@ Vmax = estimated height of bias potential (energy units) :l
Tequil = equilibration temperature (temperature units) :l
Dcut = minimum distance between boosted bonds (distance units) :l
alpha = boostostat relaxation time (time units) :l
boost = desired time boost factor (unitless) :l
Btarget = desired time boost factor (unitless) :l
zero or more keyword/value pairs may be appended :l
keyword = {histo} or {lost} or {check/bias} or {check/coeff}
{histo} values = Nevery Nbin delta Nout
Nevery = histogram bond bias coefficients every this many timesteps
Nbin = # of histogram bins
delta = width of each histogram bin
Nout = output histogram counts every this many timesteps
keyword = {lost} or {check/bias} or {check/coeff}
{lostbond} value = error/warn/ignore
{check/bias} values = Nevery error/warn/ignore
{check/coeff} values = Nevery error/warn/ignore :pre
@ -61,26 +56,26 @@ compared to running for the same number of timesteps with normal MD.
See the "hyper"_hyper.html doc page for a more general discussion of
hyperdynamics and citations that explain both GHD and LHD.
The equations and logic used by this fix to perform LHD follow the
description given in "(Voter2013)"_#Voter2013lhd. The bond-boost form
of a bias potential for HD is due to Miron and Fichthorn as described
in "(Miron)"_#Mironlhd.
The equations and logic used by this fix and described here to perform
LHD follow the description given in "(Voter2013)"_#Voter2013lhd. The
bond-boost form of a bias potential for HD is due to Miron and
Fichthorn as described in "(Miron)"_#Mironlhd.
To understand this description, you should first read the description
of the GHD algorithm on the "fix hyper/global"_fix_hyper/global.html
doc page. This description of LHD builds on the GHD description.
The definition of bonds, Eij, and Emax are the same for GHD and LHD.
The formula for Vij(max) and Fij(max) are also the same except for a
The formulas for Vij(max) and Fij(max) are also the same except for a
pre-factor Cij, explained below.
The bias energy Vij applied to a bond IJ with maximum strain is
Vij(max) = Cij * Vmax * (1 - (Eij/q)^2) for abs(Eij) < qfactor
= 0 otherwise :pre
= 0 otherwise :pre
The derivative of Vij(max) with respect to the position of each atom
in the IJ bond gives a force Fij(max) acting on the bond as
in the IJ bond gives a bias force Fij(max) acting on the bond as
Fij(max) = - dVij(max)/dEij = 2 Cij Vmax Eij / qfactor^2 for abs(Eij) < qfactor
= 0 otherwise :pre
@ -88,22 +83,76 @@ Fij(max) = - dVij(max)/dEij = 2 Cij Vmax Eij / qfactor^2 for abs(Eij) < qfacto
which can be decomposed into an equal and opposite force acting on
only the two I,J atoms in the IJ bond.
The key difference is that in GHD this bias energy and force is added
(on a particular timestep) to only one bond (pair of atoms) in the
system, the bond which has maximum strain Emax.
The key difference is that in GHD a bias energy and force is added (on
a particular timestep) to only one bond (pair of atoms) in the system,
which is the bond with maximum strain Emax.
In LHD, the bias energy and force
In LHD, a bias energy and force can be added to multiple bonds
separated by the specified {Dcut} distance or more. A bond IJ is
biased if it is the maximum strain bond within its local
"neighborhood", which is defined as the bond IJ plus any neighbor
bonds within a distance {Dcut} from IJ. The "distance" between bond
IJ and bond KL is the minimum distance between any of the IK, IL, JK,
JL pairs of atoms.
For a large system, multiple bonds will typically meet this
requirement, and thus a bias potential Vij(max) will be applied to
many bonds on the same timestep.
NOTE: Add eqs for boostostat and boost coeff.
Explain how many bonds are boosted simultaneously
and how to choose boost factor and initial Vmax.
In LHD, all bonds store a Cij prefactor which appears in the Vij(max)
and Fij(max) equations above. Note that the Cij factor scales the
strength of the bias energy and forces whenever bond IJ is the maximum
strain bond in its neighborhood.
Cij is initialized to 1.0 when a bond between the I,J atoms is first
defined. The specified {Btarget} factor is then used to adjust the
Cij prefactors for each bond every timestep in the following manner.
An instantaneous boost factor Bij is computed each timestep
for each bond, as
Bij = exp(beta * Vkl(max)) :pre
where Vkl(max) is the bias energy of the maxstrain bond KL within bond
IJ's neighborhood, beta = 1/kTequil, and {Tequil} is the temperature
of the system and an argument to this fix.
NOTE: To run LHD, the input script must also use the "fix
langevin"_fix_langevin.html command to thermostat the atoms at the
same {Tequil} as specified by this fix, so that the system is running
constant-temperature (NVT) dynamics. LAMMPS does not check that this
is done.
Note that if IJ = KL, then bond IJ is a biased bond on that timestep,
otherwise it is not. But regardless, the boost factor Bij can be
thought of an estimate of time boost currently being applied within a
local region centered on bond IJ. For LHD, we want this to be the
specified {Btarget} value everywhere in the simulation domain.
To accomplish this, if Bij < Btarget, the Cij prefactor for bond IJ is
incremented on the current timestep by an amount proportional to the
inverse of the specified {alpha} and the difference (Bij - Btarget).
Conversely if Bij > Btarget, Cij is decremented by the same amount.
This procedure is termed "boostostatting" in
"(Voter2013)"_#Voter2013lhd. It drives all of the individual Cij to
values such that when Vij{max} is applied as a bias to bond IJ, the
resulting boost factor Bij will be close to {Btarget} on average.
Thus the LHD time acceleration factor for the overall system is
effectively {Btarget}.
Note that in LHD, the boost factor {Btarget} is specified by the user.
This is in contrast to global hyperdynamics (GHD) where the boost
factor varies each timestep and is computed as a function of {Vmax},
Emax, and {Tequil}; see the "fix hyper/global"_fix_hyper_global.html
doc page for details.
:line
Here is additional information on the input parameters for GHD.
Here is additional information on the input parameters for LHD.
Note that the {cutbond}, {qfactor}, and {Tequil} arguments have the
same meaning as for GHD. The {Vmax} argument is slightly different.
The {Dcut}, {alpha}, and {Btarget} parameters are unique to LHD.
The {cutbond} argument is the cutoff distance for defining bonds
between pairs of nearby atoms. A pair of I,J atoms in their
@ -120,56 +169,92 @@ equilibrium (quenched) R0ij distance and the two atoms in the bond
could still experience a non-zero bias force.
If {qfactor} is set too large, then transitions from one energy basin
to another are affected because the bias potneial is non-zero at the
to another are affected because the bias potential is non-zero at the
transition state (e.g. saddle point). If {qfactor} is set too small
than little boost can be achieved because the Eij strain of some bond in
the system will (nearly) always exceed {qfactor). A value of 0.3 for
{qfactor} is typically a reasonable value.
What about requested boost
The {Vmax} argument is a fixed prefactor on the bias potential. There
is a also a dynamic prefactor Cij, driven by the choice of {Btarget}
as discussed above. The product of these should be a value less than
the smallest barrier height for an event to occur. Otherwise the
applied bias potential may be large enough (when added to the
interatomic potential) to produce a local energy basin with a maxima
in the center. This can produce artificial energy minima in the same
basin that trap an atom. Or if Cij*{Vmax} is even larger, it may
induce an atom(s) to rapidly transition to another energy basin. Both
cases are "bad dynamics" which violate the assumptions of LHD that
guarantee an accelerated time-accurate trajectory of the system.
The {Vmax} argument is the prefactor on the bias potential. It should
be set to a value less than the smallest barrier height for an event
to occur. Otherwise the applied bias potential may be large enough
(when added to the interatomic potential) to produce a local energy
basin with a maxima in the center. This can produce artificial energy
minima in the same basin that trap an atom. Or if {Vmax} is even
larger, it may induce an atom(s) to rapidly transition to another
energy basin. Both cases are "bad dynamics" which violate the
assumptions of GHD that guarantee an accelerated time-accurate
trajectory of the system.
NOTE: It may seem that {Vmax} can be set to any value, and Cij will
compensate to reduce the overall prefactor if necessary. However the
Cij are initialized to 1.0 and the boostostatting procedure typically
operates slowly enough that there can be a time period of bad dynamics
if {Vmax} is set too large. A better strategy is to set {Vmax} to the
smallest barrier height for an event (the same as for GHD), so that
the Cij remain near unity.
The {Tequil} argument is the temperature at which the system is
simulated; see the comment above about the "fix
langevin"_fix_langevin.html thermostatting. It is also part of the
beta term in the exponential factor that determines how much boost is
achieved as a function of the bias potential.
achieved as a function of the bias potential. See the discussion of
the {Btarget} argument below.
In general, the lower the value of {Tequil} and the higher the value
of {Vmax}, the more boost will be achieved by the GHD algorithm.
As discussed above, the {Dcut} argument is the distance required
between two locally maxstrain bonds for them to both be selected as
biased bonds on the same timestep. Computationally, the larger {Dcut}
is, the more work (computation and communication) must be done each
timestep within the LHD algorithm. And the fewer bonds can be
simultaneously biased, which may mean the specified {Btarget} time
acceleration cannot be achieved.
The {Dcut} argument is the distance required between two bonds for
them to be selected as both being boosted. The distance is between
the center points of each bond. Actually between any pair of atoms in
either bond.
Physically {Dcut} should be a long enough distance that biasing two
pairs of atoms that close together will not influence the dynamics of
each pair. E.g. something like 2x the cutoff of the interatomic
potential. In practice a {Dcut} value of ~10 Angstroms seems to work
well for many solid-state systems.
The {alpha} argument is a pre-factor on the update equation
for each atom's boostostat:
As described above the {alpha} argument is a pre-factor in the
boostostat update equation for each bond's Cij prefactor. {Alpha} is
specified in time units, similar to other thermostat or barostat
damping parameters. It is roughly the physical time it will take the
boostostat to adjust a Cij value from a too high (or too low) value to
a correct one. An {alpha} setting of a few ps is typically good for
solid-state systems. Note that the {alpha} argument here is the
inverse of the alpha parameter discussed in
"(Voter2013)"_#Voter2013lhd.
NOTE: give equation above
The {Btarget} argument is the desired time boost factor (a value > 1)
that all the atoms in the system will experience. The elapsed time
t_hyper for an LHD simulation running for {N} timesteps is simply
Note that the units for {alpha} are in time units, similar
to other thermostat or barostat damping parameters
t_hyper = Btarget * N*dt :pre
The {boost} argument is the desired boost factor (e.g. 100x)
that all the atoms in the system will experience.
where dt is the timestep size defined by the "timestep"_timestep.html
command. The effective time acceleration due to LHD is thus t_hyper /
N*dt = Btarget, where N*dt is elapsed time for a normal MD run
of N timesteps.
NOTE: explain how to choose good value for this. If this parameter is
chosen to be too large, then the bias potential applied to pairs of
bonded atoms may become unphysically large, leading to bad dynamics.
If chosen too small, the hyperdynamics run may be inefficient in the
sense that events take a long time to occur.
You cannot choose an arbitrarily large setting for {Btarget}. The
maximum value you should choose is
Btarget = exp(beta * Vsmall) :pre
where Vsmall is the smallest event barrier height in your system, beta
= 1/kTequil, and {Tequil} is the specified temperature of the system
(both by this fix and the Langevin thermostat).
Note that if {Btarget} is set smaller than this, the LHD simulation
will run correctly. There will just be fewer events because the hyper
time (t_hyper equation above) will be shorter.
NOTE: If you have no physical intuition as to the smallest barrier
height in your system, a reasonable strategy to determine the largest
{Btarget} you can use for an LHD model, is to run a sequence of
simulations with smaller and smaller {Btarget} values, until the event
rate does not change.
:line
@ -189,21 +274,20 @@ the bias potential (energy units) applied on the current timestep,
summed over all biased bonds. The vector stores the following
quantities:
1 = # of boosted bonds on this step
2 = max strain of any bond on this step (unitless)
3 = average bias potential for all bonds on this step (energy units)
4 = average bonds/atom on this step
5 = average neighbor bonds/bond on this step :ul
1 = # of biased bonds on this step
2 = max strain Eij of any bond on this step (unitless)
3 = average bias potential for all biased bonds on this step (energy units)
4 = average # of bonds/atom on this step
5 = average neighbor bonds/bond on this step within {Dcut} :ul
6 = fraction of steps and bonds with no bias during this run
7 = max drift distance of any atom during this run (distance units)
8 = max bond length during this run (distance units)
9 = average # of boosted bonds/step during this run
10 = average bias potential for all bonds during this run (energy units)
11 = max bias potential for any bond during this run (energy units)
12 = min bias potential for any bond during this run (energy units)
13 = max distance from my box of any ghost atom with maxstrain < qfactor during th
is run (distance units)
9 = average # of biased bonds/step during this run
10 = average bias potential for all biased bonds during this run (energy units)
11 = max bias potential for any biased bond during this run (energy units)
12 = min bias potential for any biased bond during this run (energy units)
13 = max distance from my sub-box of any ghost atom with maxstrain < qfactor during this run (distance units)
14 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units)
15 = count of ghost neighbor atoms not found on reneighbor steps during this run
16 = count of lost bond partners during this run
@ -222,13 +306,60 @@ hyper run is performed. Quantities 20-23 are cummulative across
multiple runs (since the fix was defined in the input script).
For value 6, the numerator is a count of all biased bonds on every
timestep whose bias value = 0.0. The denominator is the count of all
biased bonds on all timesteps.
timestep whose bias energy = 0.0 due to Eij >= {qfactor}. The
denominator is the count of all biased bonds on all timesteps.
For value 7, drift is the distance an atom moves between timesteps
when the bond list is reset, i.e. between events. Atoms involved in
an event will typically move the greatest distance since others are
typically oscillating around their lattice site.
For values 13 and 14, the maxstrain of a ghost atom is the maxstrain
of any bond it is part of, and it is checked for ghost atoms within
the bond neighbor cutoff.
Values 15-19 are mostly useful for debugging and diagnositc purposes.
For values 15-17, it is possible that a ghost atom owned by another
processor will move far enough (e.g. as part of an event-in-progress)
that it will no longer be within the communication cutoff distance for
acquiring ghost atoms. Likewise it may be a ghost atom bond partner
that cannot be found because it has moved too far. These values count
those occurrences. Because they typically involve atoms that are part
of events, they do not usually indicate bad dynamics. Value 16 is the
average bias coefficient for bonds where a partner atom was lost.
For value 18, no two bonds should be biased if they are within a
{Dcut} distance of each other. This value should be zero, indicating
that no pair of bonds "overlap", meaning they are closer than {Dcut}
from each other.
For value 19, the same bias coefficient is stored by both atoms in an
IJ bond. This value should be zero, indicating that for all bonds,
each atom in the bond stores the a bias coefficient with the same
value.
Value 20 is simply the specified {boost} factor times the number of
timestep times the timestep size.
For value 21, events are checked for by the "hyper"_hyper.html command
once every {Nevent} timesteps. This value is the count of those
timesteps on which one (or more) events was detected. It is NOT the
number of distinct events, since more than one event may occur in the
same {Nevent} time window.
For value 22, each time the "hyper"_hyper.html command checks for an
event, it invokes a compute to flag zero or more atoms as
participating in one or more events. E.g. atoms that have displaced
more than some distance from the previous quench state. Value 22 is
the cummulative count of the number of atoms participating in any of
the events that were found.
Value 23 tallies the number of new bonds created by the bond reset
operation. Bonds between a specific I,J pair of atoms may persist for
the entire hyperdynamics simulation if neither I or J are involved in
an event.
The scalar and vector values calculated by this fix are all
"intensive".
@ -256,4 +387,3 @@ LAMMPS"_Section_start.html#start_3 section for more info.
:link(Mironlhd)
[(Miron)] R. A. Miron and K. A. Fichthorn, J Chem Phys, 119, 6210 (2003).