Merge branch 'master' into imd-improvements

doc/fix_efield.html

@@ -31,11 +31,11 @@ fix kick external-field efield 0.0 0.0 v_oscillate
 external electric field being applied to the system.
 </P>
 <P>Any of the 3 quantities defining the E-field components can be
-specified as an equal-style or atom-style <A HREF = "variable">variable</A>, namely
-<I>ex</I>, <I>ey</I>, <I>ez</I>.  If the value is a variable, it should be specified
-as v_ID, where ID is the variable ID.  In this case, the variable will
-be evaluated each timestep, and its value used to determine the
-E-field component.
+specified as an equal-style or atom-style <A HREF = "variable.html">variable</A>,
+namely <I>ex</I>, <I>ey</I>, <I>ez</I>.  If the value is a variable, it should be
+specified as v_ID, where ID is the variable ID.  In this case, the
+variable will be evaluated each timestep, and its value used to
+determine the E-field component.
 </P>
 <P>Equal-style variables can specify formulas with various mathematical
 functions, and include <A HREF = "thermo_style.html">thermo_style</A> command

doc/fix_efield.txt

@@ -28,11 +28,11 @@ Add a force F = qE to each charged atom in the group due to an
 external electric field being applied to the system.
 
 Any of the 3 quantities defining the E-field components can be
-specified as an equal-style or atom-style "variable"_variable, namely
-{ex}, {ey}, {ez}.  If the value is a variable, it should be specified
-as v_ID, where ID is the variable ID.  In this case, the variable will
-be evaluated each timestep, and its value used to determine the
-E-field component.
+specified as an equal-style or atom-style "variable"_variable.html,
+namely {ex}, {ey}, {ez}.  If the value is a variable, it should be
+specified as v_ID, where ID is the variable ID.  In this case, the
+variable will be evaluated each timestep, and its value used to
+determine the E-field component.
 
 Equal-style variables can specify formulas with various mathematical
 functions, and include "thermo_style"_thermo_style.html command

doc/minimize.html

@@ -151,13 +151,22 @@ reduced.
 </P>
 <HR>
 
-<P>IMPORTANT NOTE: It is highly recommended that you use a <A HREF = "pair_style.html">pair
-style</A> that goes to 0.0 at the cutoff distance when
-performing minimization (even if you later change it when running
-dynamics).  If this is not done, the total energy of the system will
-have discontinuities when the relative distance between any pair of
-atoms changes from cutoff+epsilon to cutoff-epsilon and the minimizer
-may behave poorly.
+<P>IMPORTANT NOTE: There are several force fields in LAMMPS which have
+discontinuities or other approximations which may prevent you from
+performing an energy minimization to high tolerances.  For example,
+you should use a <A HREF = "pair_style.html">pair style</A> that goes to 0.0 at the
+cutoff distance when performing minimization (even if you later change
+it when running dynamics).  If you do not do this, the total energy of
+the system will have discontinuities when the relative distance
+between any pair of atoms changes from cutoff+epsilon to
+cutoff-epsilon and the minimizer may behave poorly.  Some of the
+manybody potentials use splines and other internal cutoffs that
+inherently have this problem.  The <A HREF = "kspace_style.html">long-range Coulombic
+styles</A> (PPPM, Ewald) are approximate to within the
+user-specified tolerance, which means their energy and forces may not
+agree to a higher precision than the Kspace-specified tolerance.  In
+all these cases, the minimizer may give up and stop before finding a
+minimum to the specified energy or force tolerance.
 </P>
 <P>Note that a cutoff Lennard-Jones potential (and others) can be shifted
 so that its energy is 0.0 at the cutoff via the

doc/minimize.txt

@@ -148,13 +148,22 @@ The iterations and force evaluation values are what is checked by the
 
 :line
 
