Merge pull request #1379 from lammps/hyper

e.g. {bond/local}.  Styles with neither "atom" or "local" in their

style produce global quantities.

Note that a single compute produces either global or per-atom or local

quantities, but never more than one of these (with only a few

exceptions, as documented by individual compute commands).

Note that a single compute can produce either global or per-atom or

local quantities, but not both global and per-atom.  It can produce

local quantities in tandem with global or per-atom quantities.  The

compute doc page will explain.

Global, per-atom, and local quantities each come in three kinds: a

single scalar value, a vector of values, or a 2d array of values.  The

each processor based on the atoms it owns, but there may be zero or

more per atoms.

Note that a single fix may produces either global or per-atom or local

quantities (or none at all), but never more than one of these.

Note that a single fix can produce either global or per-atom or local

quantities (or none at all), but not both global and per-atom.  It can

produce local quantities in tandem with global or per-atom quantities.

The fix doc page will explain.

Global, per-atom, and local quantities each come in three kinds: a

single scalar value, a vector of values, or a 2d array of values.  The

  {mode} arg = {scalar} or {vector}

    scalar = all input values are scalars

    vector = all input values are vectors

  {kind} arg = {global} or {peratom} or {local}

  {file} arg = filename

    filename = name of file to output histogram(s) to

  {ave} args = {one} or {running} or {window}

different kinds (e.g. global and per-atom) cannot be mixed.  Atom

attributes are per-atom vector values.  See the doc page for

individual "compute" and "fix" commands to see what kinds of

quantities they generate.

quantities they generate.  See the optional {kind} keyword below for

how to force the fix ave/histo command to disambiguate if necessary.

Note that the output of this command is a single histogram for all

input values combined together, not one histogram per input value.

per-atom or local vectors, or columns of global or per-atom or local

arrays.

The {kind} keyword only needs to be set if a compute or fix produces

more than one kind of output (global, per-atom, local).  If this is

not the case, then LAMMPS will determine what kind of input is

provided and whether all the input arguments are consistent.  If a

compute or fix produces more than one kind of output, the {kind}

keyword should be used to specify which output will be used.  The

remaining input arguments must still be consistent.

The {beyond} keyword determines how input values that fall outside the

{lo} to {hi} bounds are treated.  Values such that {lo} <= value <=

{hi} are assigned to one bin.  Values on a bin boundary are assigned

first bin and values > {hi} are counted in the last bin.  If {beyond}

is set to {extend} then two extra bins are created, so that there are

Nbins+2 total bins.  Values < {lo} are counted in the first bin and

values > {hi} are counted in the last bin (Nbins+1).  Values between

values > {hi} are counted in the last bin (Nbins+2).  Values between

{lo} and {hi} (inclusive) are counted in bins 2 through Nbins+1.  The

"coordinate" stored and printed for these two extra bins is {lo} and

{hi}.

[Default:] none

The option defaults are mode = scalar, ave = one, start = 0, no file

output, beyond = ignore, and title 1,2,3 = strings as described above.

The option defaults are mode = scalar, kind = figured out from input

arguments, ave = one, start = 0, no file output, beyond = ignore, and

title 1,2,3 = strings as described above.

where beta = 1/kTequil, and {Tequil} is the temperature of the system

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

NOTE: To run a GHD simulation, 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

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

{Vmax} you can use for a GHD model, is to run a sequence of

simulations with smaller and smaller {Vmax} values, until the event

rate does not change.

rate does not change (as a function of hyper time).

The {Tequil} argument is the temperature at which the system is

simulated; see the comment above about the "fix

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 achievable by the GHD algorithm.

of {Vmax}, the more time boost will be achievable by the GHD

algorithm.

:line

fix to add the energy of the bias potential to the the system's

potential energy as part of "thermodynamic output"_thermo_style.html.

This fix computes a global scalar and global vector of length 11, which

This fix computes a global scalar and global vector of length 12, which

can be accessed by various "output commands"_Howto_output.html.  The

scalar is the magnitude of the bias potential (energy units) applied on

the current timestep.  The vector stores the following quantities:

1 = boost factor on this step (unitless)

2 = max strain Eij of any bond on this step (unitless)

2 = max strain Eij of any bond on this step (absolute value, unitless)

3 = ID of first atom in the max-strain bond

4 = ID of second atom in the max-strain bond

5 = average # of bonds/atom on this step :ul

6 = fraction of timesteps with bias = 0.0 during this run

7 = max drift distance of any atom during this run (distance units)

8 = max bond length during this run (distance units) :ul

6 = fraction of timesteps where the biased bond has bias = 0.0 during this run

7 = fraction of timesteps where the biased bond has negative strain during this run

8 = max drift distance of any atom during this run (distance units)

9 = max bond length during this run (distance units) :ul

9 = cumulative hyper time since fix was defined (time units)

10 = cumulative count of event timesteps since fix was defined

11 = cumulative count of atoms in events since fix was defined :ul

