Final polish of boxrelax. Updated documentation.

doc/src/min_modify.txt

@@ -13,35 +13,45 @@ min_modify command :h3
 min_modify keyword values ... :pre
 
 one or more keyword/value pairs may be listed :ulb,l
-keyword = {dmax} or {line} or {delaystep} or {dt_grow} or {dt_shrink}
+keyword = {dmax} or {line} 
-          or {alpha0} or {alpha_shrink} or {tmax} or {tmin} or {integrator}
+          or {integrator} or {tmax} or {tmin} 
+          or {delaystep} or {dt_grow} or {dt_shrink} or {alpha0} or {alpha_shrink}
+          or {relaxbox} or {relaxbox_mod} or {relaxbox_rate} or {ptol}
+          or {halfstepback}
   {dmax} value = max
     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 = steps
+  {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} or {leapfrog}
-    steps = number of steps for dynamics temporization, adaptglok only
+    eulerimplicit,eulerexplicit,verlet,leapfrog = integration scheme (adaptglok)
-  {dtgrow} value = grow
-    grow = factor increasing the timestep, adaptglok only
-  {dtshrink} value = shrink
-    shrink = factor decreasing the timestep, adaptglok only
-  {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
+    coef = factor defining the maximum timestep (adaptglok)
   {tmin} value = coef
-    coef = factor defining the minimum timestep, adaptglok only
+    coef = factor defining the minimum timestep (adaptglok)
-  {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} 
+  {delaystep} value = steps
-                       or {leapfrog}
+    steps = number of steps for dynamics temporization (adaptglok)
-    eulerimplicit,eulerexplicit,verlet,leapfrog = integration scheme, adaptglok
+  {dtgrow} value = grow
-  {halfstepback} value = {yes} or {no}, adaptglok only :pre
+    grow = factor increasing the timestep (adaptglok)
+  {dtshrink} value = shrink
+    shrink = factor decreasing the timestep (adaptglok)
+  {alpha0} value = alpha
+    alpha = coefficient mixing velocities and forces (adaptglok)
+  {alphashrink} value = shrink
+    shrink = factor decreasing alpha (adaptglok)
+  {relaxbox} value = {no} or {iso} or {axial} or {x} or {y} or {z}
+    no = no changes in box dimension (default)
+    iso,aniso,x,y,z = type of box relaxation  (adaptglok)
+  {relaxbox_mod} value = modulus
+    modulus = bulk modulus of the system, order of magnitude (adaptglok)
+  {relaxbox_rate} value = rate
+    rate = scaling factor to relax the box dimensions (adaptglok) 
+  {halfstepback} value = {yes} or {no} (adaptglok) :pre
 :ule
 
 [Examples:]
 
-min_modify dmax 0.2 :pre
+min_modify dmax 0.2
+min_modify integrator verlet tmax 0.4 relaxbox y relaxbox_mod 2e7 :pre
 
 [Description:]
 
@@ -85,21 +95,34 @@ 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}, {dtgrow}, {dtshrink}, 
+to optimize the relaxation: {integrator}, {tmax}, {tmin}, {delaystep}, {dtgrow}, {dtshrink}, 
-{alpha0}, {alphashrink}, {tmax} and {tmin}.
+{alpha0}, and {alphashrink}.
 Different Newton {integrator} can be selected: explicit Euler, 
 semi-implicit Euler (= symplectic Euler), leapfrog, and velocity Verlet.
 The parameters {tmax} and {tmin} define the maximum and minimum timestep 
 allowed during an adaptglok minimization. Those are not in time unit, but are 
-multiplication factor applied to the "timestep"_timestep.html. For example, 
+multiplication factor applied to the "timestep"_timestep.html. Thus 
-{tmax} = 2.0 means that the maximum value the timestep can reach is twice the 
+{tmax} = 2.0 in metal "units"_units.html means that the maximum value
-value defined by "timestep"_timestep.html.
+the timestep can reached during the relaxation is 2 ps (with the default
+"timestep"_timestep.html value).
 Note that even with the default parameters being chosen to be reliable in most
 cases, adjusting "timestep"_timestep.html, {tmax} and {tmin} should be consider
 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.
 
