Commit JT 032219

- fixed bug 1: precession_spin had no min_setup
- fixed bug 2: incorrect init of spins in neb/spin
- improved doc min_spin.txt (added eqs, and connected to related
files).

doc/src/Eqs/min_spin_damping.tex

@@ -0,0 +1,13 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath, amssymb, graphics, setspace}
+
+\begin{document}
+\begin{varwidth}{50in}
+  \begin{equation}
+    \frac{d \vec{s}_{i}}{dt} = \lambda\, \vec{s}_{i} \times\left( \vec{\omega}_{i} \times\vec{s}_{i} \right)
+    \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}

doc/src/Eqs/min_spin_timestep.tex

@@ -0,0 +1,14 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath, amssymb, graphics, setspace}
+
+\begin{document}
+\begin{varwidth}{50in}
+  \begin{equation}
+    {\Delta t}_{\rm max} = \frac{2\pi}{\kappa
+    \left|\vec{\omega}_{\rm max} \right|} 
+    \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}

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/lammps.book

@@ -174,6 +174,7 @@ mass.html
 message.html
 min_modify.html
 min_style.html
+min_spin.html
 minimize.html
 molecule.html
 neb.html

doc/src/min_modify.txt

@@ -21,7 +21,7 @@ keyword = {dmax} or {line} or {alpha_damp} or {discret_factor}
   {alpha_damp} value = damping
     damping = fictitious Gilbert damping for spin minimization (adim)
   {discret_factor} value = factor
-    factor = defines a dividing factor for adaptive spin timestep (adim) :pre
+    factor = discretization factor for adaptive spin timestep (adim) :pre
 :ule
 
 [Examples:]
@@ -70,7 +70,14 @@ that difference may be smaller than machine epsilon even if atoms
 could move in the gradient direction to reduce forces further.
 
 Keywords {alpha_damp} and {discret_factor} only make sense when
-a {spinmin} minimization style is declared. 
+a "min_spin"_min_spin.html command is declared. 
+Keyword {alpha_damp} defines an analog of a magnetic Gilbert
+damping. It defines a relaxation rate toward an equilibrium for
+a given magnetic system. 
+Keyword {discret_factor} defines a discretization factor for the
+adaptive timestep used in the {spin} minimization. 
+See "min_spin"_min_spin.html for more information about those
+quantities. 
 Default values are alpha_damp = 1.0 and discret_factor = 10.0.
 
 [Restrictions:] none

doc/src/min_spin.txt

@@ -0,0 +1,63 @@
+"LAMMPS WWW Page"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Commands_all.html)
+
+:line
+
+min_style spin command :h3
+
+[Syntax:]
+
+min_style spin :pre
+
+[Examples:]
+
+min_style 	spin
+min_modify      alpha_damp 1.0 discret_factor 10.0 :pre
+
+[Description:]
+
+Apply a minimization algorithm to use when a "minimize"_minimize.html
+command is performed.
+
+Style {spin} defines a damped spin dynamics with an adaptive
+timestep, according to:
+
+:c,image(Eqs/min_spin_damping.jpg)
+
+with lambda a damping coefficient (similar to a Gilbert damping)
+Lambda can be defined by setting the {alpha_damp} keyword with the 
+"min_modify"_min_modify.html command. 
+
+The minimization procedure solves this equation using an
+adaptive timestep. The value of this timestep is conditionned 
+by the largest precession frequency that has to be solved in the 
+system:
+
+:c,image(Eqs/min_spin_timestep.jpg)
+
+with |omega|_{max} the norm of the largest precession frequency
+in the system (across all processes, and across all replicas if a
+spin/neb calculation is performed). 
+
+Kappa defines a discretization factor {discret_factor} for the 
+definition of this timestep. 
+{discret_factor} can be defined with the "min_modify"_min_modify.html
+command.
+
+NOTE: The {spin} style replaces the force tolerance by a torque
+tolerance. See "minimize"_minimize.html for more explanation. 
+
+[Restrictions:] none
+
+[Related commands:]
+
+"min_style"_min_style.html, "minimize"_minimize.html, 
+"min_modify"_min_modify.html
+
+[Default:]
+
+The option defaults are alpha_damp = 1.0 and discret_factor =
+10.0.

doc/src/min_style.txt

@@ -16,7 +16,7 @@ style = {cg} or {hftn} or {sd} or {quickmin} or {fire} or {spin} :ul
 [Examples:]
 
 min_style cg
