Commit JT 072619

- added a min_style option for norm type (euclidean or Max)
- adapted and tested spin minimizers
- adapted (net tested) regular minimizers

src/SPIN/min_spin.cpp

@@ -119,7 +119,7 @@ void MinSpin::reset_vectors()
 int MinSpin::iterate(int maxiter)
 {
   bigint ntimestep;
-  double fmdotfm;
+  double fmdotfm,fmsq,fmsqall;
   int flag,flagall;
 
   for (int iter = 0; iter < maxiter; iter++) {
@@ -166,8 +166,20 @@ int MinSpin::iterate(int maxiter)
     // magnetic torque tolerance criterion
     // sync across replicas if running multi-replica minimization
 
+    fmdotfm = fmsq = fmsqall = 0.0;
     if (update->ftol > 0.0) {
-      fmdotfm = max_torque();
+      if (normstyle == 1) {		// max torque norm
+	fmsq = max_torque();
+	fmsqall = fmsq;
+	if (update->multireplica == 0)
+	  MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_MAX,universe->uworld);
+      } else {				// Euclidean torque norm
+	fmsq = total_torque();
+	fmsqall = fmsq;
+	if (update->multireplica == 0)
+	  MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_SUM,universe->uworld);
+      }
+      fmdotfm = fmsqall*fmsqall;
       if (update->multireplica == 0) {
         if (fmdotfm < update->ftol*update->ftol) return FTOL;
       } else {
@@ -297,77 +309,3 @@ void MinSpin::advance_spins(double dts)
     // because no need for simplecticity
   }
 }
-
-/* ----------------------------------------------------------------------
-   compute and return ||mag. torque||_2^2
-------------------------------------------------------------------------- */
-
-double MinSpin::fmnorm_sqr()
-{
-  int nlocal = atom->nlocal;
-  double tx,ty,tz;
-  double **sp = atom->sp;
-  double **fm = atom->fm;
-
-  // calc. magnetic torques
-
-  double local_norm2_sqr = 0.0;
-  for (int i = 0; i < nlocal; i++) {
-    tx = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]);
-    ty = (fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]);
-    tz = (fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]);
-
-    local_norm2_sqr += tx*tx + ty*ty + tz*tz;
-  }
-
-  // no extra atom calc. for spins
-
-  if (nextra_atom)
-    error->all(FLERR,"extra atom option not available yet");
-
-  double norm2_sqr = 0.0;
-  MPI_Allreduce(&local_norm2_sqr,&norm2_sqr,1,MPI_DOUBLE,MPI_SUM,world);
-
-  return norm2_sqr;
-}
-
-/* ----------------------------------------------------------------------
-   compute and return  max_i||mag. torque_i||_2
-------------------------------------------------------------------------- */
-
-double MinSpin::max_torque()
-{
-  double fmsq,fmaxsqone,fmaxsqloc,fmaxsqall;
-  int nlocal = atom->nlocal;
-  double hbar = force->hplanck/MY_2PI;
-  double tx,ty,tz;
-  double **sp = atom->sp;
-  double **fm = atom->fm;
-
-  fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
-  for (int i = 0; i < nlocal; i++) {
-    tx = fm[i][1] * sp[i][2] - fm[i][2] * sp[i][1];
-    ty = fm[i][2] * sp[i][0] - fm[i][0] * sp[i][2];
-    tz = fm[i][0] * sp[i][1] - fm[i][1] * sp[i][0];
-    fmsq = tx * tx + ty * ty + tz * tz;
-    fmaxsqone = MAX(fmaxsqone,fmsq);
-  }
-
-  // finding max fm on this replica
-
-  fmaxsqloc = fmaxsqone;
-  MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);
-
-  // finding max fm over all replicas, if necessary
-  // this communicator would be invalid for multiprocess replicas
-
-  fmaxsqall = fmaxsqloc;
-  if (update->multireplica == 1) {
-    fmaxsqall = fmaxsqloc;
-    MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);
-  }
-
-  // multiply it by hbar so that units are in eV
-
-  return sqrt(fmaxsqall) * hbar;
-}

