more edits to hyper docs

doc/src/fix_hyper_global.txt

@@ -26,59 +26,122 @@ fix 1 all hyper/global 1.0 0.3 0.8 300.0 :pre
 [Description:]
 
 This fix is meant to be used with the "hyper"_hyper.html command to
-perform a bond-boost hyperdynamics simulation.  The role of this fix
-is a select a single pair of atoms within the system at each timestep
-to add a global bias potential to, which will alter their dynamics.
+perform a bond-boost global hyperdynamics (GHD) simulation.  The role
+of this fix is to a select a single pair of atoms in the system at
+each timestep to add a global bias potential to, which will alter the
+dynamics of the system in a manner that effectively accelerates time.
 This is in contrast to the "fix hyper/local"_fix_hyper_local.html
-command, which can add a local bias potential to multiple pairs of
-atoms at each timestep.
+command, which can be user to perform a local hyperdynamics (LHD)
+simulation, by adding a local bias potential to multiple pairs of
+atoms at each timestep.  GHD can speed-up a small simulation with up
+to a few 100 atoms.  For larger systems, LHD is needed to achieve good
+speed-ups.
 
 For a system that undergoes rare transition events, where one or more
-atoms move across an energy barrier to a new potential energy basin,
-the effect of the bias potential is to induce more rapid transitions.
+atoms move over an energy barrier to a new potential energy basin, the
+effect of the bias potential is to induce more rapid transitions.
 This can lead to a dramatic speed-up in the rate at which events
 occurs, without altering their relative frequencies, thus leading to
-an overall dramatic speed-up in the effective elapsed time of the
-simulation.
+an overall increase in the elapsed real time of the simulation as
+compared to simulating 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 are as
+follows.  This follows the description of GHD given in
+"(Voter2013)"_#Voter2013ghd.  The bond-boost form of a bias potential
+for GHD is due to Miron and Fichthorn as described in
+"(Miron)"_#Mironghd.  Here we use a simplified version of bond-boost
+GHD where a single bond is biased at any one timestep.
+
+Bonds are defined between every pair of I,J atoms whose R0ij distance
+is less than {cutbond}, when the system is in a quenched state
+(minimum) energy.  Note that these are not "bonds" in a covalent
+sense.  A bond is simply any pair of atoms that meet the distance
+criterion.  {Cutbond} is an argument to this fix; it is discussed more
+below.
+
+The current strain of bond IJ (when running dynamics) is defined as
 
-Cite various papers.
+Eij = (Rij - R0ij) / R0ij :pre
 
-The current strain of a bond IJ is defined as
+where Rij is the current distance between atoms I,J, and R0ij is the
+equilibrium distance in the quenched state.
 
-Eij = (Rij - R0ij) / R0ij :pre
+The bias energy Vij of any bond IJ is defined as
+
+Vij = Vmax * (1 - (Eij/q)^2) for abs(Eij) < qfactor
+    = 0 otherwise :pre
+
+where the prefactor {Vmax} and the cutoff {qfactor} are arguments to
+this fix; they are discussed more below.  This functional form is an
+inverse parabola centered at 0.0 with height Vmax and which goes to
+0.0 at +/- qfactor.
+
+Let Emax = the maximum of abs(Eij) for all IJ bonds in the system on a
+given timestep.  On that step, Vij is added as a bias potential to
+only the two atoms IJ in the bond with strain Emax, call it Vij(max).
+Note that Vij(max) will be 0.0 if Emax >= qfactor on that timestep.
+Also note that Vij(max) is added to the normal interatomic potential
+that is computed between all atoms in the system.
+
+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
+
+Fij(max) = - dVij(max)/dEij = 2 Vmax Eij / qfactor^2   for abs(Eij) < qfactor
+         = 0 otherwise :pre
 
-Emax = is the max of the absolute value of Eij for all IJ bonds.
+which can be decomposed into an equal and opposite force acting (only)
+on the two I,J atoms in the Emax bond.
 
-dVij = Vmax * (1 - (Eij/q)^2) for abs(Eij) < q
-     = 0  otherwise :pre
+The boost factor in time for the system is given each timestep I by
 
-Delta Vbias = minimum of dVij for all IJ bonds :pre
+Bi = exp(beta * Vij(max)) :pre
 
-The boost factor B = exp(beta * Delta Vbias)
-for a single timestep.
+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.
 