10 = cumulative hyper time since fix was defined (time units)

11 = cumulative count of event timesteps since fix was defined

12 = cumulative count of atoms in events since fix was defined :ul

The first 5 quantities are for the current timestep.  Quantities 6-8

are for the current hyper run.  Quantities 9-11 are cumulative across

multiple runs (since the fix was defined in the input script).

The first 5 quantities are for the current timestep.  Quantities 6-9

are for the current hyper run.  They are reset each time a new hyper

run is performed.  Quantities 19-12 are cumulative across multiple

runs (since the point in the input script the fix was defined).

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 8, drift is the distance an atom moves between two quenched

states when the second quench determines an event has occurred.  Atoms

involved in an event will typically move the greatest distance since

others typically remain near their original quenched position.

For value 10, events are checked for by the "hyper"_hyper.html command

For value 11, 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 11, each time the "hyper"_hyper.html command checks for an

For value 12, 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 11 is

alpha = boostostat relaxation time (time units) :l

Btarget = desired time boost factor (unitless) :l

zero or more keyword/value pairs may be appended :l

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

keyword = {check/ghost} or {check/bias} :l

  {check/ghost} values = none

  {check/bias} values = Nevery error/warn/ignore :pre

:ule

[Examples:]

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 formulas for Vij(max) and Fij(max) are also the same except for a

The definition of bonds and Eij are the same for GHD and LHD.  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

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: To run an LHD simulation, 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

potential.  In practice a {Dcut} value of ~10 Angstroms seems to work

well for many solid-state systems.

NOTE: You must also insure that ghost atom communication is performed

for a distance of at least {Dcut} + {cutevent} where {cutevent} = the

distance one or more atoms move (between quenched states) to be

considered an "event".  It is an argument to the "compute

event/displace" command used to detect events.  By default the ghost

communication distance is set by the pair_style cutoff, which will

typically be < {Dcut}.  The "comm_modify cutoff"_comm_modify.html

command can be used to set the ghost cutoff explicitly, e.g.

NOTE: You should insure that ghost atom communication is performed for

a distance of at least {Dcut} + {cutevent} = the distance one or more

atoms move (between quenched states) to be considered an "event".  It

is an argument to the "compute event/displace" command used to detect

events.  By default the ghost communication distance is set by the

pair_style cutoff, which will typically be < {Dcut}.  The "comm_modify

cutoff"_comm_modify.html command should be used to override the ghost

cutoff explicitly, e.g.

comm_modify cutoff 12.0 :pre

This fix does not know the {cutevent} parameter, but uses half the

bond length as an estimate to warn if the ghost cutoff is not long

enough.

Note that this fix does not know the {cutevent} parameter, but uses

half the {cutbond} parameter as an estimate to warn if the ghost

cutoff is not long enough.

As described above the {alpha} argument is a pre-factor in the

boostostat update equation for each bond's Cij prefactor.  {Alpha} is

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.

rate does not change (as a function of hyper time).

:line

Here is additional information on the optional keywords for this fix.

The {check/ghost} keyword turns on extra computation each timestep to

compute statistics about ghost atoms used to determine which bonds to

bias.  The output of these stats are the vector values 14 and 15,

described below.  If this keyword is not enabled, the output

of the stats will be zero.

The {check/bias} keyword turns on extra computation and communication

to check if any biased bonds are closer than {Dcut} to each other,

which should not be the case if LHD is operating correctly.  Thus it

is a debugging check.  The {Nevery} setting determines how often the

check is made.  The {error}, {warn}, or {ignore} setting determines

what is done if the count of too-close bonds is not zero.  Either the

code will exit, or issue a warning, or silently tally the count.  The

count can be output as vector value 17, as described below.  If this

keyword is not enabled, the output of that statistic will be 0.

Note that both of these computations are costly, hence they are only

enabled by these keywords.

:line

fix to add the energy of the bias potential to the the system's

potential energy as part of "thermodynamic output"_thermo_style.html.

This fix computes a global scalar and global vector of length 23,

which can be accessed by various "output

commands"_Howto_output.html.  The scalar is the magnitude of

the bias potential (energy units) applied on the current timestep,

summed over all biased bonds.  The vector stores the following

quantities:

This fix computes a global scalar and global vector of length 21,

which can be accessed by various "output commands"_Howto_output.html.

The scalar is the magnitude of the bias potential (energy units)

applied on the current timestep, summed over all biased bonds.  The

vector stores the following quantities:

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)

2 = max strain Eij of any bond on this step (absolute value, unitless)

3 = average bias coeff for all bonds on this step (unitless)

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 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

17 = average bias coeff for lost bond partners during this run

18 = count of bias overlaps found during this run

19 = count of non-matching bias coefficients found during this run :ul

20 = cumulative hyper time since fix created (time units)

21 = cumulative count of event timesteps since fix created

22 = cumulative count of atoms in events since fix created