src/SPIN/min_spin.h

@@ -35,8 +35,6 @@ class MinSpin : public Min {
   int iterate(int);
   double evaluate_dt();
   void advance_spins(double);
-  double fmnorm_sqr();
-  double max_torque();
 
  private:
 

src/SPIN/min_spin_oso_cg.cpp

@@ -29,6 +29,7 @@
 #include "universe.h"
 #include "atom.h"
 #include "citeme.h"
+#include "comm.h"
 #include "force.h"
 #include "update.h"
 #include "output.h"
@@ -99,6 +100,13 @@ void MinSpinOSO_CG::init()
 
   Min::init();
 
+  // warning if line_search combined to gneb
+
+  if ((nreplica >= 1) && (linestyle != 4) && (comm->me == 0))
+    error->warning(FLERR,"Line search incompatible gneb");
+
+  // set back use_line_search to 0 if more than one replica
+  
   if (linestyle == 3 && nreplica == 1){
     use_line_search = 1;
   }
@@ -175,7 +183,7 @@ int MinSpinOSO_CG::iterate(int maxiter)
 {
   int nlocal = atom->nlocal;
   bigint ntimestep;
-  double fmdotfm;
+  double fmdotfm,fmsq,fmsqall;
   int flag, flagall;
   double **sp = atom->sp;
   double der_e_cur_tmp = 0.0;
@@ -261,8 +269,20 @@ int MinSpinOSO_CG::iterate(int maxiter)
     // magnetic torque tolerance criterion
     // sync across replicas if running multi-replica minimization
 
+    fmdotfm = fmsq = fmsqall = 0.0;
     if (update->ftol > 0.0) {
-      fmdotfm = max_torque();
+      if (normstyle == 1) {		// max torque norm
+	fmsq = max_torque();
+	fmsqall = fmsq;
+	if (update->multireplica == 0)
+	  MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_MAX,universe->uworld);
+      } else {				// Euclidean torque norm
+	fmsq = total_torque();
+	fmsqall = fmsq;
+	if (update->multireplica == 0)
+	  MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_SUM,universe->uworld);
+      }
+      fmdotfm = fmsqall*fmsqall;
       if (update->multireplica == 0) {
         if (fmdotfm < update->ftol*update->ftol) return FTOL;
       } else {
@@ -353,6 +373,7 @@ void MinSpinOSO_CG::calc_search_direction()
     MPI_Allreduce(&g2old,&g2old_global,1,MPI_DOUBLE,MPI_SUM,world);
 
     // Sum over all replicas. Good for GNEB.
+    
     if (nreplica > 1) {
       g2 = g2_global * factor;
       g2old = g2old_global * factor;
@@ -361,7 +382,9 @@ void MinSpinOSO_CG::calc_search_direction()
     }
     if (fabs(g2_global) < 1.0e-60) beta = 0.0;
     else beta = g2_global / g2old_global;
+    
     // calculate conjugate direction
+    
     for (int i = 0; i < 3 * nlocal; i++) {
       p_s[i] = (beta * p_s[i] - g_cur[i]) * factor;
       g_old[i] = g_cur[i] * factor;
@@ -379,7 +402,7 @@ void MinSpinOSO_CG::advance_spins()
 {
   int nlocal = atom->nlocal;
   double **sp = atom->sp;
-  double rot_mat[9]; // exponential of matrix made of search direction
+  double rot_mat[9];	// exponential of matrix made of search direction
   double s_new[3];
 
   // loop on all spins on proc.
@@ -394,47 +417,6 @@ void MinSpinOSO_CG::advance_spins()
   }
 }
 
-/* ----------------------------------------------------------------------
-   compute and return  max_i||mag. torque_i||_2
-------------------------------------------------------------------------- */
-
-double MinSpinOSO_CG::max_torque()
-{
-  double fmsq,fmaxsqone,fmaxsqloc,fmaxsqall;
-  int nlocal = atom->nlocal;
-  double factor;
-  double hbar = force->hplanck/MY_2PI;
-
-  if (use_line_search) factor = 1.0;
-  else factor = hbar;
-
-  // finding max fm on this proc.
-
-  fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
-  for (int i = 0; i < nlocal; i++) {
-    fmsq = 0.0;
-    for (int j = 0; j < 3; j++)
-      fmsq += g_cur[3 * i + j] * g_cur[3 * i + j];
-    fmaxsqone = MAX(fmaxsqone,fmsq);
-  }
-
-  // finding max fm on this replica
-
-  fmaxsqloc = fmaxsqone;
-  MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);
-
-  // finding max fm over all replicas, if necessary
-  // this communicator would be invalid for multiprocess replicas
-
-  fmaxsqall = fmaxsqloc;
-  if (update->multireplica == 1) {
-    fmaxsqall = fmaxsqloc;
-    MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);
-  }
-
-  return sqrt(fmaxsqall) * factor;
-}
-
 /* ----------------------------------------------------------------------
   calculate 3x3 matrix exponential using Rodrigues' formula
   (R. Murray, Z. Li, and S. Shankar Sastry,
@@ -456,15 +438,14 @@ void MinSpinOSO_CG::rodrigues_rotation(const double *upp_tr, double *out)
       fabs(upp_tr[1]) < 1.0e-40 &&
       fabs(upp_tr[2]) < 1.0e-40){
 
-  // if upp_tr is zero, return unity matrix
+    // if upp_tr is zero, return unity matrix
-  for(int k = 0; k < 3; k++){
+    
-    for(int m = 0; m < 3; m++){
+    for(int k = 0; k < 3; k++){
-      if (m == k)
+      for(int m = 0; m < 3; m++){
-        out[3 * k + m] = 1.0;
+        if (m == k) out[3 * k + m] = 1.0;
-      else
+        else out[3 * k + m] = 0.0;
-        out[3 * k + m] = 0.0;
+      }
     }
-  }
     return;
   }
 
@@ -512,13 +493,14 @@ void MinSpinOSO_CG::rodrigues_rotation(const double *upp_tr, double *out)
 void MinSpinOSO_CG::vm3(const double *m, const double *v, double *out)
 {
   for(int i = 0; i < 3; i++){
-    //out[i] *= 0.0;
     out[i] = 0.0;
-    for(int j = 0; j < 3; j++)
+    for(int j = 0; j < 3; j++) out[i] += *(m + 3 * j + i) * v[j];
-    out[i] += *(m + 3 * j + i) * v[j];
   }
 }
 
+/* ----------------------------------------------------------------------
+  advance spins
+------------------------------------------------------------------------- */
 
 void MinSpinOSO_CG::make_step(double c, double *energy_and_der)
 {
@@ -586,7 +568,7 @@ int MinSpinOSO_CG::calc_and_make_step(double a, double b, int index)
     }
     return 1;
   }
