Commit JT 092319

- modified norm input in min_modify
- corrected doc/src/min_modify.txt
- added expression of the norms
- added a min max method in src/min.h

doc/src/Eqs/norm_max.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}
+    % \left| \left| \vec{F} \right| \right|_2 
+    || \vec{F} ||_{max} 
+    = {\rm max}\left(||\vec{F}_1||, \cdots, ||\vec{F}_N||\right)
+    \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}

doc/src/Eqs/norm_two.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}
+    % \left| \left| \vec{F} \right| \right|_2 
+    || \vec{F} ||_{2} 
+    = \sqrt{\vec{F}_1+ \cdots + \vec{F}_N}
+    \nonumber
+  \end{equation}
+\end{varwidth}
+\end{document}

doc/src/min_modify.txt

@@ -18,8 +18,9 @@ keyword = {dmax} or {line} or {norm} or {alpha_damp} or {discrete_factor}
     max = maximum distance for line search to move (distance units)
   {line} value = {backtrack} or {quadratic} or {forcezero} or {spin_cubic} or {spin_none}
     backtrack,quadratic,forcezero,spin_cubic,spin_none = style of linesearch to use 
-  {norm} value = {euclidean} or {max}
+  {norm} value = {two} or {max}
-    euclidean,max = style of norm to use 
+    two = Euclidean two-norm (length of 3N vector)
+    max = max value of across all 3-vectors
   {alpha_damp} value = damping
     damping = fictitious Gilbert damping for spin minimization (adim)
   {discrete_factor} value = factor
@@ -74,9 +75,16 @@ could move in the gradient direction to reduce forces further.
 The choice of a norm can be modified for the min styles {cg}, {sd}, 
 {quickmin}, {fire}, {spin}, {spin/cg} and {spin/lbfgs} using 
 the {norm} keyword.
-The default {euclidean} norm computes the 2-norm (length) of the
+The default {two} norm computes the 2-norm (Euclidean length) of the
-global force vector. The {max} norm computes the maximum value
+global force vector:
-of the 2-norms across all forces in the system.  
+
+:c,image(Eqs/norm_two.jpg) 
+
+The {max} norm computes the length of the 3-vector force 
+for each atom  (2-norm), and takes the maximum value of those accross 
+all atoms
+
+:c,image(Eqs/norm_max.jpg) 
 
 Keywords {alpha_damp} and {discrete_factor} only make sense when
 a "min_spin"_min_spin.html command is declared. 

src/MAKE/Makefile.serial

@@ -7,7 +7,7 @@ SHELL = /bin/sh
 # specify flags and libraries needed for your compiler
 
 CC =		g++
-CCFLAGS =	-g -O3 -Wall
+CCFLAGS =	-g -O3
 SHFLAGS =	-fPIC
 DEPFLAGS =	-M
 

src/min.cpp

@@ -49,6 +49,8 @@
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
+enum{TWO,MAX}
+
 /* ---------------------------------------------------------------------- */
 
 Min::Min(LAMMPS *lmp) : Pointers(lmp)
@@ -56,7 +58,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
   dmax = 0.1;
   searchflag = 0;
   linestyle = 1;
-  normstyle = 0;
+  normstyle = TWO;
 
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
@@ -662,8 +664,8 @@ void Min::modify_params(int narg, char **arg)
       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;
+      if (strcmp(arg[iarg+1],"two") == 0) normstyle = TWO;
-      else if (strcmp(arg[iarg+1],"max") == 0) normstyle = 1;
+      else if (strcmp(arg[iarg+1],"max") == 0) normstyle = MAX;
       else error->all(FLERR,"Illegal min_modify command");
       iarg += 2;
     } else {
@@ -827,6 +829,41 @@ double Min::fnorm_inf()
   return norm_inf;
 }
 
+/* ----------------------------------------------------------------------
+   compute and return ||force||_max (inf norm per-vector)
+------------------------------------------------------------------------- */
+
+double Min::fnorm_max()
+{
+  int i,n;
+  double fdotf,*fatom;
+
+  double local_norm_max = 0.0;
+  for (i = 0; i < nvec; i+=3) {
+    fdotf = fvec[i]*fvec[i]+fvec[i+1]*fvec[i+1]+fvec[i+2]*fvec[i+2];
+    local_norm_max = MAX(fdotf,local_norm_max);
+  }
+  if (nextra_atom) {
+    for (int m = 0; m < nextra_atom; m++) {
+      fatom = fextra_atom[m];
+      n = extra_nlen[m];
+      for (i = 0; i < n; i+=3)
+        fdotf = fvec[i]*fvec[i]+fvec[i+1]*fvec[i+1]+fvec[i+2]*fvec[i+2];
+        local_norm_max = MAX(fdotf,local_norm_max);
+    }
+  }
+
+  double norm_max = 0.0;
+  MPI_Allreduce(&local_norm_max,&norm_max,1,MPI_DOUBLE,MPI_MAX,world);
+
+  if (nextra_global)
+    for (i = 0; i < n; i+=3)
+      fdotf = fvec[i]*fvec[i]+fvec[i+1]*fvec[i+1]+fvec[i+2]*fvec[i+2];
+      norm_max = MAX(fdotf,norm_max);
+
+  return norm_max;
+}
+
 /* ----------------------------------------------------------------------
    compute and return  sum_i||mag. torque_i||_2 (in eV)
 ------------------------------------------------------------------------- */