23 = cumulative # of new bonds since fix created :ul

6 = max bond length during this run (distance units)

7 = average # of biased bonds/step during this run

8 = fraction of biased bonds with no bias during this run

9 = fraction of biased bonds with negative strain during this run

10 = average bias coeff for all bonds during this run (unitless)

11 = min bias coeff for any bond during this run (unitless)

12 = max bias coeff for any bond during this run (unitless)

13 = max drift distance of any bond atom during this run (distance units)

14 = max distance from proc subbox of any ghost atom with maxstrain < qfactor during this run (distance units)

15 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units)

16 = count of ghost atoms that could not be found on reneighbor steps during this run

17 = count of bias overlaps (< Dcut) found during this run

18 = cumulative hyper time since fix created (time units)

19 = cumulative count of event timesteps since fix created

20 = cumulative count of atoms in events since fix created

21 = cumulative # of new bonds formed since fix created :ul

The first quantities (1-5) are for the current timestep.  Quantities

6-19 are for the current hyper run.  They are reset each time a new

hyper run is performed.  Quantities 20-23 are cumulative across

multiple runs (since the fix was defined in the input script).

6-17 are for the current hyper run.  They are reset each time a new

hyper run is performed.  Quantities 18-21 are cumulative across

multiple runs (since the point in the input script the fix was

defined).

For value 6, the numerator is a count of all biased bonds on every

For value 8, the numerator is a count of all biased bonds on each

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 diagnostic 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

For value 9, the numerator is a count of all biased bonds on each

timestep with negative strain.  The denominator is the count of all

biased bonds on all timesteps.

Values 13-17 are mostly useful for debugging and diagnostic purposes.

For value 13, drift is the distance an atom moves between two quenched

states when the second quench determines an event has occurred.  Atoms

involved in an event will typically move the greatest distance since

others typically remain near their original quenched position.

For values 14-16, neighbor atoms in the full neighbor list with cutoff

{Dcut} may be ghost atoms outside a processor's sub-box.  Before the

next event occurs they may move further than {Dcut} away from the

sub-box boundary.  Value 14 is the furthest (from the sub-box) any

ghost atom in the neighbor list with maxstrain < {qfactor} was

accessed during the run.  Value 15 is the same except that the ghost

atom's maxstrain may be >= {qfactor}, which may mean it is about to

participate in an event.  Value 16 is a count of how many ghost atoms

could not be found on reneighbor steps, presumably because they moved

too far away due to their participation in an event (which will likely

be detected at the next quench).

Typical values for 14 and 15 should be slightly larger than {Dcut},

which accounts for ghost atoms initially at a {Dcut} distance moving

thermally before the next event takes place.

Note that for values 14 and 15 to be computed, the optional keyword

{check/ghost} must be specified.  Otherwise these values will be zero.

This is because computing them incurs overhead, so the values are only

computed if requested.

Value 16 should be zero or small.  As explained above a small count

likely means some ghost atoms were participating in their own events

and moved a longer distance.  If the value is large, it likely means

the communication cutoff for ghosts is too close to {Dcut} leading to

many not-found ghost atoms before the next event.  This may lead to a

reduced number of bonds being selected for biasing, since the code

assumes those atoms are part of highly strained bonds.  As explained

above, the "comm_modify cutoff"_comm_modify.html command can be used

to set a longer cutoff.

For value 17, 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.

that no pair of biased bonds are closer than {Dcut} from each other.

Note that for values 17 to be computed, the optional keyword

{check/bias} must be specified and it determines how often this check

is performed.  This is because performing the check incurs overhead,

so if only computed as often as requested.

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.

The result at the end of the run is the cumulative total from every

timestep the check was made.  Note that the value is a count of atoms

in bonds which found other atoms in bonds too close, so it is almost

always an over-count of the number of too-close bonds.

Value 20 is simply the specified {boost} factor times the number of

timestep times the timestep size.

Value 18 is simply the specified {boost} factor times the number of

timesteps times the timestep size.

For value 21, events are checked for by the "hyper"_hyper.html command

For value 19, 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

For value 20, 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

more than some distance from the previous quench state.  Value 20 is

the cumulative 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

Value 21 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".

This fix also computes a local vector of length the number of bonds

currently in the system.  The value for each bond is its Cij prefactor

(bias coefficient).  These values can be can be accessed by various

"output commands"_Howto_output.html.  A particularly useful one is the

"fix ave/histo"_fix_ave_histo.html command which can be used to

histogram the Cij values to see if they are distributed reasonably

close to 1.0, which indicates a good choice of {Vmax}.

The local values calculated by this fix are unitless.

No parameter of this fix can be used with the {start/stop} keywords of

the "run"_run.html command.  This fix is not invoked during "energy

minimization"_minimize.html.

"hyper"_hyper.html, "fix hyper/global"_fix_hyper_global.html

[Default:] None

[Default:]

The check/ghost and check/bias keywords are not enabled by default.

:line