-min_style spinmin
+min_style spin
 min_style fire :pre
 
 [Description:]
@@ -62,17 +62,13 @@ the velocity non-parallel to the current force vector.  The velocity
 of each atom is initialized to 0.0 by this style, at the beginning of
 a minimization.
 
-Style {spin} is a damped spin dynamics with a variable 
-timestep as described in "(Tranchida)"_#Tranchida.
+Style {spin} is a damped spin dynamics with an adaptive 
+timestep.
 See the "min/spin"_min_spin.html doc page for more information.
 
 Either the {quickmin} and {fire} styles are useful in the context of
 nudged elastic band (NEB) calculations via the "neb"_neb.html command.
 
-The {spin} style is useful in the context of geodesic nudged 
-elastic band (GNEB) calculations via the "neb/spin"_neb_spin.html
-command.
-
 NOTE: The damped dynamic minimizers use whatever timestep you have
 defined via the "timestep"_timestep.html command.  Often they will
 converge more quickly if you use a timestep about 10x larger than you

doc/src/minimize.txt

@@ -103,6 +103,13 @@ the line search fails because the step distance backtracks to 0.0
 the number of outer iterations or timesteps exceeds {maxiter}
 the number of total force evaluations exceeds {maxeval} :ul
 
+NOTE: the "minimization style"_min_style.html {spin} replaces
+the force tolerance {ftol} by a torque tolerance.
+The minimization procedure stops if the 2-norm (length) of the 
+global torque vector (defined as the cross product between the 
+spins and their precession vectors omega) is less than {ftol}, 
+or if any of the other criteria are met.
+
 NOTE: You can also use the "fix halt"_fix_halt.html command to specify
 a general criterion for exiting a minimization, that is a calculation
 performed on the state of the current system, as defined by an

src/SPIN/atom_vec_spin.cpp

@@ -943,6 +943,10 @@ bigint AtomVecSpin::memory_usage()
   return bytes;
 }
 
