Removing relaxbox code. Cleanup of the code. Add some parameter check. Improved documentation.

doc/src/min_modify.txt

@@ -19,19 +19,24 @@ 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
-  {delaystep} value = {int} for adaptglok
+  {delaystep} value = steps
-  {dt_grow} value = {float}  for adaptglok
+    steps = number of steps for dynamics temporization, adaptglok only
-  {dt_shrink} value = {float} for adaptglok
+  {dtgrow} value = grow
-  {alpha0} value = {float}  for adaptglok
+    grow = factor increasing the timestep, adaptglok only
-  {alpha_shrink} value = {float}  for adaptglok
+  {dtshrink} value = shrink
-  {tmax} value = {float}  for adaptglok
+    shrink = factor decreasing the timestep, adaptglok only
-  {tmin} value = {float}  for adaptglok
+  {alpha0} value = alpha
+    alpha = coefficient mixing velocities and forces, adaptglok only
+  {alphashrink} value = shrink
+    shrink = factor decreasing alpha, adaptglok only
+  {tmax} value = coef
+    coef = factor defining the maximum timestep, adaptglok only
+  {tmin} value = coef
+    coef = factor defining the minimum timestep, adaptglok only
   {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} 
-                       or {leapfrog}  for adaptglok
+                       or {leapfrog}
-  {halfstepback} value = {yes} or {no} for adaptglok :pre
+    eulerimplicit,eulerexplicit,verlet,leapfrog = integration scheme, adaptglok
-  {relaxbox} value = {no} or {iso} or {axial} for adaptglok 
+  {halfstepback} value = {yes} or {no}, adaptglok only :pre
-  {relaxbox_modulus} value = {float} for adaptglok 
-  {relaxbox_rate} value = {float} for adaptglok 
 :ule
 
 [Examples:]
@@ -80,18 +85,20 @@ 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: {integrator}, {delaystep}, {dt_grow}, {dt_shrink}, 
+to optimize the relaxation: {integrator}, {delaystep}, {dtgrow}, {dtshrink}, 
-{alpha0}, {alpha_shrink}, {tmax} and {tmin}. The default values have 
+{alpha0}, {alphashrink}, {tmax} and {tmin}.
-been chosen to be reliable in most cases, but one could be willing 
+Different Newton {integrator} can be selected: explicit Euler, 
-to adapt them for particular cases.
+semi-implicit Euler (= symplectic Euler), leapfrog, and velocity Verlet.
-In particular for debugging purpose, it is possible to switch off the 
+The parameters {tmax} and {tmin} define the maximum and minimum timestep 
-the backstep when downhill, by using and {halfstepback}.
+allowed during an adaptglok minimization. Those are not in time unit, but are 
-
+multiplication factor applied to the "timestep"_timestep.html. For example, 
-The style {adaptglok} allows for a basic box relaxation to a stress-free 
+{tmax} = 2.0 means that the maximum value the timestep can reach is twice the 
-state by choosing {relaxbox} {iso} or {axial}. One should give an
+value defined by "timestep"_timestep.html.
-approximated value for the bulk modulus {relaxbox_modulus} of the system.
+Note that even with the default parameters being chosen to be reliable in most
-The default value of the relaxation rate {relaxbox_rate} should be 
+cases, adjusting "timestep"_timestep.html, {tmax} and {tmin} should be consider
-optimum in most situations.
+to optimize the minimization, in particular for large/complex system.
+For debugging purposes, it is possible to switch off the back step when 
+downhill with {halfstepback} = 0.
 
 [Restrictions:] none
 
@@ -101,7 +108,7 @@ optimum in most situations.
 
 [Default:]
 