-Thus the accumulated hypertime is simply
+NOTE: Must use langevin set to Tequil also.  LAMMPS does not
+check that you do this.
+as maintained by a thermostat for constant NVT dynamics. 
 
-t_hyper = Sum (i = 1 to Nsteps) Bi * dt :pre
+The elapsed system 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.
+
+Note that in global hyperdynamics, 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.
+
+:line
 
 The {cutbond} argument is the cutoff distance for defining bonds
 between pairs of nearby atoms.  A pair of atoms in their equilibrium,
-minimum-energy config, which are a distance Rij < cutbond, are
-defined as a bonded pair.
+minimum-energy config, which are a distance Rij < cutbond, are defined
+as a bonded pair.  This should normally be roughly ~25% larger than
+the nearest-neighbor distance in a crystalline lattice.
 
-The {qfactor} argument is the limiting strain at which
-the Vbias (the bias potential) goes to 0.0.
+The {qfactor} argument is the limiting strain at which the bias
+potential goes to 0.0.
 
-If {qfactor} is too big, then transitions are affected b/c
-the bias energy is non-zero at the transitions.  If it is
-too small than not must boost is achieved for a large system
-with many bonds (some bonds Eij always exceeds qfactor).
+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
+transitions.  If {qfactor} is set too small than little boost is
+achieved because some the Eij strain of some bond is (nearly) always
+exceeding {qfactor).  Note
 
 The {Vmax} argument is the prefactor on the bias potential.
+It should be set to a value smaller than 
 
