Merge pull request #1421 from julient31/gneb_spin

Adding a GNEB implementation in the SPIN package

doc/src/Commands_all.txt

@@ -83,6 +83,7 @@ An alphabetic list of all general LAMMPS commands.
 "molecule"_molecule.html,
 "ndx2group"_group2ndx.html,
 "neb"_neb.html,
+"neb_spin"_neb_spin.html,
 "neigh_modify"_neigh_modify.html,
 "neighbor"_neighbor.html,
 "newton"_newton.html,

doc/src/Commands_category.txt

@@ -116,6 +116,7 @@ Actions:
 
 "minimize"_minimize.html,
 "neb"_neb.html,
+"neb_spin"_neb_spin.html,
 "prd"_prd.html,
 "rerun"_rerun.html,
 "run"_run.html,

doc/src/Commands_fix.txt

@@ -107,6 +107,7 @@ OPT.
 "mvv/edpd"_fix_mvv_dpd.html,
 "mvv/tdpd"_fix_mvv_dpd.html,
 "neb"_fix_neb.html,
+"neb_spin"_fix_neb_spin.html,
 "nph (ko)"_fix_nh.html,
 "nph/asphere (o)"_fix_nph_asphere.html,
 "nph/body"_fix_nph_body.html,

doc/src/Eqs/neb_spin_angle.tex

@@ -0,0 +1,15 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath, amssymb, graphics, setspace}
+
+\begin{document}
+\begin{varwidth}{50in}
+  \begin{equation}
+    \omega_i^{\nu} = 
+    (\nu - 1) \Delta \omega_i
+    {\rm ~~and~~} \Delta \omega_i = \frac{\omega_i}{Q-1}
+  , \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}

doc/src/Eqs/neb_spin_k.tex

@@ -0,0 +1,16 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath, amssymb, graphics, setspace}
+
+\begin{document}
+\begin{varwidth}{50in}
+  \begin{equation}
+    \vec{k}_i =  
+    \frac{\vec{m}_i^I \times \vec{m}_i^F}{\left|\vec{m}_i^I 
+    \times \vec{m}_i^F\right|} 
+    %&{\rm ~if~}& \vec{m}_i^I \times \vec{m}_i^F  
+  , \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}

doc/src/Eqs/neb_spin_rodrigues_formula.tex

@@ -0,0 +1,16 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath, amssymb, graphics, setspace}
+
+\begin{document}
+\begin{varwidth}{50in}
+  \begin{equation}
+    \vec{m}_i^{\nu} = \vec{m}_i^{I} \cos(\omega_i^{\nu})
+    + (\vec{k}_i \times \vec{m}_i^{I}) \sin(\omega_i^{\nu})
+    + (1.0-\cos(\omega_i^{\nu})) \vec{k}_i (\vec{k}_i\cdot
+    \vec{m}_i^{I})
+  , \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}

doc/src/Howto_replica.txt

@@ -17,6 +17,7 @@ periodically.
 These are the relevant commands:
 
 "neb"_neb.html for nudged elastic band calculations
+"neb_spin"_neb_spin.html for magnetic nudged elastic band calculations
 "prd"_prd.html for parallel replica dynamics
 "tad"_tad.html for temperature accelerated dynamics
 "temper"_temper.html for parallel tempering

doc/src/Packages_details.txt

@@ -918,6 +918,7 @@ src/SPIN: filenames -> commands
 "fix nve/spin"_fix_nve_spin.html
 "fix precession/spin"_fix_precession_spin.html
 "compute spin"_compute_spin.html
+"neb/spin"_neb_spin.html
 examples/SPIN :ul
 
 :line

doc/src/commands_list.txt

@@ -67,6 +67,7 @@ Commands :h1
    minimize
    molecule
    neb
+   neb_spin
    neigh_modify
    neighbor
    newton

doc/src/fix_neb.txt

@@ -97,7 +97,7 @@ Note that in this case the specified {Kspring} is in force/distance
 units.
 
 With a value of {ideal}, the spring force is computed as suggested in