-The option defaults are dmax = 0.1,  line = quadratic, integrator = eulerimplicit, delaystep = 20, 
+The option defaults are dmax = 0.1, line = quadratic, 
-dt_grow = 1.1, dt_shrink = 0.5, alpha0 = 0.25, alpha_shrink = 0.99, 
+integrator = eulerimplicit, delaystep = 20, dt_grow = 1.1, dt_shrink = 0.5, 
-tmax = 2.0, tmin = 0.02,  adaptstep = yes, halfstepback = yes,
+alpha0 = 0.25, alpha_shrink = 0.99, tmax = 2.0, tmin = 0.02,  adaptstep = yes 
-relaxbox = no, relaxbox_modulus = 1000000 and relaxbox_rate = 0.33.
+and halfstepback = yes.

doc/src/min_style.txt

@@ -66,10 +66,8 @@ fit the published and unpublished optimizations of the method described
 in "(Bitzek)"_#Bitzek. It includes a temporization of the variable
 timestep, a backstep when downhill and different integration
 schemes.
-Note that a basic scheme for having a stress-free simulation box is 
-implemented within that style.
 
-Either the {quickmin} and {fire} styles are useful in the context of
+Either the {quickmin}, {fire} and {adaptglok} styles are useful in the context of
 nudged elastic band (NEB) calculations via the "neb"_neb.html command.
 
 NOTE: The damped dynamic minimizers use whatever timestep you have
@@ -77,7 +75,7 @@ defined via the "timestep"_timestep.html command.  Often they will
 converge more quickly if you use a timestep about 10x larger than you
 would normally use for dynamics simulations.
 
-NOTE: The {quickmin}, {fire}, and {hftn} styles do not yet support the
+NOTE: The {quickmin}, {fire}, {adaptglok}, and {hftn} styles do not yet support the
 use of the "fix box/relax"_fix_box_relax.html command or minimizations
 involving the electron radius in "eFF"_pair_eff.html models.
 

doc/src/minimize.txt

@@ -64,9 +64,9 @@ backtracking method is described in Nocedal and Wright's Numerical
 Optimization (Procedure 3.1 on p 41).
 
 The "minimization styles"_min_style.html {quickmin}, {fire} and 
-{adaptglok} perform damped dynamics using an Euler integration step 
+{adaptglok} perform damped dynamics using an Euler integration step. The style
-(Velocity-Verlet for {adaptglok}). Thus they require a 
+{adaptglok} can also use a leapfrog or velocity Verlet integration step.
-"timestep"_timestep.html be defined.
+Thus they require a "timestep"_timestep.html be defined.
 
 NOTE: The damped dynamic minimizers use whatever timestep you have
 defined via the "timestep"_timestep.html command.  Often they will

src/min.cpp

@@ -56,17 +56,14 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
   linestyle = 1;
 
   delaystep = 20;
-  dt_grow = 1.1;
+  dtgrow = 1.1;
-  dt_shrink = 0.5;
+  dtshrink = 0.5;
   alpha0 = 0.25;
-  alpha_shrink = 0.99;
+  alphashrink = 0.99;
   tmax = 2.0;
   tmin = 0.02;
   integrator = 0;
-  relaxbox_modulus = 1000000;
-  relaxbox_rate = 0.33;
   halfstepback_flag = 1;
-  relaxbox_flag = 0;
 
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
@@ -660,21 +657,21 @@ void Min::modify_params(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
       delaystep = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
-    } else if (strcmp(arg[iarg],"dt_grow") == 0) {
+    } else if (strcmp(arg[iarg],"dtgrow") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      dt_grow = force->numeric(FLERR,arg[iarg+1]);
+      dtgrow = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
-    } else if (strcmp(arg[iarg],"dt_shrink") == 0) {
+    } else if (strcmp(arg[iarg],"dtshrink") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      dt_shrink = force->numeric(FLERR,arg[iarg+1]);
+      dtshrink = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
     } else if (strcmp(arg[iarg],"alpha0") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
       alpha0 = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
-    } else if (strcmp(arg[iarg],"alpha_shrink") == 0) {
+    } else if (strcmp(arg[iarg],"alphashrink") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      alpha_shrink = force->numeric(FLERR,arg[iarg+1]);
+      alphashrink = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;
     } else if (strcmp(arg[iarg],"tmax") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
@@ -690,28 +687,13 @@ void Min::modify_params(int narg, char **arg)
       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],"relaxbox") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      if (strcmp(arg[iarg+1],"no") == 0) relaxbox_flag = 0;
