diff --git a/doc/src/min_modify.txt b/doc/src/min_modify.txt index 882ef763bf..c9ad9fce35 100644 --- a/doc/src/min_modify.txt +++ b/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} - or {alpha0} or {alpha_shrink} or {tmax} or {tmin} or {integrator} +keyword = {dmax} or {line} + 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 - steps = number of steps for dynamics temporization, adaptglok only - {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 + {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} or {leapfrog} + eulerimplicit,eulerexplicit,verlet,leapfrog = integration scheme (adaptglok) {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 - {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} - or {leapfrog} - eulerimplicit,eulerexplicit,verlet,leapfrog = integration scheme, adaptglok - {halfstepback} value = {yes} or {no}, adaptglok only :pre + coef = factor defining the minimum timestep (adaptglok) + {delaystep} value = steps + steps = number of steps for dynamics temporization (adaptglok) + {dtgrow} value = grow + 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}, -{alpha0}, {alphashrink}, {tmax} and {tmin}. +to optimize the relaxation: {integrator}, {tmax}, {tmin}, {delaystep}, {dtgrow}, {dtshrink}, +{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, -{tmax} = 2.0 means that the maximum value the timestep can reach is twice the -value defined by "timestep"_timestep.html. +multiplication factor applied to the "timestep"_timestep.html. Thus +{tmax} = 2.0 in metal "units"_units.html means that the maximum value +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, -alpha0 = 0.25, alpha_shrink = 0.99, tmax = 2.0, tmin = 0.02, adaptstep = yes +integrator = eulerimplicit, tmax = 2.0, tmin = 0.02, +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. diff --git a/doc/src/min_style.txt b/doc/src/min_style.txt index 9f33bac948..9cdc29ddb7 100644 --- a/doc/src/min_style.txt +++ b/doc/src/min_style.txt @@ -77,7 +77,9 @@ 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. +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 diff --git a/src/min.cpp b/src/min.cpp index e24247aae3..8e775e5609 100644 --- a/src/min.cpp +++ b/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; diff --git a/src/min_adaptglok.cpp b/src/min_adaptglok.cpp index d9c488320f..3c54201e09 100644 --- a/src/min_adaptglok.cpp +++ b/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) - error->all(FLERR,"relaxbox_mod has to be positif"); + if (relaxbox_mod < 0.0) error->all(FLERR,"relaxbox_mod has to be positif"); // require periodicity in boxrelax dimensions @@ -64,12 +65,10 @@ void MinAdaptGlok::init() alpha = alpha0; 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]; -} + if (relaxbox_flag){ + int 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