-"(WeinenE)"_#WeinenE :
+"(WeinanE)"_#WeinanE :
 
 Fnudge_parallel = -{Kspring} * (RD-RDideal) / (2 * meanDist) :pre
 
@@ -224,8 +224,8 @@ specified (no inter-replica force on the end replicas).
 [(Henkelman2)] Henkelman, Uberuaga, Jonsson, J Chem Phys, 113,
 9901-9904 (2000).
 
-:link(WeinenE)
-[(WeinenE)] E, Ren, Vanden-Eijnden, Phys Rev B, 66, 052301 (2002).
+:link(WeinanE)
+[(WeinanE)] E, Ren, Vanden-Eijnden, Phys Rev B, 66, 052301 (2002).
 
 :link(Jonsson)
 [(Jonsson)] Jonsson, Mills and Jacobsen, in Classical and Quantum

doc/src/fix_neb_spin.txt

@@ -0,0 +1,76 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Commands_all.html)
+
+:line
+
+fix neb/spin command :h3
+
+[Syntax:]
+
+fix ID group-ID neb/spin Kspring :pre
+
+ID, group-ID are documented in "fix"_fix.html command :ulb,l
+neb/spin = style name of this fix command :l
+Kspring = spring constant for parallel nudging force
+(force/distance units or force units, see parallel keyword) :pre,ule
+
+[Examples:]
+
+fix 1 active neb/spin 1.0
+
+[Description:]
+
+Add nudging forces to spins in the group for a multi-replica
+simulation run via the "neb/spin"_neb_spin.html command to perform a 
+geodesic nudged elastic band (GNEB) calculation for finding the 
+transition state.
+Hi-level explanations of GNEB are given with the 
+"neb/spin"_neb_spin.html command and on the 
+"Howto replica"_Howto_replica.html doc page.  
+The fix neb/spin command must be used with the "neb/spin" command and 
+defines how inter-replica nudging forces are computed.  A GNEB 
+calculation is divided in two stages. In the first stage n replicas 
+are relaxed toward a MEP until convergence.  In the second stage, the 
+climbing image scheme is enabled, so that the replica having the highest 
+energy relaxes toward the saddle point (i.e. the point of highest energy 
+along the MEP), and a second relaxation is performed.
+
+The nudging forces are calculated as explained in
+"(BessarabB)"_#BessarabB).
+See this reference for more explanation about their expression.
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+No information about this fix is written to "binary restart
+files"_restart.html.  None of the "fix_modify"_fix_modify.html options
+are relevant to this fix.  No global or per-atom quantities are stored
+by this fix for access by various "output commands"_Howto_output.html.
+No parameter of this fix can be used with the {start/stop} keywords of
+the "run"_run.html command.
+
+The forces due to this fix are imposed during an energy minimization,
+as invoked by the "minimize"_minimize.html command via the
+"neb/spin"_neb_spin.html command.
+
+[Restrictions:]
+
+This command can only be used if LAMMPS was built with the SPIN
+package.  See the "Build package"_Build_package.html doc
+page for more info.
+
+[Related commands:]
+
+"neb_spin"_neb_spin.html
+
+[Default:]
+
+none
+
+:line
+
+:link(BessarabB)
+[(BessarabB)] Bessarab, Uzdin, Jonsson, Comp Phys Comm, 196,
+335-347 (2015).

doc/src/fix_precession_spin.txt

@@ -31,7 +31,7 @@ fix 1 all precession/spin zeeman 0.1 0.0 0.0 1.0 anisotropy 0.001 0.0 0.0 1.0 :p
 
 [Description:]
 
-Impose a force torque to each magnetic spin in the group.
+This fix applies a precession torque to each magnetic spin in the group.
 
 Style {zeeman} is used for the simulation of the interaction
 between the magnetic spins in the defined group and an external

doc/src/fixes.txt

@@ -84,6 +84,7 @@ Fixes :h1
    fix_msst
    fix_mvv_dpd
    fix_neb
+   fix_neb_spin
    fix_nh
    fix_nh_eff
    fix_nh_uef

doc/src/lammps.book

