Commit2 JT 040819

- finished doc (compiled and spell check)
- compiled with -Wall -Wextra, unused var. removed

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/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
@@ -308,6 +309,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

@@ -47,7 +47,7 @@ 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: "(Bessarab)"_#Bessarab.
+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
@@ -78,11 +78,11 @@ 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 precessional motion is computed in the usual way due only 
+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 precess during the GNEB 
+Their magnetic spins can be allowed to evolve during the GNEB 
 minimization procedure.
 
 The initial spin configuration for each of the replicas can be
@@ -117,7 +117,8 @@ 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 "(Bessarab)"_#Bessarab. 
+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:
@@ -142,7 +143,7 @@ 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 perpandicular to initial and final spins is chosen. 
+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
@@ -177,14 +178,14 @@ 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
-{ftol}, the same as for the "minimize"_minimize.html command.
+{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 {ftol}
+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 {ftol} include
+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
@@ -198,7 +199,7 @@ 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 spinsi involved (see equation (13) in "(Bessarab)"_#Bessarab).
+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}.
@@ -217,7 +218,7 @@ 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
-"(Bessarab)"_#Bessarab.  As before, the other replicas rearrange
+"(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,
@@ -227,37 +228,22 @@ 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
 
-A few other settings in your input script are required or advised to
-perform a NEB calculation.  See the NOTE about the choice of timestep
-at the beginning of this doc page.
-
 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.
 
-The minimizers in LAMMPS operate on all atoms in your system, even
+An initial value can be defined for the timestep. Although, the {spin} 
-non-NEB atoms, as defined above.  To prevent non-NEB atoms from moving
+minimization algorithm is an adaptive timestep methodology, so that 
-during the minimization, you should use the "fix
+this timestep is likely to evolve during the calculation.  
-setforce"_fix_setforce.html command to set the force on each of those
-atoms to 0.0.  This is not required, and may not even be desired in
-some cases, but if those atoms move too far (e.g. because the initial
-state of your system was not well-minimized), it can cause problems
-for the NEB procedure.
 
-The damped dynamics "minimizers"_min_style.html, such as {quickmin}
+The minimizers in LAMMPS operate on all spins in your system, even
-and {fire}), adjust the position and velocity of the atoms via an
+non-GNEB atoms, as defined above.  
-Euler integration step.  Thus you must define an appropriate
-"timestep"_timestep.html to use with NEB.  As mentioned above, NEB
-will often converge more quickly if you use a timestep about 10x
-larger than you would normally use for dynamics simulations.
 
 :line
 
-Each file read by the neb command containing atomic coordinates used
+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
@@ -266,130 +252,103 @@ 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 x1 y1 z1
+ID1 g1 x1 y1 z1 sx1 sy1 sz1
-ID2 x2 y2 z2
+ID2 g2 x2 y2 z2 sx2 sy2 sz2
 ...
-IDN xN yN zN :pre
+IDN gN yN zN sxN syN szN :pre
 
-The fields are the atom ID, followed by the x,y,z coordinates.  The
+The fields are the atom ID, the norm of the associated magnetic spin, 
-lines can be listed in any order.  Additional trailing information on
+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 NEB calculation you do not need to specify
+Note that for a typical GNEB calculation you do not need to specify
-initial coordinates for very many atoms to produce differing starting
+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 coordinates for atoms
+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 NEB atoms in the group defined by the "fix
+correspond to the GNEB atoms in the group defined by the "fix
-neb"_fix_neb.html command.  Not every NEB atom need be in the file,
+neb"_fix_neb.html command.  Not every GNEB atom need be in the file,
-and non-NEB atoms can be listed in the file.
+and non-GNEB atoms can be listed in the file.
 
 :line
 
