New options adaptstep and halfstepback in adaptglok. Adapted documentation.

doc/src/min_modify.txt

@@ -19,14 +19,17 @@ keyword = {dmax} or {line} or {delaystep} or {dt_grow} or {dt_shrink}
     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
   {alpha0} value = {float}  for adaptglok
   {alpha_shrink} value = {float}  for adaptglok
   {tmax} value = {float}  for adaptglok
-  {tmin} value = {float}  for adaptglok :pre
+  {tmin} value = {float}  for adaptglok
+  {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} 
+                       or {leapfrog}  for adaptglok
+  {adaptstep} arg = {yes} or {no}  for adaptglok
+  {halfstepback} value = {yes} or {no}  for adaptglok :pre
 :ule
 
 [Examples:]
@@ -79,6 +82,9 @@ 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.
+In particular for debugging purpose, it is possible to switch off the 
+variable timestep and the backstep when downhill, by using {adaptstep} 
+and {halfstepback}.
 
 [Restrictions:] none
 
@@ -90,4 +96,4 @@ to adapt them for particular cases.
 
 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.
+tmax = 2.0, tmin = 0.02, adaptstep = yes and halfstepback = yes.

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 different integration
+timestep, a backstep when downhill and different integration
-schemes (Euler, Leap Frog, Velocity-Verlet).
+schemes.
 
 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

@@ -63,6 +63,8 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
   tmax = 2.0;
   tmin = 0.02;
   integrator = 0;
+  halfstepback_flag = 1;
+  adaptstep_flag = 1;
 
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
@@ -680,6 +682,18 @@ void Min::modify_params(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
       tmin = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;       
+     } else if (strcmp(arg[iarg],"halfstepback") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
+      if (strcmp(arg[iarg+1],"yes") == 0) halfstepback_flag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) halfstepback_flag = 0;
+      else error->all(FLERR,"Illegal min_modify command");
+      iarg += 2;       
+     } else if (strcmp(arg[iarg],"adaptstep") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
+      if (strcmp(arg[iarg+1],"yes") == 0) adaptstep_flag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) adaptstep_flag = 0;
+      else error->all(FLERR,"Illegal min_modify command");
+      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],"eulerimplicit") == 0) integrator = 0;

src/min.h

@@ -63,6 +63,8 @@ class Min : protected Pointers {
   double alpha0,alpha_shrink; // mixing velocities+forces coefficient (adaptglok)
   double tmax,tmin;           // timestep max, min (adaptglok)
   int integrator;             // choose the style of time integrator (adaptglok)
+  int halfstepback_flag;      // 1: half step backward when v.f <= 0.0 (adaptglok)
+  int adaptstep_flag;         // 1: variable timestep, 0: constant (adaptglok)
 
   int nelist_global,nelist_atom;    // # of PE,virial computes to check
   int nvlist_global,nvlist_atom;

src/min_adaptglok.cpp

@@ -184,7 +184,7 @@ int MinAdaptGlok::iterate(int maxiter)
       if (fdotfall <= 1e-20) scale2 = 0.0;
       else scale2 = alpha * sqrt(vdotvall/fdotfall);
 
-      if (ntimestep - last_negative > delaystep) {
+      if (ntimestep - last_negative > delaystep && adaptstep_flag) {
         dt = MIN(dt*dt_grow,dtmax);
         update->dt = dt;
         alpha *= alpha_shrink;
@@ -201,7 +201,7 @@ int MinAdaptGlok::iterate(int maxiter)
     } else {
       last_negative = ntimestep;
       // Limit decrease of timestep
-      if (ntimestep - ntimestep_fire > delaystep) {
+      if (ntimestep - ntimestep_fire > delaystep && adaptstep_flag) {
         alpha = alpha0;
         if (dt > dtmin) {
           dt *= dt_shrink;
@@ -209,11 +209,13 @@ int MinAdaptGlok::iterate(int maxiter)
         }
       }
       double **x = atom->x;
+      if (halfstepback_flag) {
         for (int i = 0; i < nlocal; i++) {
           x[i][0] -= 0.5 * dtv * v[i][0];
           x[i][1] -= 0.5 * dtv * v[i][1];
           x[i][2] -= 0.5 * dtv * v[i][2];
         }
+      }
       for (int i = 0; i < nlocal; i++)
         v[i][0] = v[i][1] = v[i][2] = 0.0;
     }