+/* ----------------------------------------------------------------------
+   clear all forces (mech and mag)
+------------------------------------------------------------------------- */
+
 void AtomVecSpin::force_clear(int /*n*/, size_t nbytes)
 {
   memset(&atom->f[0][0],0,3*nbytes);

src/SPIN/fix_precession_spin.cpp

@@ -170,7 +170,14 @@ void FixPrecessionSpin::setup(int vflag)
 
 /* ---------------------------------------------------------------------- */
 
-void FixPrecessionSpin::post_force(int /*vflag*/)
+void FixPrecessionSpin::min_setup(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPrecessionSpin::post_force(int vflag)
 {
 
   // update mag field with time (potential improvement)

src/SPIN/fix_precession_spin.h

@@ -33,6 +33,7 @@ class FixPrecessionSpin : public Fix {
   int setmask();
   void init();
   void setup(int);
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/SPIN/min_spin.cpp

@@ -102,7 +102,7 @@ void MinSpin::reset_vectors()
 {
   // atomic dof
 
-  // not really good size => sp is 4N vector
+  // size sp is 4N vector
   nvec = 4 * atom->nlocal;
   if (nvec) spvec = atom->sp[0];
   
@@ -132,7 +132,9 @@ int MinSpin::iterate(int maxiter)
     niter++;
 
     // optimize timestep accross processes / replicas
-    
+    // need a force calculation for timestep optimization
+   
+    energy_force(0);
     dts = evaluate_dt();
    
     // apply damped precessional dynamics to the spins

src/SPIN/neb_spin.cpp

@@ -44,7 +44,6 @@
 #include "timer.h"
 #include "memory.h"
 #include "error.h"
-#include "force.h"
 #include "math_const.h"
 
 using namespace LAMMPS_NS;
@@ -101,14 +100,22 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
     spfinal[1] = buf_final[ii+1];
     spfinal[2] = buf_final[ii+2];
 
-    // circular initialization
-    // a better procedure may be developed
-
-    initial_rotation(spinit,spfinal,fraction);
-
-    sp[i][0] = spfinal[0];
-    sp[i][1] = spfinal[1];
-    sp[i][2] = spfinal[2];
+    // interpolate intermediate spin states
+  
+    if (fraction == 0.0) {
+      sp[i][0] = spinit[0];
+      sp[i][1] = spinit[1];
+      sp[i][2] = spinit[2];
+    } else if (fraction == 1.0) {
+      sp[i][0] = spfinal[0];
+      sp[i][1] = spfinal[1];
+      sp[i][2] = spfinal[2];
+    } else {
+      initial_rotation(spinit,spfinal,fraction);
+      sp[i][0] = spfinal[0];
+      sp[i][1] = spfinal[1];
+      sp[i][2] = spfinal[2];
+    }
 
     ii += 3;
   }
@@ -499,7 +506,7 @@ void NEB_spin::readfile(char *file, int flag)
       for (j = 1; j < nwords; j++)
         values[j] = strtok(NULL," \t\n\r\f");
 
-      // adjust atom coord based on replica fraction
+      // adjust spin coord based on replica fraction
       // for flag = 0, interpolate for intermediate and final replicas
       // for flag = 1, replace existing coord with new coord
       // ignore image flags of final x
@@ -530,19 +537,24 @@ void NEB_spin::readfile(char *file, int flag)
 	  spfinal[0] = spx;
 	  spfinal[1] = spy;
 	  spfinal[2] = spz;
-  
-	  printf("test spinit[0]:%g \n",sp[m][0]);
-
+ 
 	  // interpolate intermediate spin states
 
-	  initial_rotation(spinit,spfinal,fraction);
-  
-	  printf("test spfinal[0]:%g \n",spfinal[0]);
-
-	  sp[m][0] = spfinal[0];
-	  sp[m][1] = spfinal[1];
-	  sp[m][2] = spfinal[2];
 	  sp[m][3] = musp;
+	  if (fraction == 0.0) {
+	    sp[m][0] = spinit[0];
+	    sp[m][1] = spinit[1];
+	    sp[m][2] = spinit[2];
+	  } else if (fraction == 1.0) {
+	    sp[m][0] = spfinal[0];
+	    sp[m][1] = spfinal[1];
+	    sp[m][2] = spfinal[2];
+	  } else {
+	    initial_rotation(spinit,spfinal,fraction);
+	    sp[m][0] = spfinal[0];
+	    sp[m][1] = spfinal[1];
+	    sp[m][2] = spfinal[2];
+	  }
         } else {
           sp[m][3] = musp;
 	  x[m][0] = xx;
@@ -560,8 +572,6 @@ void NEB_spin::readfile(char *file, int flag)
     nread += nchunk;
   }
 
-  printf("test sp[1][2]:%g \n",sp[1][2]);
-
   // check that all atom IDs in file were found by a proc
 
   if (flag == 0) {
@@ -605,6 +615,10 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
   // implementing initial rotation using atan2
   // this may not be a sufficient routine, need more accurate verifications
 
+  // interpolation only for intermediate replica
+
+  if (fraction == 0.0 || fraction == 1.0) return;
+
   // initial, final and inter ang. values 
   
   double itheta,iphi,ftheta,fphi,ktheta,kphi;
@@ -636,7 +650,7 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
   spky = sin(ktheta)*sin(kphi);
   spkz = cos(kphi);
 
-  double knormsq = spkx*spkx+spky*spky+spkz*spkz;
+  double knormsq = spkx*spkx + spky*spky + spkz*spkz;
   if (knormsq != 0.0)
     iknorm = 1.0/sqrt(knormsq);
 
@@ -696,9 +710,6 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
   //spky *= iknorm;
   //spkz *= iknorm;
  
-  printf("init: %g %g %g \n",spix,spiy,spiz);
-  printf("fina: %g %g %g \n",spkx,spky,spkz);
-
   sploc[0] = spkx;
   sploc[1] = spky;
   sploc[2] = spkz;

src/SPIN/pair_spin_dmi.cpp

@@ -432,9 +432,9 @@ void PairSpinDmi::compute_dmi(int i, int j, double eij[3], double fmi[3], double
   dmiy = eij[2]*v_dmx[itype][jtype] - eij[0]*v_dmz[itype][jtype];
   dmiz = eij[0]*v_dmy[itype][jtype] - eij[1]*v_dmx[itype][jtype];
 
-  fmi[0] -= (spj[1]*dmiz - spj[2]*dmiy);
-  fmi[1] -= (spj[2]*dmix - spj[0]*dmiz);
-  fmi[2] -= (spj[0]*dmiy - spj[1]*dmix);
+  fmi[0] -= (dmiy*spj[2] - dmiz*spj[1]);
+  fmi[1] -= (dmiz*spj[0] - dmix*spj[2]);
+  fmi[2] -= (dmix*spj[1] - dmiy*spj[0]);
 }
 
 /* ----------------------------------------------------------------------