-The {Tequil} argument is part of the beta term in the
-exponential factor that determines how much boost is
-achieved as a function of the bias potential.
+The {Tequil} argument is part of the beta term in the exponential
+factor that determines how much boost is achieved as a function of the
+bias potential.
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
@@ -96,12 +159,12 @@ 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 of any bond on this step (unitless)
+2 = max strain Eij of any bond on this step (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 step with no bias during this run
+6 = fraction of steps 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
 
@@ -135,12 +198,19 @@ minimization"_minimize.html.
 
 [Restrictions:]
 
-This fix is part of the REPLICA package.  It is only enabled if LAMMPS
-was built with that package.  See the "Making
-LAMMPS"_Section_start.html#start_3 section for more info.
+This command can only be used if LAMMPS was built with the REPLICA
+package.  See the "Build package"_Build_package.html doc
+page for more info.
 
 [Related commands:]
 
-"fix hyper/local"_fix_hyper_local.html
+"hyper"_hyper.html, "fix hyper/local"_fix_hyper_local.html
 
 [Default:] None
+
+:link(Voter2013)
+[(Voter2013)] S. Y. Kim, D. Perez, A. F. Voter, J Chem Phys, 139,
+144110 (2013).
+
+:link(Miron)
+[(Miron)] R. A. Miron and K. A. Fichthorn, J Chem Phys, 119, 6210 (2003).

doc/src/hyper.txt

@@ -10,14 +10,14 @@ hyper command :h3
 
 [Syntax:]
 
-hyper N Nevent fixID computeID keyword values ... :pre
+hyper N Nevent fix-ID compute-ID keyword values ... :pre
 
 N = # of timesteps to run :ulb,l
 Nevent = check for events every this many steps :l
-fixID = ID of a fix that applies a global or local bias potential :l
-computeID = ID of a compute that identifies when an event has occurred :l
+fix-ID = ID of a fix that applies a global or local bias potential, can be NULL :l
+compute-ID = ID of a compute that identifies when an event has occurred :l
 zero or more keyword/value pairs may be appended :l
-keyword = {min} or {time} :l
+keyword = {min} or {time} or {dump} or {rebond} :l
   {min} values = etol ftol maxiter maxeval
     etol = stopping tolerance for energy, used in quenching
     ftol = stopping tolerance for force, used in quenching
@@ -25,82 +25,169 @@ keyword = {min} or {time} :l
     maxeval = max number of force/energy evaluations, used in quenching
   {time} value = {steps} or {clock}
     {steps} = simulation runs for N timesteps (default)
-    {clock} = simulation runs until hyper time exceeds N timesteps :pre
+    {clock} = simulation runs until hyper time exceeds N timesteps
+  {dump} value = dump-ID
+    dump-ID = ID of dump to trigger whenever an event takes place
+  {rebond} value = Nrebond
+    Nrebond = frequency at which to reset bonds, even if no event has occurred
+ :pre
 :ule
 
 [Examples:]
 
-hyper 5000 100 global event :pre
+compute event all event/displace 1.0
+fix HG mobile hyper/global 3.0 0.3 0.4 800.0
+hyper 5000 100 HG event min 1.0e-6 1.0e-6 100 100 dump 1 dump 5 :pre
 
 [Description:]
 
 Run a bond-boost hyperdynamics (HD) simulation where time is
-accelerated by application of a bias potential to a one or more pairs
-of atoms in the system.  This command can be used to run both global
-and local hyperdyamics.  In global HD a single bond (nearby pair of
-atoms) within the system is biased on a given timestep.  In local HD
-multiple bonds (separated by a sufficient distance) can be biased
-simultaneously at each timestep.
-
-Make next paragraph more HD-specific.
-
-Both global and local HD are described in "this paper"_#Voter2013 by
-Art Voter and collaborators.  They are methods for performing
-accelerated dynamics that is suitable for infrequent-event systems
-that obey first-order kinetics.  A good overview of accelerated
-dynamics methods for such systems in given in "this review
-paper"_#Voter2002prd from the same group.  To quote from the paper:
-"The dynamical evolution is characterized by vibrational excursions
-within a potential basin, punctuated by occasional transitions between
-basins."  The transition probability is characterized by p(t) =
-k*exp(-kt) where k is the rate constant.  Running multiple replicas
-gives an effective enhancement in the timescale spanned by the
-multiple simulations, while waiting for an event to occur.
-
-How different than PRD.
+accelerated by application of a bias potential to one or more pairs of
+nearby atoms in the system.  This command can be used to run both
+global and local hyperdyamics.  In global HD a single bond within the
+system is biased on each timestep.  In local HD multiple bonds
+(separated by a sufficient distance) can be biased simultaneously at
+each timestep.  In the bond-boost hyperdynamics context, a "bond" is
+not a covalent bond between a pair of atoms in a molecule.  Rather it
+is simply a pair of nearby atoms as discussed below.
+
+Both global and local HD are described in "(Voter2013)"_#Voter2013 by
+Art Voter and collaborators.  Similar to parallel replica dynamics
+(PRD), global and local HD are methods for performing accelerated
+dynamics that are suitable for infrequent-event systems that obey
+first-order kinetics.  A good overview of accelerated dynamics methods
+for such systems in given in "(Voter2002)"_#Voter2002hd from the same
+group.  To quote from the review paper: "The dynamical evolution is
+characterized by vibrational excursions within a potential basin,
+punctuated by occasional transitions between basins."  The transition
+probability is characterized by p(t) = k*exp(-kt) where k is the rate
+constant.  Running multiple replicas gives an effective enhancement in
+the timescale spanned by the multiple simulations, while waiting for
+an event to occur.
+
+Both HD and PRD produce a time-accurate trajectory that effectively
+extends the timescale over which a system can be simulated, but they
+do it differently.  HD uses a single replica of the system and
+accelerates time by biasing the interaction potential in a manner such
+that each timestep is effectively longer.  PRD creates Nr replicas of
+the system and runs dynamics on each independently with a normal
+unbiased potential until an event occurs in one of the replicas.  The
+time between events is reduced by a factor of Nr replicas.  For both
+methods, per wall-clock second, more physical time elapses and more
+events occur.  See the "prd"_prd.html doc page for more info about
+PRD.
 
 An HD run has several stages, which are repeated each time an "event"
-occurs, as defined below.  The logic for a HD run is as follows:
+occurs, as explained below.  The logic for an HD run is as follows:
 
 quench
-reset list of bonds :pre
+create initial list of bonds :pre
 
 while (time remains):
-  run dynamics for Nevery steps
+  run dynamics for Nevent steps
   quench
   check for an event
   if event occurred: reset list of bonds
   restore pre-quench state :pre
 
-Explain each of the steps above.  Explain what list of bonds is and
-how event is detected.
+The list of bonds is the list of atom pairs of atoms that are within a
+short cutoff distance of each other after the system energy is
+minimized (quenched).  This list is created and reset by a "fix
+hyper/global"_fix_hyper_global.html or "fix
+hyper/local"_fix_hyper_local.html command specified as {fix-ID}.  At
+every dynamics timestep, the same fix selects one of more bonds to
+apply a bias potential to.
+
+IMPORTANT NOTE: The style of fix associated with the specified
+{fix-ID} determines whether you are running the global versus local
+hyperdynamics algorithm.
+
+Dynamics (with the bias potential) is run continuously, stopping every
+{Nevent} steps to check if a transition event has occurred.  The
+specified {N} for total steps must be a multiple of {Nevent}.  check
+is performed by quenching the system and comparing the resulting atom
+coordinates to the coordinates from the previous basin.
+
+A quench is an energy minimization and is performed by whichever
+algorithm has been defined by the "min_style"_min_style.html command.
+Minimization parameters may be set via the
+"min_modify"_min_modify.html command and by the {min} keyword of the
+hyper command.  The latter are the settings that would be used with
+the "minimize"_minimize.html command.  Note that typically, you do not
+need to perform a highly-converged minimization to detect a transition
+event, though you may need to in order to prevent a set of atoms in
+the system from relaxing to a saddle point.
+
+The event check is performed by a compute with the specified
+{compute-ID}.  Currently there is only one compute that works with the
+hyper command, which is the "compute
+event/displace"_compute_event_displace.html command.  Other
+event-checking computes may be added.  "Compute
+event/displace"_compute_event_displace.html checks whether any atom in
+the compute group has moved further than a specified threshold
+distance.  If so, an "event" has occurred.
+
+If this happens, the list of bonds is reset, since some bond pairs
+are likely now too far apart, and new pairs are likely close enough
+to be considered a bond.  The pre-quenched state of the
+system (coordinates and velocities) is restored, and dynamics continue.
+
+At the end of the hyper run, a variety of statistics are output to the
+screen and logfile.  These include info relevant to both global and
+local hyperdynamics, such as the number of events and the elapsed
+hyper time (acclerated time), And it includes info specific to one or
+the other, depending on which style of fix was specified by {fix-ID}.
 
-Explain the specified fix and compute.
+:line
 
-Statistics about the number of events, the number of bonds that a bias
-potential is applied to, the accumulated hyper time, etc are stored by
-the fix and can be output with thermodynamic output.
+The optional keywords operate as follows.
+
+As explained above, the {min} keyword can be used to specify
+parameters for the quench.  Their meaning is the same
+as for the "minimize"_minimize.html command
+
+The {time} keyword determines how the {N} timesteps argument is
+interpreted.  If {time} is set to {steps}, then hyperdynamics is run
+for {N} timesteps.  If the time acceleration provided by the bias
+potential is 10x, then that is the same as running a normal (unbiased)
+MD simulation for 10N steps.  If {time} is set to {clock}, then
+hyperdynamics would run for N/10 steps, since that would correspond to
+an elapsed real time of N*dt.
+
+The {dump} keyword can be used to trigger a specific dump command with
+the specified {dump-ID} to output a snapshot each time an event is
+detected.  It can be specified multiple times with different {dump-ID}
+values, as in the example above.  These snapshots will be for the
+quenched state of the system on a timestep that is a multiple of
+{Nevent}, i.e. a timestep after the event has occurred.  Note that any
+dump command in the input script will also output snapshots at
+whatever timestep interval it defines via its {N} argument; see the
+"dump"_dump.html command for details.  This means if you only want a
+particular dump to output snapshots when events are detected, you
+should specify its {N} as a value larger than the length of the
+hyperdynamics run.
+
+As in the code logic above, the bond list is normally only reset when
+an event occurs.  The {rebond} keyword will force a reset of the bond
+list every {Nrebond} steps, even if an event has not occurred.
+{Nrebond} must be a multiple of {Nevent}.  This can be useful to check
+if more frequent resets alter event statistics, perhaps because the
+parameters chosen for defining what is a bond and what is an event are
+producing bad dynamics in the presence of the bias potential.
 
 :line
 
 [Restrictions:]
 
 This command can only be used if LAMMPS was built with the REPLICA
-package.  See the "Making LAMMPS"_Section_start.html#start_3 section
-for more info on packages.
-
-NOTE: is this true?
-This command cannot be used when any fixes are defined that keep track
-of elapsed time to perform time-dependent operations.  Examples
-include the "ave" fixes such as "fix ave/chunk"_fix_ave_chunk.html.
-Also "fix dt/reset"_fix_dt_reset.html and "fix
-deposit"_fix_deposit.html.
+package.  See the "Build package"_Build_package.html doc
+page for more info.
 
 [Related commands:]
 
 "fix hyper/global"_fix_hyper_global.html, "fix
 hyper/local"_fix_hyper_local.html, "compute
-event/displace"_compute_event_displace.html
+event/displace"_compute_event_displace.html, "prd"_prd.html
 
 [Default:]