-Four kinds of output can be generated during a NEB calculation: energy
+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 force per replica, the maximum
+contains the timestep, the maximum torque per replica, the maximum
-force per atom (in any replica), potential gradients in the initial,
+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 force per replica" is the two-norm of the 3N-length force
+The "maximum torque per replica" is the two-norm of the
-vector for the atoms in each replica, maximized across replicas, which
+3N-length vector given by the cross product of a spin by its
-is what the {ftol} setting is checking against.  In this case, N is
+precession vector omega, in each replica, maximized across replicas, 
-all the atoms in each replica.  The "maximum force per atom" is the
+which is what the {ttol} setting is checking against.  In this case, N is
-maximum force component of any atom in any replica.  The potential
+all the atoms in each replica.  The "maximum torque per atom" is the
-gradients are the two-norm of the 3N-length force vector solely due to
+maximum torque component of any atom in any replica.  The potential
-the interaction potential i.e.  without adding in inter-replica
+gradients are the two-norm of the 3N-length magnetic precession vector 
-forces.
+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 distances between its atoms and the preceding
+3N-length vector of geodesic distances between its spins and the preceding
-replica's atoms, added to the RD of the preceding replica. The RD of
+replica's spins (see equation (13) of "(BessarabA)"_#BessarabA), added to 
-the first replica RD1 = 0.0; the RD of the final replica RDN = RDT,
+the RD of the preceding replica. The RD of the first replica RD1 = 0.0; 
-the total reaction coordinate.  The normalized RDs are divided by RDT,
+the RD of the final replica RDN = RDT, the total reaction coordinate.  
-so that they form a monotonically increasing sequence from zero to
+The normalized RDs are divided by RDT, so that they form a monotonically 
-one. When computing RD, N only includes the atoms being operated on by
+increasing sequence from zero to one. When computing RD, N only includes 
-the fix neb command.
+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 "path angle" (pathangle) for
+information include the following. 
-the replica i which is the angle between the 3N-length vectors (Ri-1 -
+The "GradVidottan" are the projections of the potential gradient for 
-Ri) and (Ri+1 - Ri) (where Ri is the atomic coordinates of replica
+the replica i on its tangent vector (as detailed in Appendix D of 
-i). A "path angle" of 180 indicates that replicas i-1, i and i+1 are
+"(BessarabA)"_#BessarabA).
-aligned.  "angletangrad" is the angle between the 3N-length tangent
+The "DNi" are the non normalized geodesic distances (see equation (13) 
-vector and the 3N-length force vector at image i. The tangent vector
+of "(BessarabA)"_#BessarabA), between a replica i and the next replica 
-is calculated as in "(HenkelmanA)"_#HenkelmanA for all intermediate
+i+1. For the last replica, this distance is not defined and a "NAN"
-replicas and at R2 - R1 and RM - RM-1 for the first and last replica,
+value is the corresponding output. 
-respectively.  "anglegrad" is the angle between the 3N-length energy
-gradient vector of replica i and that of replica i+1. It is not
-defined for the final replica and reads nan.  gradV is the norm of the
-energy gradient of image i.  ReplicaForce is the two-norm of the
-3N-length force vector (including nudging forces) for replica i.
-MaxAtomForce is the maximum force component of any atom in replica i.
 
 When a NEB calculation does not converge properly, the supplementary
-information can help understanding what is going wrong. For instance
+information can help understanding what is going wrong. 
-when the path angle becomes acute, the definition of tangent used in
-the NEB calculation is questionable and the NEB cannot may diverge
-"(Maras)"_#Maras2.
-
 
 When running on multiple partitions, LAMMPS produces additional log
 files for each partition, e.g. log.lammps.0, log.lammps.1, etc.  For a
-NEB calculation, these contain the thermodynamic output for each
+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 NEB calculation, the final snapshot in
+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 NEB
+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.
 
-There are 2 Python scripts provided in the tools/python directory,
+A c file script in provided in the examples/SPIN/gneb/interpolate
-neb_combine.py and neb_final.py, which are useful in analyzing output
+directory, that interpolates the MEP given the information provided 
-from a NEB calculation.  Assume a NEB simulation with M replicas, and
+by the verbose output option (as detailed in Appendix D of
-the NEB atoms labeled with a specific atom type.
+"(BessarabA)"_#BessarabA). 
-
-The neb_combine.py script extracts atom coords for the NEB atoms from
-all M dump files and creates a single dump file where each snapshot
-contains the NEB atoms from all the replicas and one copy of non-NEB
-atoms from the first replica (presumed to be identical in other
-replicas).  This can be visualized/animated to see how the NEB atoms
-relax as the NEB calculation proceeds.
-
-The neb_final.py script extracts the final snapshot from each of the M
-dump files to create a single dump file with M snapshots.  This can be
-visualized to watch the system make its transition over the energy
-barrier.
-
-To illustrate, here are images from the final snapshot produced by the
-neb_combine.py script run on the dump files produced by the two
-example input scripts in examples/neb.  Click on them to see a larger
-image.
-
-:image(JPG/hop1_small.jpg,JPG/hop1.jpg)
-:image(JPG/hop2_small.jpg,JPG/hop2.jpg)
 
 :line
 
@@ -411,6 +370,6 @@ none
 
 :line
 
-:link(Bessarab)
+:link(BessarabA)
-[(Bessarab)] Bessarab, Uzdin, Jonsson, Comp Phys Comm, 196,
+[(BessarabA)] Bessarab, Uzdin, Jonsson, Comp Phys Comm, 196,
 335-347 (2015).

src/SPIN/fix_neb_spin.cpp

@@ -310,7 +310,6 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 
   int nlocal = atom->nlocal;
   int *mask = atom->mask;
-  double **x = atom->x;
   double **sp = atom->sp;
   double dot = 0.0;
   double prefactor = 0.0;
@@ -322,7 +321,8 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
   nlen = 0.0;
   double tlen = 0.0;
   double gradnextlen = 0.0;
-  double delndots, delpdots;
+  double delndots = 0.0;
+  double delpdots = 0.0;
 
   dotgrad = gradlen = dotpath = dottangrad = 0.0;
 
@@ -584,10 +584,7 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
 
   if (ireplica == 0 || ireplica == nreplica-1) return ;
 
-  double AngularContr;
   dotpath = dotpath/(plen*nlen);
-  AngularContr = 0.5 *(1+cos(MY_PI * dotpath));
-
 
   for (int i = 0; i < nlocal; i++) 
     if (mask[i] & groupbit) {
@@ -620,7 +617,6 @@ void FixNEB_spin::min_post_force(int /*vflag*/)
           fm[i][2] = 0;
         }
       prefactor =  kspring*(nlen-plen);
-      AngularContr=0;
     }
   }
 

src/SPIN/neb_spin.cpp

@@ -88,7 +88,6 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
          int n2steps_in, int nevery_in, double *buf_init, double *buf_final)
   : Pointers(lmp)
 {
-  double delspx,delspy,delspz;
 
   etol = etol_in;
   ttol = ftol_in;
@@ -450,7 +449,7 @@ void NEB_spin::readfile(char *file, int flag)
   tagint tag;
   char *eof,*start,*next,*buf;
   char line[MAXLINE];
-  double xx,yy,zz,delx,dely,delz;
+  double xx,yy,zz;
   double musp,spx,spy,spz;
 
   if (me_universe == 0 && screen)
@@ -701,6 +700,20 @@ int NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
     }
   }
 
+  // knormsq should not be 0 anymore
+
+  if (knormsq == 0.0)
+    error->all(FLERR,"Incorrect initial rotation operation");
+
+  // normalize k vector
+
+  iknorm = 1.0/sqrt(knormsq);
+  kx *= iknorm;
+  ky *= iknorm;
+  kz *= iknorm;
+
+  // calc. k x spi and total rotation angle
+
   kcrossx = ky*spiz - kz*spiy;
   kcrossy = kz*spix - kx*spiz;
   kcrossz = kx*spiy - ky*spix;
@@ -710,7 +723,7 @@ int NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
   omega = acos(sidotsf);
   omega *= fraction;
 
-  // applying Rodrigues' formula
+  // apply Rodrigues' formula
 
   spkx = spix*cos(omega); 
   spky = spiy*cos(omega); 
@@ -873,7 +886,6 @@ void NEB_spin::print_status()
   }
 
   if (me_universe == 0) {
-    const double todeg=180.0/MY_PI;
     FILE *uscreen = universe->uscreen;
     FILE *ulogfile = universe->ulogfile;
     if (uscreen) {