New stopping criterion based on v.f number. Timestep reset moved to min.*

src/min.cpp

@@ -65,6 +65,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
   integrator = 0;
   halfstepback_flag = 1;
   delaystep_start_flag = 1;
+  max_vdotf_negatif = 2000;
   relaxbox_mod = 1000000;
   relaxbox_rate = 0.33;
   relaxbox_flag = 0;
@@ -190,6 +191,10 @@ void Min::init()
   neighbor->dist_check = 1;
 
   niter = neval = 0;
+
+  // store timestep size (important for variable timestep minimizer)
+
+  dtinit = update->dt;
 }
 
 /* ----------------------------------------------------------------------
@@ -477,6 +482,10 @@ void Min::cleanup()
 
   modify->delete_fix("MINIMIZE");
   domain->box_too_small_check();
+
+  // reset timestep size (important for variable timestep minimizer)
+
+  update->dt = dtinit;
 }
 
 /* ----------------------------------------------------------------------
@@ -699,6 +708,10 @@ void Min::modify_params(int narg, char **arg)
       else if (strcmp(arg[iarg+1],"no") == 0) delaystep_start_flag = 0;
       else error->all(FLERR,"Illegal min_modify command");
       iarg += 2;       
+    } else if (strcmp(arg[iarg],"vdfmax") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
+      max_vdotf_negatif = 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;
@@ -924,6 +937,7 @@ char *Min::stopstrings(int n)
                            "quadratic factors are zero",
                            "trust region too small",
                            "HFTN minimizer error",
-                           "walltime limit reached"};
+                           "walltime limit reached",
+                           "max iterations with v.f negative"};
   return (char *) strings[n];
 }

src/min.h

@@ -48,7 +48,7 @@ class Min : protected Pointers {
 
   // possible return values of iterate() method
   enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,
-       ZEROQUAD,TRSMALL,INTERROR,TIMEOUT};
+       ZEROQUAD,TRSMALL,INTERROR,TIMEOUT,MAXVDOTF};
 
  protected:
   int eflag,vflag;            // flags for energy/virial computation
@@ -57,6 +57,7 @@ class Min : protected Pointers {
 
   double dmax;                // max dist to move any atom in one step
   int linestyle;              // 0 = backtrack, 1 = quadratic, 2 = forcezero
+  double dtinit;              // store the default timestep
 
   // only for minimize style adaptglok
   int delaystep;              // minium steps of dynamics
@@ -66,6 +67,7 @@ class Min : protected Pointers {
   int integrator;             // Newton integration: euler, leapfrog, verlet...
   int halfstepback_flag;      // half step backward when v.f <= 0.0
   int delaystep_start_flag;   // delay the initial dt_shrink
+  int max_vdotf_negatif;      // maximum iteration with v.f > 0.0
   double relaxbox_mod;        // Bulk modulus used for box relax
   double relaxbox_rate;       // for box relaxation to 0 pressure
   int relaxbox_flag;          // 1: box relaxation iso; 2: aniso

src/min_adaptglok.cpp

@@ -59,11 +59,12 @@ void MinAdaptGlok::init()
   if (p_flag[2] && domain->zperiodic == 0)
     error->all(FLERR,"Cannot use boxrelax on a non-periodic dimension");
 
-  dt = dtinit = update->dt;
+  dt = update->dt;
   dtmax = tmax * dt;
   dtmin = tmin * dt;
   alpha = alpha0;
   last_negative = ntimestep_start = update->ntimestep;
+  vdotf_negatif = 0;
 
   if (relaxbox_flag){
     int icompute = modify->find_compute("thermo_press");
@@ -269,6 +270,7 @@ int MinAdaptGlok::iterate(int maxiter)
 
     if (vdotfall > 0.0) {
       vdotv = 0.0;
+      vdotf_negatif = 0;
       for (int i = 0; i < nlocal; i++)
         vdotv += v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
       MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world);
@@ -325,6 +327,15 @@ int MinAdaptGlok::iterate(int maxiter)
           update->dt = dt;
         }
       }
+
+      // stopping criterion while stuck in a local bassin
+
+      vdotf_negatif++;
+      if (max_vdotf_negatif > 0 && vdotf_negatif > max_vdotf_negatif)
+        return MAXVDOTF;
+
+      // inertia correction
+      
       if (halfstepback_flag) {
         for (int i = 0; i < nlocal; i++) {
           x[i][0] -= 0.5 * dtv * v[i][0];
@@ -529,10 +540,8 @@ int MinAdaptGlok::iterate(int maxiter)
       if (update->multireplica == 0) {
         if (fabs(ecurrent-eprevious) <
             update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)
-            && pflag) {
-          update->dt = dtinit;
+            && pflag)
           return ETOL;
-        }
       } else {
         if (fabs(ecurrent-eprevious) <
             update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)
@@ -540,10 +549,8 @@ int MinAdaptGlok::iterate(int maxiter)
           flag = 0;
         else flag = 1;
         MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
-        if (flagall == 0) {
-          update->dt = dtinit;
+        if (flagall == 0)
           return ETOL;
-        }
       }
     }
 
@@ -554,18 +561,14 @@ int MinAdaptGlok::iterate(int maxiter)
     if (update->ftol > 0.0) {
       fdotf = fnorm_sqr();
       if (update->multireplica == 0) {
-        if (fdotf < update->ftol*update->ftol && pflag) {
-          update->dt = dtinit;
+        if (fdotf < update->ftol*update->ftol && pflag)
           return FTOL;
-        }
       } else {
         if (fdotf < update->ftol*update->ftol && pflag) flag = 0;
         else flag = 1;
         MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
-        if (flagall == 0) {
-          update->dt = dtinit;
+        if (flagall == 0)
           return FTOL;
-        }
       }
     }
 
@@ -578,8 +581,5 @@ int MinAdaptGlok::iterate(int maxiter)
     }
   }
 
-  // reset the timestep to the initial value
-
-  update->dt = dtinit;
   return MAXITER;
 }

src/min_adaptglok.h

@@ -36,10 +36,10 @@ class MinAdaptGlok : public Min {
   int iterate(int);
 
  private:
-  double dt,dtmax,dtmin,dtinit;
+  double dt,dtmax,dtmin;
   double alpha;
   bigint last_negative,ntimestep_start;
-  int pflag;
+  int pflag,vdotf_negatif;
   class Compute *temperature,*pressure;
   double boxlo[3],h_inv[6];
 };