-  else{
+  else {
     double r,f0,f1,df0,df1;
     r = b - a;
     f0 = eprevious;

src/SPIN/min_spin_oso_cg.h

@@ -25,44 +25,44 @@ MinimizeStyle(spin_oso_cg, MinSpinOSO_CG)
 namespace LAMMPS_NS {
 
 class MinSpinOSO_CG: public Min {
-  public:
+ public:
-    MinSpinOSO_CG(class LAMMPS *);
+  MinSpinOSO_CG(class LAMMPS *);
-    virtual ~MinSpinOSO_CG();
+  virtual ~MinSpinOSO_CG();
-    void init();
+  void init();
-    void setup_style();
+  void setup_style();
-    int modify_param(int, char **);
+  void reset_vectors();
-    void reset_vectors();
+  int modify_param(int, char **);
-    int iterate(int);
+  int iterate(int);
-  private:
-    double dt;        // global timestep
-    double dts;       // spin timestep
-    int ireplica,nreplica; // for neb
-    double *spvec;		// variables for atomic dof, as 1d vector
-    double *fmvec;		// variables for atomic dof, as 1d vector
-    double *g_old;  	// gradient vector at previous step
-    double *g_cur;  	// current gradient vector
-    double *p_s;  		// search direction vector
-    double **sp_copy;   // copy of the spins
-    int local_iter;     // for neb
-    int nlocal_max;		// max value of nlocal (for size of lists)
-    double discrete_factor;	// factor for spin timestep evaluation
 
-    double evaluate_dt();
+ private:
-    void advance_spins();
+  int local_iter;		// for neb
-    void calc_gradient();
+  int nlocal_max;		// max value of nlocal (for size of lists)
-    void calc_search_direction();
+  int use_line_search;		// use line search or not.
-    double maximum_rotation(double *);
+  int ireplica,nreplica;	// for neb
-    void vm3(const double *, const double *, double *);
+  double dt;			// global timestep
-    void rodrigues_rotation(const double *, double *);
+  double dts;			// spin timestep
-    int calc_and_make_step(double, double, int);
+  double discrete_factor;	// factor for spin timestep evaluation
-    int awc(double, double, double, double);
+  double der_e_cur;		// current derivative along search dir.
-    void make_step(double, double *);
+  double der_e_pr;		// previous derivative along search dir.
-    double max_torque();
+  double *spvec;		// variables for atomic dof, as 1d vector
-    double der_e_cur;   // current derivative along search dir.
+  double *fmvec;		// variables for atomic dof, as 1d vector
-    double der_e_pr;    // previous derivative along search dir.
+  double *g_old;  		// gradient vector at previous step
-    int use_line_search; // use line search or not.
+  double *g_cur;  		// current gradient vector
+  double *p_s;  		// search direction vector
+  double **sp_copy;		// copy of the spins
 
-    bigint last_negative;
+  void advance_spins();
+  void calc_gradient();
+  void calc_search_direction();
+  void vm3(const double *, const double *, double *);
+  void rodrigues_rotation(const double *, double *);
+  void make_step(double, double *);
+  int calc_and_make_step(double, double, int);
+  int awc(double, double, double, double);
+  double evaluate_dt();
+  double maximum_rotation(double *);
+
+  bigint last_negative;
 };
 
 }

src/SPIN/min_spin_oso_lbfgs.cpp

@@ -26,9 +26,9 @@
 #include <cstdlib>
 #include <cstring>
 #include "min_spin_oso_lbfgs.h"
-#include "universe.h"
 #include "atom.h"
 #include "citeme.h"
+#include "comm.h"
 #include "force.h"
 #include "update.h"
 #include "output.h"
@@ -107,6 +107,13 @@ void MinSpinOSO_LBFGS::init()
 
   Min::init();
 
+  // warning if line_search combined to gneb
+
+  if ((nreplica >= 1) && (linestyle != 4) && (comm->me == 0))
+    error->warning(FLERR,"Line search incompatible gneb");
+
+  // set back use_line_search to 0 if more than one replica
+
   if (linestyle != 4 && nreplica == 1){
     use_line_search = 1;
   }
@@ -188,7 +195,7 @@ int MinSpinOSO_LBFGS::iterate(int maxiter)
 {
   int nlocal = atom->nlocal;
   bigint ntimestep;
-  double fmdotfm;
+  double fmdotfm,fmsq,fmsqall;
   int flag, flagall;
   double **sp = atom->sp;
   double der_e_cur_tmp = 0.0;
@@ -280,8 +287,20 @@ int MinSpinOSO_LBFGS::iterate(int maxiter)
     // magnetic torque tolerance criterion
     // sync across replicas if running multi-replica minimization
 
+    fmdotfm = fmsq = fmsqall = 0.0;
     if (update->ftol > 0.0) {
-      fmdotfm = max_torque();
+      if (normstyle == 1) {		// max torque norm
+	fmsq = max_torque();
+	fmsqall = fmsq;
+	if (update->multireplica == 0)
+	  MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_MAX,universe->uworld);
+      } else {				// Euclidean torque norm
+	fmsq = total_torque();
+	fmsqall = fmsq;
+	if (update->multireplica == 0)
+	  MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_SUM,universe->uworld);
+      }
+      fmdotfm = fmsqall*fmsqall;
       if (update->multireplica == 0) {
         if (fmdotfm < update->ftol*update->ftol) return FTOL;
       } else {
@@ -534,42 +553,6 @@ void MinSpinOSO_LBFGS::advance_spins()
   }
 }
 
-/* ----------------------------------------------------------------------
-   compute and return  max_i||mag. torque_i||_2
-------------------------------------------------------------------------- */
-
-double MinSpinOSO_LBFGS::max_torque()
-{
-  double fmsq,fmaxsqone,fmaxsqloc,fmaxsqall;
-  int nlocal = atom->nlocal;
-
-  // finding max fm on this proc.
-
-  fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
-  for (int i = 0; i < nlocal; i++) {
-    fmsq = 0.0;
-    for (int j = 0; j < 3; j++)
-      fmsq += g_cur[3 * i + j] * g_cur[3 * i + j];
-    fmaxsqone = MAX(fmaxsqone,fmsq);
-  }
-
-  // finding max fm on this replica
-
-  fmaxsqloc = fmaxsqone;
-  MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);
-
-  // finding max fm over all replicas, if necessary
-  // this communicator would be invalid for multiprocess replicas
-
-  fmaxsqall = fmaxsqloc;
-  if (update->multireplica == 1) {
-    fmaxsqall = fmaxsqloc;
-    MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);
-  }
-
-  return sqrt(fmaxsqall);
-}
-
 /* ----------------------------------------------------------------------
   calculate 3x3 matrix exponential using Rodrigues' formula
   (R. Murray, Z. Li, and S. Shankar Sastry,

src/SPIN/min_spin_oso_lbfgs.h

@@ -25,45 +25,45 @@ MinimizeStyle(spin_oso_lbfgs, MinSpinOSO_LBFGS)
 namespace LAMMPS_NS {
 
 class MinSpinOSO_LBFGS: public Min {
-  public:
+ public:
-    MinSpinOSO_LBFGS(class LAMMPS *);
+  MinSpinOSO_LBFGS(class LAMMPS *);
-    virtual ~MinSpinOSO_LBFGS();
+  virtual ~MinSpinOSO_LBFGS();
-    void init();
+  void init();
-    void setup_style();
+  void setup_style();
-    int modify_param(int, char **);
+  int modify_param(int, char **);
-    void reset_vectors();
+  void reset_vectors();
-    int iterate(int);
+  int iterate(int);
-  private:
-    int ireplica,nreplica; 	// for neb
-    double *spvec;		// variables for atomic dof, as 1d vector
-    double *fmvec;		// variables for atomic dof, as 1d vector
-    double *g_old;  		// gradient vector at previous step
-    double *g_cur;  		// current gradient vector
-    double *p_s;  		// search direction vector
-    int local_iter;     	// for neb
-    int nlocal_max;		// max value of nlocal (for size of lists)
 
-    void advance_spins();
+ private:
-    void calc_gradient();
+  int local_iter;     		// for neb
-    void calc_search_direction();
+  int use_line_search; 		// use line search or not.
-    double maximum_rotation(double *);
+  int nlocal_max;		// max value of nlocal (for size of lists)
-    void vm3(const double *, const double *, double *);
+  int ireplica,nreplica; 	// for neb
-    void rodrigues_rotation(const double *, double *);
+  double der_e_cur;		// current derivative along search dir.
-    int calc_and_make_step(double, double, int);
+  double der_e_pr;    		// previous derivative along search dir.
-    int awc(double, double, double, double);
+  double maxepsrot;
-    void make_step(double, double *);
+  double *spvec;		// variables for atomic dof, as 1d vector
-    double max_torque();
+  double *fmvec;		// variables for atomic dof, as 1d vector
-    double der_e_cur;		// current derivative along search dir.
+  double *g_old;  		// gradient vector at previous step
-    double der_e_pr;    	// previous derivative along search dir.
+  double *g_cur;  		// current gradient vector
-    int use_line_search; 	// use line search or not.
+  double *p_s;  		// search direction vector
-    double maxepsrot;
 
-    double *rho;        // estimation of curvature
+  void advance_spins();
-    double **ds;	// change in rotation matrix between two iterations, da
+  void calc_gradient();
-    double **dy;        // change in gradients between two iterations, dg
+  void calc_search_direction();
-    double **sp_copy;	// copy of the spins
+  void vm3(const double *, const double *, double *);
-    int num_mem;        // number of stored steps
+  void rodrigues_rotation(const double *, double *);
-    bigint last_negative;
+  void make_step(double, double *);
+  int calc_and_make_step(double, double, int);
+  int awc(double, double, double, double);
+  double maximum_rotation(double *);
+
+  double *rho;		// estimation of curvature
+  double **ds;		// change in rotation matrix between two iterations, da
+  double **dy;		// change in gradients between two iterations, dg
+  double **sp_copy;	// copy of the spins
+  int num_mem;		// number of stored steps
+  bigint last_negative;
 };
 
 }

src/min.cpp

@@ -42,10 +42,12 @@
 #include "output.h"
 #include "thermo.h"
 #include "timer.h"
+#include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -54,6 +56,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
   dmax = 0.1;
   searchflag = 0;
   linestyle = 1;
+  normstyle = 0;
 
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
@@ -653,6 +656,14 @@ void Min::modify_params(int narg, char **arg)
       if (strcmp(arg[iarg+1],"backtrack") == 0) linestyle = 0;
       else if (strcmp(arg[iarg+1],"quadratic") == 0) linestyle = 1;
       else if (strcmp(arg[iarg+1],"forcezero") == 0) linestyle = 2;
+      else if (strcmp(arg[iarg+1],"spin_cubic") == 0) linestyle = 3;
+      else if (strcmp(arg[iarg+1],"spin_none") == 0) linestyle = 4;
+      else error->all(FLERR,"Illegal min_modify command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"norm") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
+      if (strcmp(arg[iarg+1],"euclidean") == 0) normstyle = 0;
+      else if (strcmp(arg[iarg+1],"max") == 0) normstyle = 1;
       else error->all(FLERR,"Illegal min_modify command");
       iarg += 2;
     } else {
@@ -816,6 +827,69 @@ double Min::fnorm_inf()
   return norm_inf;
 }
 
+/* ----------------------------------------------------------------------
+   compute and return  sum_i||mag. torque_i||_2 (in eV)
+------------------------------------------------------------------------- */
+
+double Min::total_torque()
+{
+  double fmsq,ftotsqone,ftotsqall;
+  int nlocal = atom->nlocal;
+  double hbar = force->hplanck/MY_2PI;
+  double tx,ty,tz;
+  double **sp = atom->sp;
+  double **fm = atom->fm;
+
+  fmsq = ftotsqone = ftotsqall = 0.0;
+  for (int i = 0; i < nlocal; i++) {
+    tx = fm[i][1] * sp[i][2] - fm[i][2] * sp[i][1];
+    ty = fm[i][2] * sp[i][0] - fm[i][0] * sp[i][2];
+    tz = fm[i][0] * sp[i][1] - fm[i][1] * sp[i][0];
+    fmsq = tx * tx + ty * ty + tz * tz;
+    ftotsqone += fmsq;
+  }
+
+  // summing all fmsqtot on this replica
+
+  MPI_Allreduce(&ftotsqone,&ftotsqall,1,MPI_DOUBLE,MPI_SUM,world);
+
+  // multiply it by hbar so that units are in eV
+  
+  return sqrt(ftotsqall) * hbar;
+}
+
+/* ----------------------------------------------------------------------
+   compute and return max_i ||mag. torque_i|| (in eV)
+------------------------------------------------------------------------- */
+
+double Min::max_torque()
+{
+  double fmsq,fmaxsqone,fmaxsqall;
+  int nlocal = atom->nlocal;
+  double hbar = force->hplanck/MY_2PI;
+  double tx,ty,tz;
+  double **sp = atom->sp;
+  double **fm = atom->fm;
+
+  fmsq = fmaxsqone = fmaxsqall = 0.0;
+  for (int i = 0; i < nlocal; i++) {
+    tx = fm[i][1] * sp[i][2] - fm[i][2] * sp[i][1];
+    ty = fm[i][2] * sp[i][0] - fm[i][0] * sp[i][2];
+    tz = fm[i][0] * sp[i][1] - fm[i][1] * sp[i][0];
+    fmsq = tx * tx + ty * ty + tz * tz;
+    fmaxsqone = MAX(fmaxsqone,fmsq);
+  }
+
+  // finding max fm on this replica
+
+  fmaxsqall = fmaxsqone;
+  MPI_Allreduce(&fmaxsqone,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,world);
+
+  // multiply it by hbar so that units are in eV
+
+  return sqrt(fmaxsqall) * hbar;
+}
+
 /* ----------------------------------------------------------------------
    possible stop conditions
 ------------------------------------------------------------------------- */