@@ -842,10 +879,10 @@ double Min::total_torque()
 
   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];
+    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];
+    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];
+    tz = fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0];
-    fmsq = tx * tx + ty * ty + tz * tz;
+    fmsq = tx*tx + ty*ty + tz*tz;
     ftotsqone += fmsq;
   }
 
@@ -873,10 +910,10 @@ double Min::max_torque()
 
   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];
+    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];
+    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];
+    tz = fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0];
-    fmsq = tx * tx + ty * ty + tz * tz;
+    fmsq = tx*tx + ty*ty + tz*tz;
     fmaxsqone = MAX(fmaxsqone,fmsq);
   }
 

src/min.h

@@ -41,6 +41,7 @@ class Min : protected Pointers {
   virtual int modify_param(int, char **) {return 0;}
   double fnorm_sqr();
   double fnorm_inf();
+  double fnorm_max();
 
   // methods for spin minimizers
   double max_torque(); 
@@ -64,7 +65,8 @@ class Min : protected Pointers {
   int linestyle;		// 0 = backtrack, 1 = quadratic, 2 = forcezero 
   				// 3 = spin_cubic, 4 = spin_none
 
-  int normstyle;		// 0 = Euclidean norm, 1 = inf. norm
+  int normstyle;		// TWO or MAX flag for force norm evaluation
+  // int normstyle;		// 0 = Euclidean norm, 1 = inf. norm
 
   int nelist_global,nelist_atom;    // # of PE,virial computes to check
   int nvlist_global,nvlist_atom;

src/min_cg.cpp

@@ -35,7 +35,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],fmax,fmaxall;
+  double beta,gg,dot[2],dotall[2],fmax;
   double *fatom,*gatom,*hatom;
 
   // nlimit = max # of CG iterations before restarting
@@ -85,13 +85,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++) {
         fatom = fextra_atom[m];
@@ -104,16 +103,17 @@ int MinCG::iterate(int maxiter)
         }
       }
     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 (normstyle == 1) {	// max force norm
+    fmax = 0.0;
+    if (normstyle == MAX) {	// max force norm
+      fmax = fnorm_max();
       if (fmax < update->ftol*update->ftol) return FTOL;
-    } else {			// Euclidean force norm
+    } else {			// Euclidean force 2-norm
       if (dotall[0] < update->ftol*update->ftol) return FTOL;
     }
 

src/min_fire.cpp

@@ -250,15 +250,8 @@ int MinFire::iterate(int maxiter)
     // sync across replicas if running multi-replica minimization
 
     if (update->ftol > 0.0) {
-      if (normstyle == 1) {		// max force norm
+      if (normstyle == MAX) fdotf = fnorm_inf();	// max force norm
-	fdotf = fnorm_inf();
+      else fdotf = fnorm_sqr();				// Euclidean force 2-norm
-	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

@@ -113,7 +113,7 @@ void MinHFTN::init()
 {
   Min::init();
 
-  if (normstyle == 1)
+  if (normstyle == MAX)
     error->all(FLERR,"Incorrect min_modify option");
 
   for (int  i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) {

src/min_quickmin.cpp

@@ -215,15 +215,8 @@ int MinQuickMin::iterate(int maxiter)
     // sync across replicas if running multi-replica minimization
 
     if (update->ftol > 0.0) {
-      if (normstyle == 1) {		// max force norm
+      if (normstyle == MAX) fdotfloc = fnorm_max();	// max force norm
-	fdotf = fnorm_inf();
+      else fdotf = fnorm_sqr();				// Euclidean force 2-norm
-	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

@@ -34,7 +34,7 @@ MinSD::MinSD(LAMMPS *lmp) : MinLineSearch(lmp) {}
 int MinSD::iterate(int maxiter)
 {
   int i,m,n,fail,ntimestep;
-  double fdotf;
+  double fdotf,fdotfloc;
   double *fatom,*hatom;
 
   // initialize working vectors
@@ -77,8 +77,8 @@ int MinSD::iterate(int maxiter)
 
     // force tolerance criterion
 
-    if (normstyle == 1) fdotf = fnorm_inf();	// max force norm
+    if (normstyle == MAX) fdotf = fnorm_max();	// max force norm
-    else fdotf = fnorm_sqr();			// Euclidean force norm
+    else fdotf = fnorm_sqr();			// Euclidean force 2-norm
     if (fdotf < update->ftol*update->ftol) return FTOL;
 
     // set new search direction h to f = -Grad(x)