+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}.
+Note 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.
+
+
 [Restrictions:] none
 
 [Related commands:]
@@ -109,6 +132,7 @@ downhill with {halfstepback} = 0.
 [Default:]
 
 The option defaults are dmax = 0.1, line = quadratic, 
-integrator = eulerimplicit, delaystep = 20, dt_grow = 1.1, dt_shrink = 0.5, 
+integrator = eulerimplicit, tmax = 2.0, tmin = 0.02,
-alpha0 = 0.25, alpha_shrink = 0.99, tmax = 2.0, tmin = 0.02,  adaptstep = yes 
+delaystep = 20, dt_grow = 1.1, dt_shrink = 0.5, alpha0 = 0.25, alpha_shrink = 0.99,
+relaxbox = no, relaxbox_mod = 1e6 and relaxbox_rate = 0.33, ptol = 0.1
 and halfstepback = yes.

doc/src/min_style.txt

@@ -78,6 +78,8 @@ 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
 

src/min.cpp

@@ -67,7 +67,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
   relaxbox_mod = 1000000;
   relaxbox_rate = 0.33;
   relaxbox_flag = 0;
-  ptol = 1e-5;
+  ptol = 0.1;
   p_flag[0] = p_flag[1] = p_flag[2] = 0;
 
   elist_global = elist_atom = NULL;

src/min_adaptglok.cpp

@@ -24,6 +24,8 @@
 #include "modify.h"
 #include "compute.h"
 #include "domain.h"
+#include "neighbor.h"
+#include "comm.h"
 
 using namespace LAMMPS_NS;
 
@@ -46,8 +48,7 @@ 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) 
+  if (relaxbox_mod < 0.0) error->all(FLERR,"relaxbox_mod has to be positif");
-      error->all(FLERR,"relaxbox_mod has to be positif");
 
   // require periodicity in boxrelax dimensions
 
@@ -65,9 +66,7 @@ void MinAdaptGlok::init()
   last_negative = ntimestep_start = update->ntimestep;
 
   if (relaxbox_flag){
-  int icompute = modify->find_compute("thermo_temp");
+    int icompute = modify->find_compute("thermo_press");
-  temperature = modify->compute[icompute];
-  icompute = modify->find_compute("thermo_press");
     pressure = modify->compute[icompute];
   }
 }
@@ -79,6 +78,24 @@ void MinAdaptGlok::setup_style()
   double **v = atom->v;
   int nlocal = atom->nlocal;
 
+  // print the parameters used within adaptglok into the log
+
+  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",
+      dmax, delaystep, dtgrow, dtshrink, alpha0, alphashrink, tmax, tmin, 
+      s1[integrator], s2[halfstepback_flag], s3[relaxbox_flag], relaxbox_mod,
+      relaxbox_rate, ptol);
+  }
+
+  // initialize the velocities
+
   for (int i = 0; i < nlocal; i++)
     v[i][0] = v[i][1] = v[i][2] = 0.0;
 }
@@ -122,7 +139,7 @@ void MinAdaptGlok::relax_box()
   // rescale simulation box and scale atom coords for all atoms
 
   double **x = atom->x;
-  double epsilon;
+  double epsilon,disp;
   double *current_pressure_v;
   int *mask = atom->mask;
   n = atom->nlocal + atom->nghost;
@@ -144,7 +161,9 @@ void MinAdaptGlok::relax_box()
   epsilon = pressure->scalar / relaxbox_mod;
   for (int i = 0; i < 3; i++) {
     if (relaxbox_flag == 2) epsilon = pressure->vector[i] / relaxbox_mod;
-    domain->boxhi[i] += p_flag[i] * domain->boxhi[i] * epsilon * relaxbox_rate;
+    disp = domain->boxhi[i] * epsilon * relaxbox_rate;
+    if (fabs(disp) > dmax) disp > 0.0 ? disp = dmax : disp = -1 * dmax;
+    domain->boxhi[i] += p_flag[i] * disp;
   }
 
   // reset global and local box to new size/shape