diff --git a/doc/src/Eqs/norm_max.jpg b/doc/src/Eqs/norm_max.jpg new file mode 100644 index 0000000000..c10db9a531 Binary files /dev/null and b/doc/src/Eqs/norm_max.jpg differ diff --git a/doc/src/Eqs/norm_max.tex b/doc/src/Eqs/norm_max.tex new file mode 100644 index 0000000000..3b2198bdf0 --- /dev/null +++ b/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} diff --git a/doc/src/Eqs/norm_two.jpg b/doc/src/Eqs/norm_two.jpg new file mode 100644 index 0000000000..5554de32f9 Binary files /dev/null and b/doc/src/Eqs/norm_two.jpg differ diff --git a/doc/src/Eqs/norm_two.tex b/doc/src/Eqs/norm_two.tex new file mode 100644 index 0000000000..d428081a49 --- /dev/null +++ b/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} diff --git a/doc/src/min_modify.txt b/doc/src/min_modify.txt index 06c1f7514f..b941a33559 100644 --- a/doc/src/min_modify.txt +++ b/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} - euclidean,max = style of norm to use + {norm} value = {two} or {max} + 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 -global force vector. The {max} norm computes the maximum value -of the 2-norms across all forces in the system. +The default {two} norm computes the 2-norm (Euclidean length) of the +global force vector: + +: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. diff --git a/src/MAKE/Makefile.serial b/src/MAKE/Makefile.serial index 8628d2bb73..5954d97761 100644 --- a/src/MAKE/Makefile.serial +++ b/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 diff --git a/src/min.cpp b/src/min.cpp index 31adf58525..33643d4837 100644 --- a/src/min.cpp +++ b/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; - else if (strcmp(arg[iarg+1],"max") == 0) normstyle = 1; + if (strcmp(arg[iarg+1],"two") == 0) normstyle = TWO; + 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]; - 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; + 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; } @@ -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]; - 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; + 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); } diff --git a/src/min.h b/src/min.h index e18d0dd677..ac7a3c1e9b 100644 --- a/src/min.h +++ b/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; diff --git a/src/min_cg.cpp b/src/min_cg.cpp index 2267a1ebb6..d98ec0ef97 100644 --- a/src/min_cg.cpp +++ b/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; } diff --git a/src/min_fire.cpp b/src/min_fire.cpp index 5b047ccd0e..dbb7f36148 100644 --- a/src/min_fire.cpp +++ b/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 - 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 (normstyle == MAX) fdotf = fnorm_inf(); // max force norm + else fdotf = fnorm_sqr(); // Euclidean force 2-norm if (update->multireplica == 0) { if (fdotf < update->ftol*update->ftol) return FTOL; } else { diff --git a/src/min_hftn.cpp b/src/min_hftn.cpp index 3c2cff9205..63432aab63 100644 --- a/src/min_hftn.cpp +++ b/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++) { diff --git a/src/min_quickmin.cpp b/src/min_quickmin.cpp index 3450f7785e..6846aaba0a 100644 --- a/src/min_quickmin.cpp +++ b/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 - 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 (normstyle == MAX) fdotfloc = fnorm_max(); // max force norm + else fdotf = fnorm_sqr(); // Euclidean force 2-norm if (update->multireplica == 0) { if (fdotf < update->ftol*update->ftol) return FTOL; } else { diff --git a/src/min_sd.cpp b/src/min_sd.cpp index dd59c9d2d6..3ded85990d 100644 --- a/src/min_sd.cpp +++ b/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 - else fdotf = fnorm_sqr(); // Euclidean force norm + if (normstyle == MAX) fdotf = fnorm_max(); // max 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)