src/min.h

@@ -42,6 +42,10 @@ class Min : protected Pointers {
   double fnorm_sqr();
   double fnorm_inf();
 
+  // methods for spin minimizers
+  double max_torque(); 
+  double total_torque();
+
   virtual void init_style() {}
   virtual void setup_style() = 0;
   virtual void reset_vectors() = 0;
@@ -56,8 +60,11 @@ class Min : protected Pointers {
   int virial_style;           // compute virial explicitly or implicitly
   int external_force_clear;   // clear forces locally or externally
 
-  double dmax;                // max dist to move any atom in one step
+  double dmax;			// max dist to move any atom in one step
-  int linestyle;              // 0 = backtrack, 1 = quadratic, 2 = forcezero
+  int linestyle;		// 0 = backtrack, 1 = quadratic, 2 = forcezero 
+  				// 3 = spin_cubic, 4 = spin_none
+
+  int normstyle;		// 0 = Euclidean norm, 1 = inf. norm
 
   int nelist_global,nelist_atom;    // # of PE,virial computes to check
   int nvlist_global,nvlist_atom;
@@ -102,9 +109,6 @@ class Min : protected Pointers {
   double energy_force(int);
   void force_clear();
 
-  double compute_force_norm_sqr();
-  double compute_force_norm_inf();
-
   void ev_setup();
   void ev_set(bigint);
 

src/min_cg.cpp

@@ -37,7 +37,7 @@ MinCG::MinCG(LAMMPS *lmp) : MinLineSearch(lmp) {}
 int MinCG::iterate(int maxiter)
 {
   int i,m,n,fail,ntimestep;
-  double beta,gg,dot[2],dotall[2];
+  double beta,gg,dot[2],dotall[2],fmax,fmaxall;
   double *fatom,*gatom,*hatom;
 
   // nlimit = max # of CG iterations before restarting
@@ -87,10 +87,12 @@ int MinCG::iterate(int maxiter)
 
     // force tolerance criterion
 
+    fmax = fmaxall = 0.0;
     dot[0] = dot[1] = 0.0;
     for (i = 0; i < nvec; i++) {
       dot[0] += fvec[i]*fvec[i];
       dot[1] += fvec[i]*g[i];
+      fmax = MAX(fmax,fvec[i]*fvec[i]);
     }
     if (nextra_atom)
       for (m = 0; m < nextra_atom; m++) {
@@ -100,16 +102,22 @@ int MinCG::iterate(int maxiter)
         for (i = 0; i < n; i++) {
           dot[0] += fatom[i]*fatom[i];
           dot[1] += fatom[i]*gatom[i];
+	  fmax = MAX(fmax,fatom[i]*fatom[i]);
         }
       }
     MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(&fmax,&fmaxall,2,MPI_DOUBLE,MPI_MAX,world);
     if (nextra_global)
       for (i = 0; i < nextra_global; i++) {
         dotall[0] += fextra[i]*fextra[i];
         dotall[1] += fextra[i]*gextra[i];
       }
 
-    if (dotall[0] < update->ftol*update->ftol) return FTOL;
+    if (normstyle == 1) {	// max force norm
+      if (fmax < update->ftol*update->ftol) return FTOL;
+    } else {			// Euclidean force norm
+      if (dotall[0] < update->ftol*update->ftol) return FTOL;
+    }
 
     // update new search direction h from new f = -Grad(x) and old g
     // this is Polak-Ribieri formulation

src/min_fire.cpp

@@ -80,7 +80,7 @@ void MinFire::reset_vectors()
 int MinFire::iterate(int maxiter)
 {
   bigint ntimestep;
-  double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfall;
+  double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfloc,fdotfall;
   double scale1,scale2;
   double dtvone,dtv,dtf,dtfm;
   int flag,flagall;
@@ -250,7 +250,15 @@ int MinFire::iterate(int maxiter)
     // sync across replicas if running multi-replica minimization
 
     if (update->ftol > 0.0) {
-      fdotf = fnorm_sqr();
+      if (normstyle == 1) {		// max force norm
+	fdotf = fnorm_inf();
+	fdotfloc = fdotf;
+	MPI_Allreduce(&fdotfloc,&fdotf,1,MPI_INT,MPI_MAX,universe->uworld);
+      } else {			// Euclidean force norm
+	fdotf = fnorm_sqr();
+	fdotfloc = fdotf;
+	MPI_Allreduce(&fdotfloc,&fdotf,1,MPI_INT,MPI_SUM,universe->uworld);
+      }
       if (update->multireplica == 0) {
         if (fdotf < update->ftol*update->ftol) return FTOL;
       } else {

src/min_hftn.cpp

@@ -20,6 +20,7 @@
 #include <cmath>
 #include <cstring>
 #include "atom.h"
+#include "error.h"
 #include "fix_minimize.h"
 #include "min_hftn.h"
 #include "modify.h"
@@ -111,6 +112,9 @@ void MinHFTN::init()
 {
   Min::init();
 
+  if (normstyle == 1)
+    error->all(FLERR,"Incorrect min_modify option");
+
   for (int  i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) {
     if (_daExtraGlobal[i] != NULL)
       delete [] _daExtraGlobal[i];

src/min_quickmin.cpp

@@ -76,7 +76,7 @@ void MinQuickMin::reset_vectors()
 int MinQuickMin::iterate(int maxiter)
 {
   bigint ntimestep;
-  double vmax,vdotf,vdotfall,fdotf,fdotfall,scale;
+  double vmax,vdotf,vdotfall,fdotf,fdotfloc,fdotfall,scale;
   double dtvone,dtv,dtf,dtfm;
   int flag,flagall;
 
@@ -216,7 +216,15 @@ int MinQuickMin::iterate(int maxiter)
     // sync across replicas if running multi-replica minimization
 
     if (update->ftol > 0.0) {
-      fdotf = fnorm_sqr();
+      if (normstyle == 1) {		// max force norm
+	fdotf = fnorm_inf();
+	fdotfloc = fdotf;
+	MPI_Allreduce(&fdotfloc,&fdotf,1,MPI_INT,MPI_MAX,universe->uworld);
+      } else {			// Euclidean force norm
+	fdotf = fnorm_sqr();
+	fdotfloc = fdotf;
+	MPI_Allreduce(&fdotfloc,&fdotf,1,MPI_INT,MPI_SUM,universe->uworld);
+      }
       if (update->multireplica == 0) {
         if (fdotf < update->ftol*update->ftol) return FTOL;
       } else {

src/min_sd.cpp

@@ -79,7 +79,8 @@ int MinSD::iterate(int maxiter)
 
     // force tolerance criterion
 
-    fdotf = fnorm_sqr();
+    if (normstyle == 1) fdotf = fnorm_inf();	// max force norm
+    else fdotf = fnorm_sqr();			// Euclidean force norm
     if (fdotf < update->ftol*update->ftol) return FTOL;
 
     // set new search direction h to f = -Grad(x)