Add relax_box function. Adapted documentation.

doc/src/min_modify.txt

@@ -28,7 +28,10 @@ keyword = {dmax} or {line} or {delaystep} or {dt_grow} or {dt_shrink}
   {tmin} value = {float}  for adaptglok
   {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} 
                        or {leapfrog}  for adaptglok
-  {halfstepback} value = {yes} or {no}  for adaptglok :pre
+  {halfstepback} value = {yes} or {no} for adaptglok :pre
+  {relaxbox} value = {no} or {iso} or {axial} for adaptglok 
+  {relaxbox_modulus} value = {float} for adaptglok 
+  {relaxbox_rate} value = {float} for adaptglok 
 :ule
 
 [Examples:]
@@ -84,6 +87,12 @@ to adapt them for particular cases.
 In particular for debugging purpose, it is possible to switch off the 
 the backstep when downhill, by using and {halfstepback}.
 
+The style {adaptglok} allows for a basic box relaxation to a stress-free 
+state by choosing {relaxbox} {iso} or {axial}. One should give an
+approximated value for the bulk modulus {relaxbox_modulus} of the system.
+The default value of the relaxation rate {relaxbox_rate} should be 
+optimum in most situations.
+
 [Restrictions:] none
 
 [Related commands:]
@@ -94,4 +103,5 @@ the backstep when downhill, by using and {halfstepback}.
 
 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, tmin = 0.02, adaptstep = yes and halfstepback = yes.
+tmax = 2.0, tmin = 0.02,  adaptstep = yes, halfstepback = yes,
+relaxbox = no, relaxbox_modulus = 1000000 and relaxbox_rate = 0.33.

doc/src/min_style.txt

@@ -66,6 +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
 nudged elastic band (NEB) calculations via the "neb"_neb.html command.

src/min.cpp

@@ -63,7 +63,10 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
   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;
@@ -687,6 +690,21 @@ 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;

src/min.h

@@ -63,7 +63,10 @@ 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)
+  double relaxbox_modulus;        // Bulk modulus used for box relax (for adaptglok)
+  double relaxbox_rate;      // for box relaxation to 0 pressure (for adaptglok)
   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,6 +20,10 @@
 #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;
 
@@ -42,6 +46,13 @@ void MinAdaptGlok::init()
   dtmin = tmin * dt;
   alpha = alpha0;
   last_negative = ntimestep_fire = 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];
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -69,6 +80,68 @@ 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)
@@ -123,7 +196,11 @@ int MinAdaptGlok::iterate(int maxiter)
     ntimestep = ++update->ntimestep;
     niter++;
 
-    // vdotfall = v dot f
+    // Relax the simulation box
+
+    if (relaxbox_flag) relax_box();
+
+   // vdotfall = v dot f
 
     // Euler || Leap Frog integration
 
@@ -132,6 +209,7 @@ int MinAdaptGlok::iterate(int maxiter)
       double **f = atom->f;
     }
     int nlocal = atom->nlocal;
+    double **x = atom->x;
 
     vdotf = 0.0;
     for (int i = 0; i < nlocal; i++)
@@ -243,9 +321,7 @@ int MinAdaptGlok::iterate(int maxiter)
       MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld);
     }
 
-  double **x = atom->x;
-
-  // Semi-implicit Euler integration
+    // Semi-implicit Euler integration
 
   if (integrator == 0) {
 

src/min_adaptglok.h

@@ -31,12 +31,16 @@ 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;
+  class Compute *temperature,*pressure;
+  double boxlo[3],h_inv[6];
 };
 
 }