@@ -179,6 +179,7 @@ min_spin.html
 minimize.html
 molecule.html
 neb.html
+neb_spin.html
 neigh_modify.html
 neighbor.html
 newton.html
@@ -309,6 +310,7 @@ fix_mscg.html
 fix_msst.html
 fix_mvv_dpd.html
 fix_neb.html
+fix_neb_spin.html
 fix_nh.html
 fix_nh_eff.html
 fix_nph_asphere.html

doc/src/neb_spin.txt

@@ -0,0 +1,375 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Commands_all.html)
+
+:line
+
+neb command :h3
+
+[Syntax:]
+
+neb/spin etol ttol N1 N2 Nevery file-style arg keyword :pre
+
+etol = stopping tolerance for energy (energy units) :ulb,l
+ttol = stopping tolerance for torque ( units) :l
+N1 = max # of iterations (timesteps) to run initial NEB :l
+N2 = max # of iterations (timesteps) to run barrier-climbing NEB :l
+Nevery = print replica energies and reaction coordinates every this many timesteps :l
+file-style = {final} or {each} or {none} :l
+  {final} arg = filename
+    filename = file with initial coords for final replica
+      coords for intermediate replicas are linearly interpolated
+      between first and last replica
+  {each} arg = filename
+    filename = unique filename for each replica (except first)
+      with its initial coords
+  {none} arg = no argument all replicas assumed to already have
+      their initial coords :pre
+keyword = {verbose}
+:ule
+
+[Examples:]
+
+neb/spin 0.1 0.0 1000 500 50 final coords.final
+neb/spin 0.0 0.001 1000 500 50 each coords.initial.$i
+neb/spin 0.0 0.001 1000 500 50 none verbose :pre
+
+[Description:]
+
+Perform a geodesic nudged elastic band (GNEB) calculation using multiple
+replicas of a system.  Two or more replicas must be used; the first
+and last are the end points of the transition path.
+
+GNEB is a method for finding both the spin configurations and height
+of the energy barrier associated with a transition state, e.g.
+spins to perform a collective rotation from one energy basin to
+another. 
+The implementation in LAMMPS follows the discussion in the
+following paper: "(BessarabA)"_#BessarabA.
+
+Each replica runs on a partition of one or more processors.  Processor
+partitions are defined at run-time using the "-partition command-line
+switch"_Run_options.html.  Note that if you have MPI installed, you
+can run a multi-replica simulation with more replicas (partitions)
+than you have physical processors, e.g you can run a 10-replica
+simulation on just one or two processors.  You will simply not get the
+performance speed-up you would see with one or more physical
+processors per replica.  See the "Howto replica"_Howto_replica.html
+doc page for further discussion.
+
+NOTE: As explained below, a GNEB calculation performs a damped dynamics
+minimization across all the replicas.  The "spin"_min_spin.html
+style minimizer has to be defined in your input script. 
+
+When a GNEB calculation is performed, it is assumed that each replica
+is running the same system, though LAMMPS does not check for this.
+I.e. the simulation domain, the number of magnetic atoms, the 
+interaction potentials, and the starting configuration when the neb 
+command is issued should be the same for every replica.
+
+In a GNEB calculation each replica is connected to other replicas by
+inter-replica nudging forces.  These forces are imposed by the "fix
+neb/spin"_fix_neb_spin.html command, which must be used in conjunction 
+with the neb command.  
+The group used to define the fix neb/spin command defines the
+GNEB magnetic atoms which are the only ones that inter-replica springs 
+are applied to.  
+If the group does not include all magnetic atoms, then non-GNEB
+magnetic atoms have no inter-replica springs and the torques they feel 
+and their precession motion is computed in the usual way due only 
+to other magnetic atoms within their replica.  
+Conceptually, the non-GNEB atoms provide a background force field for 
+the GNEB atoms.  
+Their magnetic spins can be allowed to evolve during the GNEB 
+minimization procedure.
+
+The initial spin configuration for each of the replicas can be
+specified in different manners via the {file-style} setting, as
+discussed below.  Only atomic spins whose initial coordinates should 
+differ from the current configuration need to be specified.
+
+Conceptually, the initial and final configurations for the first
+replica should be states on either side of an energy barrier.
+
+As explained below, the initial configurations of intermediate
+replicas can be spin coordinates interpolated in a linear fashion
+between the first and last replicas.  This is often adequate for
+simple transitions.  For more complex transitions, it may lead to slow
+convergence or even bad results if the minimum energy path (MEP, see
+below) of states over the barrier cannot be correctly converged to
+from such an initial path.  In this case, you will want to generate
+initial states for the intermediate replicas that are geometrically
+closer to the MEP and read them in.
+
+:line
+
+For a {file-style} setting of {final}, a filename is specified which
+contains atomic and spin coordinates for zero or more atoms, in the 
+format described below.  
+For each atom that appears in the file, the new coordinates are 
+assigned to that atom in the final replica.  Each intermediate replica 
+also assigns a new spin to that atom in an interpolated manner.  
+This is done by using the current direction of the spin at the starting 
+point and the read-in direction as the final point.  
+The "angular distance" between them is calculated, and the new direction 
+is assigned to be a fraction of the angular distance.
+
+NOTE: The "angular distance" between the starting and final point is 
+evaluated in the geodesic sense, as described in 
+"(BessarabA)"_#BessarabA. 
+
+NOTE: The angular interpolation between the starting and final point 
+is achieved using Rodrigues formula:
+
+:c,image(Eqs/neb_spin_rodrigues_formula.jpg)
+
+where m_i^I is the initial spin configuration for the spin i,
+omega_i^nu is a rotation angle defined as:
+
+:c,image(Eqs/neb_spin_angle.jpg)
+
+with nu the image number, Q the total number of images, and 
+omega_i the total rotation between the initial and final spins.
+k_i defines a rotation axis such as:
+
+:c,image(Eqs/neb_spin_k.jpg)
+
+if the initial and final spins are not aligned.
+If the initial and final spins are aligned, then their cross
+product is null, and the expression above does not apply.
+If they point toward the same direction, the intermediate images 
+conserve the same orientation.
+If the initial and final spins are aligned, but point toward
+opposite directions, an arbitrary rotation vector belonging to
+the plane perpendicular to initial and final spins is chosen. 
+In this case, a warning message is displayed.
+
+For a {file-style} setting of {each}, a filename is specified which is
+assumed to be unique to each replica. 
+See the "neb"_neb.html documentation page for more information about this 
+option.
+
+For a {file-style} setting of {none}, no filename is specified.  Each
+replica is assumed to already be in its initial configuration at the
+time the neb command is issued.  This allows each replica to define
+its own configuration by reading a replica-specific data or restart or
+dump file, via the "read_data"_read_data.html,
+"read_restart"_read_restart.html, or "read_dump"_read_dump.html
+commands.  The replica-specific names of these files can be specified
+as in the discussion above for the {each} file-style.  Also see the
+section below for how a NEB calculation can produce restart files, so
+that a long calculation can be restarted if needed.
+
+NOTE: None of the {file-style} settings change the initial
+configuration of any atom in the first replica.  The first replica
+must thus be in the correct initial configuration at the time the neb
+command is issued.
+
+:line
+
+A NEB calculation proceeds in two stages, each of which is a
+minimization procedure, performed via damped dynamics.  To enable
+this, you must first define a damped spin dynamics
+"min_style"_min_style.html, using the {spin} style (see
+"min_spin"_min_spin.html for more information).  
+The other styles cannot be used, since they relax the lattice
+degrees of freedom instead of the spins.
+
+The minimizer tolerances for energy and force are set by {etol} and
+{ttol}, the same as for the "minimize"_minimize.html command.
+
+A non-zero {etol} means that the GNEB calculation will terminate if the
+energy criterion is met by every replica.  The energies being compared
+to {etol} do not include any contribution from the inter-replica
+nudging forces, since these are non-conservative.  A non-zero {ttol}
+means that the GNEB calculation will terminate if the torque criterion
+is met by every replica.  The torques being compared to {ttol} include
+the inter-replica nudging forces.
+
+The maximum number of iterations in each stage is set by {N1} and
+{N2}.  These are effectively timestep counts since each iteration of
+damped dynamics is like a single timestep in a dynamics
+"run"_run.html.  During both stages, the potential energy of each
+replica and its normalized distance along the reaction path (reaction
+coordinate RD) will be printed to the screen and log file every
+{Nevery} timesteps.  The RD is 0 and 1 for the first and last replica.
+For intermediate replicas, it is the cumulative angular distance 
+(normalized by the total cumulative angular distance) between adjacent 
+replicas, where "distance" is defined as the length of the 3N-vector of 
+the geodesic distances in spin coordinates, with N the number of
+GNEB spins involved (see equation (13) in "(BessarabA)"_#BessarabA).
+These outputs allow you to monitor NEB's progress in
+finding a good energy barrier.  {N1} and {N2} must both be multiples
+of {Nevery}.
+
+In the first stage of GNEB, the set of replicas should converge toward
+a minimum energy path (MEP) of conformational states that transition
+over a barrier.  The MEP for a transition is defined as a sequence of
+3N-dimensional spin states, each of which has a potential energy 
+gradient parallel to the MEP itself.  
+The configuration of highest energy along a MEP corresponds to a saddle 
+point.  The replica states will also be roughly equally spaced along 
+the MEP due to the inter-replica nudging force added by the 
+"fix neb"_fix_neb.html command.
+
+In the second stage of GNEB, the replica with the highest energy is
+selected and the inter-replica forces on it are converted to a force
+that drives its spin coordinates to the top or saddle point of the
+barrier, via the barrier-climbing calculation described in
+"(BessarabA)"_#BessarabA.  As before, the other replicas rearrange
+themselves along the MEP so as to be roughly equally spaced.
+
+When both stages are complete, if the GNEB calculation was successful,
+the configurations of the replicas should be along (close to) the MEP
+and the replica with the highest energy should be a spin
+configuration at (close to) the saddle point of the transition. The
+potential energies for the set of replicas represents the energy
+profile of the transition along the MEP.
+
+:line
+
+An atom map must be defined which it is not by default for "atom_style
+atomic"_atom_style.html problems.  The "atom_modify
+map"_atom_modify.html command can be used to do this.
+
+An initial value can be defined for the timestep. Although, the {spin} 
+minimization algorithm is an adaptive timestep methodology, so that 
+this timestep is likely to evolve during the calculation.  
+
+The minimizers in LAMMPS operate on all spins in your system, even
+non-GNEB atoms, as defined above.  
+
+:line
+
+Each file read by the neb/spin command containing spin coordinates used
+to initialize one or more replicas must be formatted as follows.
+
+The file can be ASCII text or a gzipped text file (detected by a .gz
+suffix).  The file can contain initial blank lines or comment lines
+starting with "#" which are ignored.  The first non-blank, non-comment
+line should list N = the number of lines to follow.  The N successive
+lines contain the following information:
+
+ID1 g1 x1 y1 z1 sx1 sy1 sz1
+ID2 g2 x2 y2 z2 sx2 sy2 sz2
+...
+IDN gN yN zN sxN syN szN :pre
+
+The fields are the atom ID, the norm of the associated magnetic spin, 
+followed by the {x,y,z} coordinates and the {sx,sy,sz} spin coordinates.
+The lines can be listed in any order.  Additional trailing information on
+the line is OK, such as a comment.
+
+Note that for a typical GNEB calculation you do not need to specify
+initial spin coordinates for very many atoms to produce differing starting
+and final replicas whose intermediate replicas will converge to the
+energy barrier.  Typically only new spin coordinates for atoms
+geometrically near the barrier need be specified.
+
+Also note there is no requirement that the atoms in the file
+correspond to the GNEB atoms in the group defined by the "fix
+neb"_fix_neb.html command.  Not every GNEB atom need be in the file,
+and non-GNEB atoms can be listed in the file.
+
+:line
+
+Four kinds of output can be generated during a GNEB calculation: energy
+barrier statistics, thermodynamic output by each replica, dump files,
+and restart files.
+
+When running with multiple partitions (each of which is a replica in
+this case), the print-out to the screen and master log.lammps file
+contains a line of output, printed once every {Nevery} timesteps.  It
+contains the timestep, the maximum torque per replica, the maximum
+torque per atom (in any replica), potential gradients in the initial,
+final, and climbing replicas, the forward and backward energy
+barriers, the total reaction coordinate (RDT), and the normalized
+reaction coordinate and potential energy of each replica.
+
+The "maximum torque per replica" is the two-norm of the
+3N-length vector given by the cross product of a spin by its
+precession vector omega, in each replica, maximized across replicas, 
+which is what the {ttol} setting is checking against.  In this case, N is
+all the atoms in each replica.  The "maximum torque per atom" is the
+maximum torque component of any atom in any replica.  The potential
+gradients are the two-norm of the 3N-length magnetic precession vector 
+solely due to the interaction potential i.e. without adding in 
+inter-replica forces, and projected along the path tangent (as detailed 
+in Appendix D of "(BessarabA)"_#BessarabA).
+
+The "reaction coordinate" (RD) for each replica is the two-norm of the
+3N-length vector of geodesic distances between its spins and the preceding
+replica's spins (see equation (13) of "(BessarabA)"_#BessarabA), added to 
+the RD of the preceding replica. The RD of the first replica RD1 = 0.0; 
+the RD of the final replica RDN = RDT, the total reaction coordinate.  
+The normalized RDs are divided by RDT, so that they form a monotonically 
+increasing sequence from zero to one. When computing RD, N only includes 
+the spins being operated on by the fix neb/spin command.
+
+The forward (reverse) energy barrier is the potential energy of the
+highest replica minus the energy of the first (last) replica.
+
+Supplementary information for all replicas can be printed out to the
+screen and master log.lammps file by adding the verbose keyword. This
+information include the following. 
+The "GradVidottan" are the projections of the potential gradient for 
+the replica i on its tangent vector (as detailed in Appendix D of 
+"(BessarabA)"_#BessarabA).
+The "DNi" are the non normalized geodesic distances (see equation (13) 
+of "(BessarabA)"_#BessarabA), between a replica i and the next replica 
+i+1. For the last replica, this distance is not defined and a "NAN"
+value is the corresponding output. 
+
+When a NEB calculation does not converge properly, the supplementary
+information can help understanding what is going wrong. 
+
+When running on multiple partitions, LAMMPS produces additional log
+files for each partition, e.g. log.lammps.0, log.lammps.1, etc.  For a
+GNEB calculation, these contain the thermodynamic output for each
+replica.
+
+If "dump"_dump.html commands in the input script define a filename
+that includes a {universe} or {uloop} style "variable"_variable.html,
+then one dump file (per dump command) will be created for each
+replica.  At the end of the GNEB calculation, the final snapshot in
+each file will contain the sequence of snapshots that transition the
+system over the energy barrier.  Earlier snapshots will show the
+convergence of the replicas to the MEP.
+
+Likewise, "restart"_restart.html filenames can be specified with a
+{universe} or {uloop} style "variable"_variable.html, to generate
+restart files for each replica.  These may be useful if the GNEB
+calculation fails to converge properly to the MEP, and you wish to
+restart the calculation from an intermediate point with altered
+parameters.
+
+A c file script in provided in the tool/spin/interpolate_gneb
+directory, that interpolates the MEP given the information provided 
+by the verbose output option (as detailed in Appendix D of
+"(BessarabA)"_#BessarabA). 
+
+:line
+
+[Restrictions:]
+
+This command can only be used if LAMMPS was built with the SPIN
+package.  See the "Build package"_Build_package.html doc
+page for more info.
+
+:line
+
+[Related commands:]
+
+"min/spin"_min_spin.html, "fix neb/spin"_fix_neb_spin.html
+
+[Default:]
+
+none
+
+:line
+
+:link(BessarabA)
+[(BessarabA)] Bessarab, Uzdin, Jonsson, Comp Phys Comm, 196,
+335-347 (2015).