-      else if (strcmp(arg[iarg+1],"iso") == 0) relaxbox_flag = 1;
-      else if (strcmp(arg[iarg+1],"axial") == 0) relaxbox_flag = 2;
-      else error->all(FLERR,"Illegal min_modify command");
-      iarg += 2;       
-    } else if (strcmp(arg[iarg],"relaxbox_modulus") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      relaxbox_modulus = force->numeric(FLERR,arg[iarg+1]);
-      iarg += 2;       
-    } else if (strcmp(arg[iarg],"relaxbox_rate") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      relaxbox_rate = force->numeric(FLERR,arg[iarg+1]);
-      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;
       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],"eulerexplicit") == 0) integrator = 3;
-     else error->all(FLERR,"Illegal min_modify command");
+      else error->all(FLERR,"Illegal min_modify command");
       iarg += 2;       
    } else if (strcmp(arg[iarg],"line") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");

src/min.h

@@ -58,15 +58,13 @@ class Min : protected Pointers {
   double dmax;                // max dist to move any atom in one step
   int linestyle;              // 0 = backtrack, 1 = quadratic, 2 = forcezero
 
-  int delaystep;              // minium steps of dynamics (adaptglok)
+  // only for min_adaptglok
-  double dt_grow,dt_shrink;   // timestep increase, decrease (adaptglok)
+  int delaystep;              // minium steps of dynamics
-  double alpha0,alpha_shrink; // mixing velocities+forces coefficient (adaptglok)
+  double dtgrow,dtshrink;     // timestep increase, decrease 
-  double tmax,tmin;           // timestep max, min (adaptglok)
+  double alpha0,alphashrink;  // mixing velocities+forces coefficient
-  int integrator;             // choose the style of time integrator (adaptglok)
+  double tmax,tmin;           // timestep multiplicators max, min
-  double relaxbox_modulus;        // Bulk modulus used for box relax (for adaptglok)
+  int integrator;             // Newton integration: euler, leapfrog, verlet...
-  double relaxbox_rate;      // for box relaxation to 0 pressure (for adaptglok)
+  int halfstepback_flag;      // half step backward when v.f <= 0.0
-  int halfstepback_flag;      // 1: half step backward when v.f <= 0.0 (adaptglok)
-  int relaxbox_flag;          // 1: box relaxation iso; 2: axial
 
   int nelist_global,nelist_atom;    // # of PE,virial computes to check
   int nvlist_global,nvlist_atom;

src/min_adaptglok.cpp

@@ -20,10 +20,6 @@
 #include "output.h"
 #include "timer.h"
 #include "error.h"
-#include "variable.h"
-#include "modify.h"
-#include "compute.h"
-#include "domain.h"
 
 using namespace LAMMPS_NS;
 
@@ -41,18 +37,15 @@ void MinAdaptGlok::init()
 {
   Min::init();
 
+  if (tmax < tmin) error->all(FLERR,"tmax cannot be smaller than tmin");
+  if (dtgrow < 1.0) error->all(FLERR,"dtgrow cannot be smaller than 1.0");
+  if (dtshrink > 1.0) error->all(FLERR,"dtshrink cannot be greater than 1.0");
+
   dt = dtinit = update->dt;
   dtmax = tmax * dt;
   dtmin = tmin * dt;
   alpha = alpha0;
-  last_negative = ntimestep_fire = update->ntimestep;
+  last_negative = ntimestep_start = update->ntimestep;
-  
-  if (relaxbox_flag){
-    int icompute = modify->find_compute("thermo_temp");
-    temperature = modify->compute[icompute];
-    icompute = modify->find_compute("thermo_press");
-    pressure = modify->compute[icompute];
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -80,68 +73,6 @@ void MinAdaptGlok::reset_vectors()
   if (nvec) fvec = atom->f[0];
 }
 
-/* ----------------------------------------------------------------------
-   save current box state for converting atoms to lamda coords
-------------------------------------------------------------------------- */
-
-void MinAdaptGlok::save_box_state()
-{
-  boxlo[0] = domain->boxlo[0];
-  boxlo[1] = domain->boxlo[1];
-  boxlo[2] = domain->boxlo[2];
-
-  for (int i = 0; i < 6; i++)
-    h_inv[i] = domain->h_inv[i];
-}
-
-/* ----------------------------------------------------------------------
-   deform the simulation box and remap the particles
-------------------------------------------------------------------------- */
-
-void MinAdaptGlok::relax_box()
-{
-  int i,n;
-      
-  // rescale simulation box from linesearch starting point
-  // scale atom coords for all atoms or only for fix group atoms
-
-  double **x = atom->x;
-  double epsilon;
-  double *current_pressure_v;
-  int *mask = atom->mask;
-  n = atom->nlocal + atom->nghost;
-  save_box_state();
-
-  // convert pertinent atoms and rigid bodies to lamda coords
-
-  domain->x2lamda(n);
-
-  // ensure the virial is tallied
-
-  update->vflag_global = update->ntimestep;
-
-  // compute pressure and change simulation box
-
-  pressure->compute_scalar();
-  pressure->compute_vector();
-  epsilon = pressure->scalar / relaxbox_modulus;
-  current_pressure_v = pressure->vector;
-  for (int i = 0; i < 3; i++) {
-    if (relaxbox_flag == 2) epsilon = current_pressure_v[i] / relaxbox_modulus;
-    domain->boxhi[i] += domain->boxhi[i] * epsilon * relaxbox_rate;;
-  }
-
-  // reset global and local box to new size/shape
-
-  domain->set_global_box();
-  domain->set_local_box();
-
-  // convert atoms and back to box coords
-
-  domain->lamda2x(n);
-  save_box_state();
-}
-
 /* ---------------------------------------------------------------------- */
 
 int MinAdaptGlok::iterate(int maxiter)
@@ -154,13 +85,12 @@ int MinAdaptGlok::iterate(int maxiter)
 
   alpha_final = 0.0;
 
-  double **f = atom->f;
-  double **v = atom->v;
-
   // Leap Frog integration initialization
 
   if (integrator == 2) {
 
+    double **f = atom->f;
+    double **v = atom->v;
     double *rmass = atom->rmass;
     double *mass = atom->mass;
     int *type = atom->type;
@@ -170,6 +100,7 @@ int MinAdaptGlok::iterate(int maxiter)
     neval++;
 
     dtf = 0.5 * dt * force->ftm2v;
+
     if (rmass) {
       for (int i = 0; i < nlocal; i++) {
         dtfm = dtf / rmass[i];
@@ -185,7 +116,6 @@ int MinAdaptGlok::iterate(int maxiter)
         v[i][2] = dtfm * f[i][2];
       }
     }
-
   }
 
   for (int iter = 0; iter < maxiter; iter++) {
@@ -196,21 +126,18 @@ int MinAdaptGlok::iterate(int maxiter)
     ntimestep = ++update->ntimestep;
     niter++;
 
-    // Relax the simulation box
+    // pointers
 
-    if (relaxbox_flag) relax_box();
+    int nlocal = atom->nlocal;
+    double **v = atom->v;
+    double **f = atom->f;
+    double **x = atom->x;
+    double *rmass = atom->rmass;
+    double *mass = atom->mass;
+    int *type = atom->type;
 
    // vdotfall = v dot f
 
-    // Euler || Leap Frog integration
-
-    if (integrator == 0 || integrator == 2) {
-      double **v = atom->v;
-      double **f = atom->f;
-    }
-    int nlocal = atom->nlocal;
-    double **x = atom->x;
-
     vdotf = 0.0;
     for (int i = 0; i < nlocal; i++)
       vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2];
@@ -227,9 +154,10 @@ int MinAdaptGlok::iterate(int maxiter)
     // if (v dot f) > 0:
     // v = (1-alpha) v + alpha |v| Fhat
     // |v| = length of v, Fhat = unit f
-    // The modificatin of v is made wihtin the Verlet integration, after v update
+    // Only: (1-alpha) and alpha |v| Fhat is calculated here
+    // the modificatin of v is made wihtin the integration, after v update
     // if more than delaystep since v dot f was negative:
-    // increase timestep and decrease alpha
+    // increase timestep, update global timestep and decrease alpha
 
     if (vdotfall > 0.0) {
       vdotv = 0.0;
@@ -263,30 +191,29 @@ int MinAdaptGlok::iterate(int maxiter)
       else scale2 = alpha * sqrt(vdotvall/fdotfall);
 
       if (ntimestep - last_negative > delaystep) {
-        dt = MIN(dt*dt_grow,dtmax);
+        dt = MIN(dt*dtgrow,dtmax);
         update->dt = dt;
-        alpha *= alpha_shrink;
+        alpha *= alphashrink;
       }
 
     // else (v dot f) <= 0
     // if more than delaystep since starting the relaxation:
     // reset alpha
-    //    if dt > dtmin:
+    //    if dt*dtshrink > dtmin:
     //    decrease timestep
-    // set x(t) = x(t-0.5*dt)
+    //    update global timestep (for thermo output)
-    // set v = 0
+    // half step back within the dynamics: x(t) = x(t-0.5*dt)
+    // reset velocities: v = 0
 
     } else {
       last_negative = ntimestep;
-      // Limit decrease of timestep
+      if (ntimestep - ntimestep_start > delaystep) {
-      if (ntimestep - ntimestep_fire > delaystep) {
         alpha = alpha0;
-        if (dt > dtmin) {
+        if (dt*dtshrink >= dtmin) {
-          dt *= dt_shrink;
+          dt *= dtshrink;
           update->dt = dt;
         }
       }
-      double **x = atom->x;
       if (halfstepback_flag) {
         for (int i = 0; i < nlocal; i++) {
           x[i][0] -= 0.5 * dtv * v[i][0];
@@ -300,10 +227,6 @@ int MinAdaptGlok::iterate(int maxiter)
 
     // limit timestep so no particle moves further than dmax
 
-    double *rmass = atom->rmass;
-    double *mass = atom->mass;
-    int *type = atom->type;
-
     dtvone = dt;
 
     for (int i = 0; i < nlocal; i++) {
@@ -321,57 +244,20 @@ int MinAdaptGlok::iterate(int maxiter)
       MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld);
     }
 
-    // Semi-implicit Euler integration
+    // Dynamic integration scheme:
+    // 0: semi-implicit Euler 
+    // 1: velocity Verlet
+    // 2: leapfrog (initial half step before the iteration loop)
+    // 3: explicit Euler 
 
-  if (integrator == 0) {
+    // Semi-implicit Euler OR Leap Frog integration
 
-    dtf = dtv * force->ftm2v; 
+    if (integrator == 0 || integrator == 2) {
 
-    if (rmass) {
+      dtf = dtv * force->ftm2v; 
-      for (int i = 0; i < nlocal; i++) {
-        dtfm = dtf / rmass[i];
-        v[i][0] += dtfm * f[i][0];
-        v[i][1] += dtfm * f[i][1];
-        v[i][2] += dtfm * f[i][2];
-        if (vdotfall > 0.0) {
-          v[i][0] = scale1*v[i][0] + scale2*f[i][0];
-          v[i][1] = scale1*v[i][1] + scale2*f[i][1];
-          v[i][2] = scale1*v[i][2] + scale2*f[i][2];
-        }
-        x[i][0] += dtv * v[i][0];
-        x[i][1] += dtv * v[i][1];
-        x[i][2] += dtv * v[i][2];
-      }
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-        dtfm = dtf / mass[type[i]];
-        v[i][0] += dtfm * f[i][0];
-        v[i][1] += dtfm * f[i][1];
-        v[i][2] += dtfm * f[i][2];
-        if (vdotfall > 0.0) {
-          v[i][0] = scale1*v[i][0] + scale2*f[i][0];
-          v[i][1] = scale1*v[i][1] + scale2*f[i][1];
-          v[i][2] = scale1*v[i][2] + scale2*f[i][2];
-        }
-        x[i][0] += dtv * v[i][0];
-        x[i][1] += dtv * v[i][1];
-        x[i][2] += dtv * v[i][2];
-      }
-    }
 
-    eprevious = ecurrent;
+      if (rmass) {
-    ecurrent = energy_force(0);
+        for (int i = 0; i < nlocal; i++) {
-    neval++;
-
-
-  // Velocity Verlet integration
-
-  }else if (integrator == 1) {
-
-    dtf = 0.5 * dtv * force->ftm2v;
-
-    if (rmass) {
-      for (int i = 0; i < nlocal; i++) {
           dtfm = dtf / rmass[i];
           v[i][0] += dtfm * f[i][0];
           v[i][1] += dtfm * f[i][1];
@@ -384,9 +270,9 @@ int MinAdaptGlok::iterate(int maxiter)
           x[i][0] += dtv * v[i][0];
           x[i][1] += dtv * v[i][1];
           x[i][2] += dtv * v[i][2];
-      }
+        }
-    } else {
+      } else {
-      for (int i = 0; i < nlocal; i++) {
+        for (int i = 0; i < nlocal; i++) {
           dtfm = dtf / mass[type[i]];
           v[i][0] += dtfm * f[i][0];
           v[i][1] += dtfm * f[i][1];
@@ -399,123 +285,124 @@ int MinAdaptGlok::iterate(int maxiter)
           x[i][0] += dtv * v[i][0];
           x[i][1] += dtv * v[i][1];
           x[i][2] += dtv * v[i][2];
+        }
       }
-    }
-    
-    eprevious = ecurrent;
-    ecurrent = energy_force(0);
-    neval++;
 
-    if (rmass) {
+      eprevious = ecurrent;
-      for (int i = 0; i < nlocal; i++) {
+      ecurrent = energy_force(0);
+      neval++;
+
+    // Velocity Verlet integration
+
+    }else if (integrator == 1) {
+
+      dtf = 0.5 * dtv * force->ftm2v;
+
+      if (rmass) {
+        for (int i = 0; i < nlocal; i++) {
           dtfm = dtf / rmass[i];
           v[i][0] += dtfm * f[i][0];
           v[i][1] += dtfm * f[i][1];
           v[i][2] += dtfm * f[i][2];
+          if (vdotfall > 0.0) {
+            v[i][0] = scale1*v[i][0] + scale2*f[i][0];
+            v[i][1] = scale1*v[i][1] + scale2*f[i][1];
+            v[i][2] = scale1*v[i][2] + scale2*f[i][2];
+          }
+          x[i][0] += dtv * v[i][0];
+          x[i][1] += dtv * v[i][1];
+          x[i][2] += dtv * v[i][2];
+        }
+      } else {
+        for (int i = 0; i < nlocal; i++) {
+          dtfm = dtf / mass[type[i]];
+          v[i][0] += dtfm * f[i][0];
+          v[i][1] += dtfm * f[i][1];
+          v[i][2] += dtfm * f[i][2];
+          if (vdotfall > 0.0) {
+            v[i][0] = scale1*v[i][0] + scale2*f[i][0];
+            v[i][1] = scale1*v[i][1] + scale2*f[i][1];
+            v[i][2] = scale1*v[i][2] + scale2*f[i][2];
+          }
+          x[i][0] += dtv * v[i][0];
+          x[i][1] += dtv * v[i][1];
+          x[i][2] += dtv * v[i][2];
         }
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-        dtfm = dtf / mass[type[i]];
-        v[i][0] += dtfm * f[i][0];
-        v[i][1] += dtfm * f[i][1];
-        v[i][2] += dtfm * f[i][2];
       }
+      
+      eprevious = ecurrent;
+      ecurrent = energy_force(0);
+      neval++;
+
+      if (rmass) {
+        for (int i = 0; i < nlocal; i++) {
+          dtfm = dtf / rmass[i];
+          v[i][0] += dtfm * f[i][0];
+          v[i][1] += dtfm * f[i][1];
+          v[i][2] += dtfm * f[i][2];
+          }
+      } else {
+        for (int i = 0; i < nlocal; i++) {
+          dtfm = dtf / mass[type[i]];
+          v[i][0] += dtfm * f[i][0];
+          v[i][1] += dtfm * f[i][1];
+          v[i][2] += dtfm * f[i][2];
+        }
+      }
+
+    // Standard Euler integration
+
+    }else if (integrator == 3) {
+
+      dtf = dtv * force->ftm2v; 
+
+      if (rmass) {
+        for (int i = 0; i < nlocal; i++) {
+          dtfm = dtf / rmass[i];
+          if (vdotfall > 0.0) {
+            v[i][0] = scale1*v[i][0] + scale2*f[i][0];
+            v[i][1] = scale1*v[i][1] + scale2*f[i][1];
+            v[i][2] = scale1*v[i][2] + scale2*f[i][2];
+          }
+          x[i][0] += dtv * v[i][0];
+          x[i][1] += dtv * v[i][1];
+          x[i][2] += dtv * v[i][2];
+          v[i][0] += dtfm * f[i][0];
+          v[i][1] += dtfm * f[i][1];
+          v[i][2] += dtfm * f[i][2];
+        }
+      } else {
+        for (int i = 0; i < nlocal; i++) {
+          dtfm = dtf / mass[type[i]];
+          if (vdotfall > 0.0) {
+            v[i][0] = scale1*v[i][0] + scale2*f[i][0];
+            v[i][1] = scale1*v[i][1] + scale2*f[i][1];
+            v[i][2] = scale1*v[i][2] + scale2*f[i][2];
+          }
+          x[i][0] += dtv * v[i][0];
+          x[i][1] += dtv * v[i][1];
+          x[i][2] += dtv * v[i][2];
+          v[i][0] += dtfm * f[i][0];
+          v[i][1] += dtfm * f[i][1];
+          v[i][2] += dtfm * f[i][2];
+        }
+      }
+
+      eprevious = ecurrent;
+      ecurrent = energy_force(0);
+      neval++;
+
     }
 
-  // Leap Frog integration
-
-  }else if (integrator == 2) {
-
-    dtf = dtv * force->ftm2v; 
-
-    if (rmass) {
-      for (int i = 0; i < nlocal; i++) {
-        dtfm = dtf / rmass[i];
-        v[i][0] += dtfm * f[i][0];
-        v[i][1] += dtfm * f[i][1];
-        v[i][2] += dtfm * f[i][2];
-        if (vdotfall > 0.0) {
-          v[i][0] = scale1*v[i][0] + scale2*f[i][0];
-          v[i][1] = scale1*v[i][1] + scale2*f[i][1];
-          v[i][2] = scale1*v[i][2] + scale2*f[i][2];
-        }
-        x[i][0] += dtv * v[i][0];
-        x[i][1] += dtv * v[i][1];
-        x[i][2] += dtv * v[i][2];
-      }
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-        dtfm = dtf / mass[type[i]];
-        v[i][0] += dtfm * f[i][0];
-        v[i][1] += dtfm * f[i][1];
-        v[i][2] += dtfm * f[i][2];
-        if (vdotfall > 0.0) {
-          v[i][0] = scale1*v[i][0] + scale2*f[i][0];
-          v[i][1] = scale1*v[i][1] + scale2*f[i][1];
-          v[i][2] = scale1*v[i][2] + scale2*f[i][2];
-        }
-        x[i][0] += dtv * v[i][0];
-        x[i][1] += dtv * v[i][1];
-        x[i][2] += dtv * v[i][2];
-      }
-    }
-
-    eprevious = ecurrent;
-    ecurrent = energy_force(0);
-    neval++;
-
-  // Standard Euler integration
-
-  }else if (integrator == 3) {
-
-    dtf = dtv * force->ftm2v; 
-
-    if (rmass) {
-      for (int i = 0; i < nlocal; i++) {
-        dtfm = dtf / rmass[i];
-        if (vdotfall > 0.0) {
-          v[i][0] = scale1*v[i][0] + scale2*f[i][0];
-          v[i][1] = scale1*v[i][1] + scale2*f[i][1];
-          v[i][2] = scale1*v[i][2] + scale2*f[i][2];
-        }
-        x[i][0] += dtv * v[i][0];
-        x[i][1] += dtv * v[i][1];
-        x[i][2] += dtv * v[i][2];
-        v[i][0] += dtfm * f[i][0];
-        v[i][1] += dtfm * f[i][1];
-        v[i][2] += dtfm * f[i][2];
-      }
-    } else {
-      for (int i = 0; i < nlocal; i++) {
-        dtfm = dtf / mass[type[i]];
-        if (vdotfall > 0.0) {
-          v[i][0] = scale1*v[i][0] + scale2*f[i][0];
-          v[i][1] = scale1*v[i][1] + scale2*f[i][1];
-          v[i][2] = scale1*v[i][2] + scale2*f[i][2];
-        }
-        x[i][0] += dtv * v[i][0];
-        x[i][1] += dtv * v[i][1];
-        x[i][2] += dtv * v[i][2];
-        v[i][0] += dtfm * f[i][0];
-        v[i][1] += dtfm * f[i][1];
-        v[i][2] += dtfm * f[i][2];
-      }
-    }
-
-    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
+    // reset the timestep to the initial value
 
     if (update->etol > 0.0 && ntimestep-last_negative > delaystep) {
       if (update->multireplica == 0) {
         if (fabs(ecurrent-eprevious) <
-            update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)) {
+            update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)){
           update->dt = dtinit;
           return ETOL;
         }
@@ -534,6 +421,7 @@ int MinAdaptGlok::iterate(int maxiter)
 
     // force tolerance criterion
     // sync across replicas if running multi-replica minimization
+    // reset the timestep to the initial value
 
     if (update->ftol > 0.0) {
       fdotf = fnorm_sqr();
@@ -562,6 +450,8 @@ int MinAdaptGlok::iterate(int maxiter)
     }
   }
 
+  // reset the timestep to the initial value
+
   update->dt = dtinit;
   return MAXITER;
 }

src/min_adaptglok.h

@@ -31,16 +31,12 @@ class MinAdaptGlok : public Min {
   void init();
   void setup_style();
   void reset_vectors();
-  void save_box_state();
-  void relax_box();
   int iterate(int);
 
  private:
   double dt,dtmax,dtmin,dtinit;
   double alpha;
-  bigint last_negative,ntimestep_fire;
+  bigint last_negative,ntimestep_start;
-  class Compute *temperature,*pressure;
-  double boxlo[3],h_inv[6];
 };
 
 }