-IMPORTANT NOTE: It is highly recommended that you use a "pair
-style"_pair_style.html that goes to 0.0 at the cutoff distance when
-performing minimization (even if you later change it when running
-dynamics).  If this is not done, the total energy of the system will
-have discontinuities when the relative distance between any pair of
-atoms changes from cutoff+epsilon to cutoff-epsilon and the minimizer
-may behave poorly.
+IMPORTANT NOTE: There are several force fields in LAMMPS which have
+discontinuities or other approximations which may prevent you from
+performing an energy minimization to high tolerances.  For example,
+you should use a "pair style"_pair_style.html that goes to 0.0 at the
+cutoff distance when performing minimization (even if you later change
+it when running dynamics).  If you do not do this, the total energy of
+the system will have discontinuities when the relative distance
+between any pair of atoms changes from cutoff+epsilon to
+cutoff-epsilon and the minimizer may behave poorly.  Some of the
+manybody potentials use splines and other internal cutoffs that
+inherently have this problem.  The "long-range Coulombic
+styles"_kspace_style.html (PPPM, Ewald) are approximate to within the
+user-specified tolerance, which means their energy and forces may not
+agree to a higher precision than the Kspace-specified tolerance.  In
+all these cases, the minimizer may give up and stop before finding a
+minimum to the specified energy or force tolerance.
 
 Note that a cutoff Lennard-Jones potential (and others) can be shifted
 so that its energy is 0.0 at the cutoff via the

src/fix_box_relax.cpp

@@ -24,6 +24,7 @@
 #include "update.h"
 #include "comm.h"
 #include "force.h"
+#include "kspace.h"
 #include "modify.h"
 #include "compute.h"
 #include "error.h"
@@ -318,6 +319,9 @@ void FixBoxRelax::init()
 
   pv2e = 1.0 / force->nktv2p;
 
+  if (force->kspace) kspace_flag = 1;
+  else kspace_flag = 0;
+
   // detect if any rigid fixes exist so rigid bodies move when box is remapped
   // rfix[] = indices to each fix rigid
 
@@ -528,6 +532,7 @@ void FixBoxRelax::min_step(double alpha, double *hextra)
     }
   }
   remap();
+  if (kspace_flag) force->kspace->setup();
 }
 
 /* ----------------------------------------------------------------------

src/fix_box_relax.h

@@ -50,6 +50,7 @@ class FixBoxRelax : public Fix {
   double p_target[6],p_current[6];
   double vol0,xprdinit,yprdinit,zprdinit;
   double vmax,pv2e,pflagsum;
+  int kspace_flag;
 
   int current_lifo;              // LIFO stack pointer
   double boxlo0[2][3];           // box bounds at start of line search

src/min_linesearch.cpp

@@ -409,7 +409,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha,
   // double alphatmp = alphamax*1.0e-4;
   // etmp = alpha_step(alphatmp,1,nfunc);
   // printf("alpha = %g dele = %g dele_force = %g err = %g\n",
-  //        alphatmp,etmp-eoriginal,-alphatmp*fdothall,
+  //   	    alphatmp,etmp-eoriginal,-alphatmp*fdothall,
   //        etmp-eoriginal+alphatmp*fdothall);
   // alpha_step(0.0,1,nfunc);
 
@@ -462,13 +462,15 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha,
     relerr = fabs(1.0-(0.5*(alpha-alphaprev)*(fh+fhprev)+ecurrent)/engprev);
     alpha0 = alpha - (alpha-alphaprev)*fh/delfh;
 
-    if (relerr <= QUADRATIC_TOL && alpha0 > 0.0) {
+    if (relerr <= QUADRATIC_TOL && alpha0 > 0.0 && alpha0 < alphamax) {
       ecurrent = alpha_step(alpha0,1,nfunc);
-      if (nextra_global) {
-	int itmp = modify->min_reset_ref();
-	if (itmp) ecurrent = energy_force(1);
+      if (ecurrent < eoriginal) {
+	if (nextra_global) {
+	  int itmp = modify->min_reset_ref();
+	  if (itmp) ecurrent = energy_force(1);
+	}
+	return 0;
       }
-      return 0;
     }
 
     // if backtracking energy change is better than ideal, exit with success

src/version.h

@@ -1 +1 @@
-#define LAMMPS_VERSION "26 Apr 2010"
+#define LAMMPS_VERSION "27 Apr 2010"