Remove option relaxbox from adaptglok: wrong behavior with non-P boundaries. Code cleanup.

doc/src/min_modify.txt

@@ -26,8 +26,7 @@ keyword = {dmax} or {line} or {alpha_damp} or {discrete_factor}
 
 these keywords apply only to the "min_style"_min_style.html {adaptglok} :ulb,l
 keyword = {integrator} or {tmax} or {tmin} or {delaystep} or {dtgrow} or {dtshrink} 
-          or {alpha0} or {alphashrink} or {relaxbox} or {relaxbox_mod} or {relaxbox_rate} 
+          or {alpha0} or {alphashrink} or {halfstepback} or {initialdelay} or {vdfmax}
-          or {ptol} or {halfstepback} or {initialdelay} or {vdfmax}
   {integrator} arg = {eulerimplicit} or {eulerexplicit} or {verlet} or {leapfrog}
     eulerimplicit,eulerexplicit,verlet,leapfrog = integration scheme
   {tmax} arg = coef
@@ -44,15 +43,6 @@ keyword = {integrator} or {tmax} or {tmin} or {delaystep} or {dtgrow} or {dtshri
     alpha = coefficient mixing velocities and forces
   {alphashrink} arg = shrink
     shrink = factor decreasing alpha
-  {relaxbox} arg = {no} or {iso} or {aniso} or {x} or {y} or {z}
-    no = no changes in box dimension (default)
-    iso, aniso, x, y, z = type of box relaxation
-  {relaxbox_mod} arg = modulus
-    modulus = bulk modulus of the system, order of magnitude (pressure units)
-  {relaxbox_rate} arg = rate
-    rate = scaling factor to relax the box dimensions
-  {ptol} arg = pressure
-    pressure = threshold below which the box pressure is considered as null (pressure units)
   {vdfmax} arg = max 
     max = maximum number of consecutive iterations with P(t) < 0 before forced interruption 
   {halfstepback} arg = {yes} or {no}
@@ -62,7 +52,7 @@ keyword = {integrator} or {tmax} or {tmin} or {delaystep} or {dtgrow} or {dtshri
 [Examples:]
 
 min_modify dmax 0.2
-min_modify integrator verlet tmax 0.4 relaxbox y relaxbox_mod 2e7 :pre
+min_modify integrator verlet tmax 0.4 :pre
 
 [Description:]
 
@@ -127,25 +117,6 @@ happen when the system comes to be stuck in a local basin of the phase space.
 For debugging purposes, it is possible to switch off the inertia correction 
 ({halfstepback} = {no}) and the initial delay ({initialdelay} = {no}).
 
-The style {adaptglok} performing a damped dynamics, it is not possible to use
-"fix box/relax"_fix_box_relax.html to relax the simulation box.
-Thus {adaptglok} implements a rudimentary box relaxation procedure that can be
-controlled by the keywords {relaxbox}, {relaxbox_mod}, {relaxbox_rate} 
-and {ptol}.
-The argument {relaxbox} control in which directions the pressure is relaxed.
-Note that the corresponding directions have to be periodic.
-Note also that {relaxbox_mod} doesn't requires the exact value for the bulk modulus,
-but rather the order of magnitude (in pressure "units"_units.html).
-Ultimately, {relaxbox_mod} and {relaxbox_rate} 
-control how fast the relaxation of the box is performed: lower values will slow
-down the box relaxation but improve the stability of the procedure.
-
-NOTE: the option {relaxbox} is currently experimental and often 
-requires to tune the communication cutoff for ghost atoms with the command 
-"comm_modify cutoff"_comm_modify.html. The value will depend on the expected
-box variation and the number of cpu. A value from 2 to 3 times the current cutoff
-(largest pair cutoff + neighbor skin) is often enough.
-
 Keywords {alpha_damp} and {discrete_factor} only make sense when
 a "min_spin"_min_spin.html command is declared. 
 Keyword {alpha_damp} defines an analog of a magnetic Gilbert
@@ -168,5 +139,4 @@ Default values are {alpha_damp} = 1.0 and {discrete_factor} = 10.0.
 The option defaults are dmax = 0.1, line = quadratic, 
 integrator = eulerimplicit, tmax = 10.0, tmin = 0.02,
 delaystep = 20, dtgrow = 1.1, dtshrink = 0.5, alpha0 = 0.25, alphashrink = 0.99,
-relaxbox = no, relaxbox_mod = 1e6 and relaxbox_rate = 0.33, ptol = 0.1
 vdfmax = 2000, halfstepback = yes and initialdelay = yes.

doc/src/min_style.txt

@@ -64,7 +64,7 @@ a minimization.
 
 Style {adaptglok} is a re-implementation of the style {fire} with
 additional optimizations of the method described
-in "(Bitzek)"_#Bitzek.
+in "(Bitzek)"_#Bitzek, including different time integration schemes.
 
 Style {spin} is a damped spin dynamics with an adaptive 
 timestep.
@@ -81,8 +81,6 @@ would normally use for dynamics simulations.
 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. 
-The {adaptglok} style actually support box relaxation by the implementation of
-a basic relaxation scheme, see "min_modify"_min_modify.html.
 
 [Restrictions:] none
 

doc/src/minimize.txt

@@ -64,8 +64,7 @@ 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. The style
+{adaptglok} perform damped dynamics using an Euler integration step.
-{adaptglok} can also use a leapfrog or velocity Verlet integration step.
 Thus they require a "timestep"_timestep.html be defined.
 
 NOTE: The damped dynamic minimizers use whatever timestep you have

src/min.cpp

@@ -66,11 +66,6 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
   halfstepback_flag = 1;
   delaystep_start_flag = 1;
   max_vdotf_negatif = 2000;
-  relaxbox_mod = 1000000;
-  relaxbox_rate = 0.33;
-  relaxbox_flag = 0;
-  ptol = 0.1;
-  p_flag[0] = p_flag[1] = p_flag[2] = 0;
 
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
@@ -724,43 +719,6 @@ void Min::modify_params(int narg, char **arg)
       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");
-      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;
-        p_flag[0] = p_flag[1] = p_flag[2] = 1;
-        if (dimension == 2) p_flag[2] = 0;
-      } else if (strcmp(arg[iarg+1],"aniso") == 0) {
-        relaxbox_flag = 2;
-        p_flag[0] = p_flag[1] = p_flag[2] = 1;
-        if (dimension == 2) p_flag[2] = 0;
-      } else if (strcmp(arg[iarg+1],"x") == 0) {
-        relaxbox_flag = 2;
-        p_flag[0] = 1;
-      } else if (strcmp(arg[iarg+1],"y") == 0) {
-        relaxbox_flag = 2;
-        p_flag[1] = 1;
-      } else if (strcmp(arg[iarg+1],"z") == 0) {
-        relaxbox_flag = 2;
-        p_flag[2] = 1;
-        if (dimension == 2)
-          error->all(FLERR,"Invalid min_modify command for a 2d simulation");
-      } else error->all(FLERR,"Illegal min_modify command");
-      iarg += 2;       
-    } else if (strcmp(arg[iarg],"relaxbox_mod") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      relaxbox_mod = 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],"ptol") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
-      ptol = force->numeric(FLERR,arg[iarg+1]);
       iarg += 2;        
     } else if (strcmp(arg[iarg],"line") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");

src/min.h

@@ -69,12 +69,6 @@ class Min : protected Pointers {
   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
-  double ptol;                // pressure threshold for boxrelax
-  int p_flag[3];
-  int dimension;
 
   int nelist_global,nelist_atom;    // # of PE,virial computes to check
   int nvlist_global,nvlist_atom;

src/min_adaptglok.cpp

@@ -11,7 +11,7 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#include <math.h>
+#include <cmath>
 #include "min_adaptglok.h"
 #include "universe.h"
 #include "atom.h"
@@ -48,17 +48,6 @@ void MinAdaptGlok::init()
   if (tmax < tmin) error->all(FLERR,"tmax has to be larger than tmin");
   if (dtgrow < 1.0) error->all(FLERR,"dtgrow has to be larger than 1.0");
   if (dtshrink > 1.0) error->all(FLERR,"dtshrink has to be smaller than 1.0");
-  if (relaxbox_mod < 0.0) error->all(FLERR,"relaxbox_mod has to be positif");
-  if (relaxbox_rate < 0.0 || relaxbox_rate > 1.0) error->all(FLERR,"relaxbox_rate has to be positif, lower than 1.0");
-
-  // require periodicity in boxrelax dimensions
-
-  if (p_flag[0] && domain->xperiodic == 0)
-    error->all(FLERR,"Cannot use boxrelax on a non-periodic dimension");
-  if (p_flag[1] && domain->yperiodic == 0)
-    error->all(FLERR,"Cannot use boxrelax on a non-periodic dimension");
-  if (p_flag[2] && domain->zperiodic == 0)
-    error->all(FLERR,"Cannot use boxrelax on a non-periodic dimension");
 
   dt = update->dt;
   dtmax = tmax * dt;
@@ -67,18 +56,6 @@ void MinAdaptGlok::init()
   last_negative = ntimestep_start = update->ntimestep;
   vdotf_negatif = 0;
 
-  if (relaxbox_flag){
-
-    // require the box to be orthogonal
-
-    if (domain->triclinic)
-      error->all(FLERR,"Cannot use boxrelax with triclinic box");
-
-    int icompute = modify->find_compute("thermo_press");
-    pressure = modify->compute[icompute];
-
-    pflag = 1;
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -92,16 +69,14 @@ void MinAdaptGlok::setup_style()
 
   const char *s1[] = {"eulerimplicit","verlet","leapfrog","eulerexplicit"};
   const char *s2[] = {"no","yes"};
-  const char *s3[] = {"no","iso","aniso"};
 
   if (comm->me == 0 && logfile) {
       fprintf(logfile,"  Parameters for adaptglok: \n"
       "    dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin "
       "   integrator halfstepback relaxbox relaxbox_mod relaxbox_rate ptol \n"
-      "    %4g %9i %6g %8g %6g %11g %4g %4g %13s %12s %8s %12g %13g %4g \n",
+      "    %4g %9i %6g %8g %6g %11g %4g %4g %13s %12s \n",
       dmax, delaystep, dtgrow, dtshrink, alpha0, alphashrink, tmax, tmin, 
-      s1[integrator], s2[halfstepback_flag], s3[relaxbox_flag], relaxbox_mod,
+      s1[integrator], s2[halfstepback_flag]);
-      relaxbox_rate, ptol);
   }
 
   // initialize the velocities
@@ -124,98 +99,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()
-{
-  // rescale simulation box and scale atom coords for all atoms
-  // inspired by change_box.cpp
-
-  int i,n;
-  double **x = atom->x;
-  double **v = atom->v;
-  double epsilon,disp;
-  int nlocal = atom->nlocal;
-
-  domain->pbc();
-  save_box_state();
-  neighbor->setup_bins();
-  comm->exchange();
-  comm->borders();
-  if (neighbor->decide()) neighbor->build(1);
-
-  // ensure the virial is tallied, update the flag
-
-  pressure->addstep(update->ntimestep);
-  update->vflag_global = update->ntimestep;
-
-  // Only when the presure is not yet free:
-  // - compute and apply box re-scaling
-  // - freez atoms
-
-  if (pflag != 1){
-
-    // compute pressure and change simulation box
-
-    pressure->compute_scalar();
-    pressure->compute_vector();
-    epsilon = pressure->scalar / relaxbox_mod;
-    for (int i = 0; i < 3; i++) {
-      if (relaxbox_flag == 2) epsilon = pressure->vector[i] / relaxbox_mod;
-      disp = domain->h[i] * epsilon * relaxbox_rate;
-      if (fabs(disp) > dmax) disp > 0.0 ? disp = dmax : disp = -1 * dmax;
-      domain->boxlo[i] -= p_flag[i] * disp * 0.5;
-      domain->boxhi[i] += p_flag[i] * disp * 0.5;
-    }
-
-    // reset global and local box to new size/shape
-
-    domain->set_initial_box();
-    domain->set_global_box();
-    domain->set_local_box();
-
-    // convert atoms to lamda coords, using last box state
-    // convert atoms back to box coords, using current box state
-    // save current box state
-    
-    for (i = 0; i < nlocal; i++)
-      domain->x2lamda(x[i],x[i],boxlo,h_inv);
-    for (i = 0; i < nlocal; i++)
-      domain->lamda2x(x[i],x[i]);
-    save_box_state();
-
-    // move atoms back inside simulation box and to new processors
-    // use remap() instead of pbc()
-    //   in case box moved a long distance relative to atoms
-
-    imageint *image = atom->image;
-    for (i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
-    domain->reset_box();
-
-    // freez atoms velocities when the box is rescaled
-    // prevent atoms getting unintended extra velocity
-
-    for (int i = 0; i < nlocal; i++)
-      v[i][0] = v[i][1] = v[i][2] = 0.0;
-  }
-}
-
 /* ---------------------------------------------------------------------- */
 
 int MinAdaptGlok::iterate(int maxiter)
@@ -269,10 +152,6 @@ int MinAdaptGlok::iterate(int maxiter)
     ntimestep = ++update->ntimestep;
     niter++;
 
-    // Relax the simulation box
-
-    if (relaxbox_flag) relax_box();
-
     // pointers
 
     int nlocal = atom->nlocal;
@@ -367,10 +246,9 @@ int MinAdaptGlok::iterate(int maxiter)
       }
 
       // stopping criterion while stuck in a local bassin of the PES
-      // only checked when the box dimesions are not modified bu relax_box()
 
       vdotf_negatif++;
-      if (pflag == 1 && max_vdotf_negatif > 0 && vdotf_negatif > max_vdotf_negatif)
+      if (max_vdotf_negatif > 0 && vdotf_negatif > max_vdotf_negatif)
         return MAXVDOTF;
 
       // inertia correction
@@ -555,21 +433,6 @@ int MinAdaptGlok::iterate(int maxiter)
 
     }
 
-    // Pressure relaxation criterion
-    // set pflag = 0 if relaxbox is activated and pressure > ptol
-    // pflag = 0 will hinder the energy or force criterion
-
-    pflag = 1;
-    if (relaxbox_flag == 1){
-      pressure->compute_scalar();
-      if (fabs(pressure->scalar) > ptol) pflag = 0;
-      
-    }else if (relaxbox_flag == 2){
-      pressure->compute_vector();
-      for (int i = 0; i < 3; i++)
-        if (fabs(pressure->vector[i]) * p_flag[i] > ptol) pflag = 0;
-    }
-
     // energy tolerance criterion
     // only check after delaystep elapsed since velocties reset to 0
     // sync across replicas if running multi-replica minimization
@@ -578,13 +441,11 @@ int MinAdaptGlok::iterate(int maxiter)
     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))
-            && pflag)
           return ETOL;
       } else {
         if (fabs(ecurrent-eprevious) <
-            update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)
+            update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
-            && pflag)
           flag = 0;
         else flag = 1;
         MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
@@ -600,10 +461,10 @@ int MinAdaptGlok::iterate(int maxiter)
     if (update->ftol > 0.0) {
       fdotf = fnorm_sqr();
       if (update->multireplica == 0) {
-        if (fdotf < update->ftol*update->ftol && pflag)
+        if (fdotf < update->ftol*update->ftol)
           return FTOL;
       } else {
-        if (fdotf < update->ftol*update->ftol && pflag) flag = 0;
+        if (fdotf < update->ftol*update->ftol) flag = 0;
         else flag = 1;
         MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
         if (flagall == 0)

src/min_adaptglok.h

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