From c2e400910652defcff31ebc4559ae053db58d97a Mon Sep 17 00:00:00 2001 From: jguenole Date: Wed, 29 May 2019 10:37:53 +0200 Subject: [PATCH] Remove option relaxbox from adaptglok: wrong behavior with non-P boundaries. Code cleanup. --- doc/src/min_modify.txt | 34 +-------- doc/src/min_style.txt | 4 +- doc/src/minimize.txt | 3 +- src/min.cpp | 42 ----------- src/min.h | 6 -- src/min_adaptglok.cpp | 155 +++-------------------------------------- src/min_adaptglok.h | 6 +- 7 files changed, 13 insertions(+), 237 deletions(-) diff --git a/doc/src/min_modify.txt b/doc/src/min_modify.txt index 5a90066738..5e5b71501f 100644 --- a/doc/src/min_modify.txt +++ b/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 {ptol} or {halfstepback} or {initialdelay} or {vdfmax} + or {alpha0} or {alphashrink} 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. diff --git a/doc/src/min_style.txt b/doc/src/min_style.txt index 9e26493d3a..3a6eec25dc 100644 --- a/doc/src/min_style.txt +++ b/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 diff --git a/doc/src/minimize.txt b/doc/src/minimize.txt index 82a3587086..6f7a3cfd2c 100644 --- a/doc/src/minimize.txt +++ b/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} can also use a leapfrog or velocity Verlet integration step. +{adaptglok} perform damped dynamics using an Euler integration step. Thus they require a "timestep"_timestep.html be defined. NOTE: The damped dynamic minimizers use whatever timestep you have diff --git a/src/min.cpp b/src/min.cpp index b309e119a8..b665d74b53 100644 --- a/src/min.cpp +++ b/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"); diff --git a/src/min.h b/src/min.h index cc1a59a03e..97dfff5d20 100644 --- a/src/min.h +++ b/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; diff --git a/src/min_adaptglok.cpp b/src/min_adaptglok.cpp index 173ca6c38b..72eddc8ecc 100644 --- a/src/min_adaptglok.cpp +++ b/src/min_adaptglok.cpp @@ -11,7 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include +#include #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, - relaxbox_rate, ptol); + s1[integrator], s2[halfstepback_flag]); } // 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) - && pflag) + update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)) return ETOL; } else { if (fabs(ecurrent-eprevious) < - update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY) - && pflag) + update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)) 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) diff --git a/src/min_adaptglok.h b/src/min_adaptglok.h index b3d8968578..5159805c00 100644 --- a/src/min_adaptglok.h +++ b/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]; }; }