Correct calculation of the Forces within the Verlet integrator. Updated doc.

doc/src/min_modify.txt

@@ -14,11 +14,12 @@ min_modify keyword values ... :pre
 
 one or more keyword/value pairs may be listed :ulb,l
 keyword = {dmax} or {line} or {delaystep} or {dt_grow} or {dt_shrink}
-          or {alpha0} or {alpha_shrink} or {tmax} or {tmin}
+          or {alpha0} or {alpha_shrink} or {tmax} or {tmin} or {integrator}
   {dmax} value = max
     max = maximum distance for line search to move (distance units)
   {line} value = {backtrack} or {quadratic} or {forcezero}
     backtrack,quadratic,forcezero = style of linesearch to use
+  {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} or {leapfrog}  for adaptglok :pre
   {delaystep} value = {int} for adaptglok
   {dt_grow} value = {float}  for adaptglok
   {dt_shrink} value = {float}  for adaptglok
@@ -74,10 +75,10 @@ that difference may be smaller than machine epsilon even if atoms
 could move in the gradient direction to reduce forces further.
 
 The style {adaptglok} has several parameters that can be tuned in order
-to optimize the relaxation: {delaystep}, {dt_grow}, {dt_shrink}, 
+to optimize the relaxation: {integrator}, {delaystep}, {dt_grow}, {dt_shrink}, 
 {alpha0}, {alpha_shrink}, {tmax} and {tmin}. The default values have 
 been chosen to be reliable in most cases, but one could be willing 
-to adapt them for particular cases. 
+to adapt them for particular cases.
 
 [Restrictions:] none
 
@@ -87,6 +88,6 @@ to adapt them for particular cases.
 
 [Default:]
 
-The option defaults are dmax = 0.1,  line = quadratic, delaystep = 20, 
+The option defaults are dmax = 0.1,  line = quadratic, integrator = eulerimplicit, delaystep = 20, 
 dt_grow = 1.1, dt_shrink = 0.5, alpha0 = 0.25, alpha_shrink = 0.99, 
 tmax = 2.0 and tmin = 0.02.

doc/src/min_style.txt

@@ -64,8 +64,8 @@ a minimization.
 Style {adaptglok} is a re-implementation of the style {fire} to better 
 fit the published and unpublished optimizations of the method described
 in "(Bitzek)"_#Bitzek. It includes a temporization of the variable
-timestep, a backstep when uphill and a Velocity-Verlet integration
+timestep, a backstep when uphill and different integration
-scheme.
+schemes (Euler, Leap Frog, Velocity-Verlet).
 
 Either the {quickmin} and {fire} styles are useful in the context of
 nudged elastic band (NEB) calculations via the "neb"_neb.html command.

src/min.cpp

@@ -682,10 +682,10 @@ void Min::modify_params(int narg, char **arg)
       iarg += 2;       
      } else if (strcmp(arg[iarg],"integrator") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      if (strcmp(arg[iarg+1],"euler") == 0) integrator = 0;
+      if (strcmp(arg[iarg+1],"eulerimplicit") == 0) integrator = 0;
       else if (strcmp(arg[iarg+1],"verlet") == 0) integrator = 1;
       else if (strcmp(arg[iarg+1],"leapfrog") == 0) integrator = 2;
-      else if (strcmp(arg[iarg+1],"eulerstandard") == 0) integrator = 3;
+      else if (strcmp(arg[iarg+1],"eulerexplicit") == 0) integrator = 3;
      else error->all(FLERR,"Illegal min_modify command");
       iarg += 2;       
    } else if (strcmp(arg[iarg],"line") == 0) {

src/min_adaptglok.cpp

@@ -92,6 +92,10 @@ int MinAdaptGlok::iterate(int maxiter)
     double *mass = atom->mass;
     int *type = atom->type;
     int nlocal = atom->nlocal;
+
+    energy_force(0);
+    neval++;
+
     dtf = 0.5 * dt * force->ftm2v;
     if (rmass) {
       for (int i = 0; i < nlocal; i++) {
@@ -277,6 +281,10 @@ int MinAdaptGlok::iterate(int maxiter)
       }
     }
 
+    eprevious = ecurrent;
+    ecurrent = energy_force(0);
+    neval++;
+
 
   // Velocity Verlet integration
 
@@ -316,8 +324,9 @@ int MinAdaptGlok::iterate(int maxiter)
       }
     }
     
-    double **f = atom->f;
+    eprevious = ecurrent;
-    double **v = atom->v;
+    ecurrent = energy_force(0);
+    neval++;
 
     if (rmass) {
       for (int i = 0; i < nlocal; i++) {
@@ -373,6 +382,10 @@ int MinAdaptGlok::iterate(int maxiter)
       }
     }
 
+    eprevious = ecurrent;
+    ecurrent = energy_force(0);
+    neval++;
+
   // Standard Euler integration
 
   }else if (integrator == 3) {
@@ -411,12 +424,12 @@ int MinAdaptGlok::iterate(int maxiter)
       }
     }
 
-  }
-
     eprevious = ecurrent;
     ecurrent = energy_force(0);
     neval++;
 
+  }
+
     // energy tolerance criterion
     // only check after delaystep elapsed since velocties reset to 0
     // sync across replicas if running multi-replica minimization