diff --git a/doc/src/min_modify.txt b/doc/src/min_modify.txt index a97bf8eba5..882ef763bf 100644 --- a/doc/src/min_modify.txt +++ b/doc/src/min_modify.txt @@ -19,19 +19,24 @@ keyword = {dmax} or {line} or {delaystep} or {dt_grow} or {dt_shrink} 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 = {int} for adaptglok - {dt_grow} value = {float} for adaptglok - {dt_shrink} value = {float} for adaptglok - {alpha0} value = {float} for adaptglok - {alpha_shrink} value = {float} for adaptglok - {tmax} value = {float} for adaptglok - {tmin} value = {float} for adaptglok + {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 + {tmax} value = coef + coef = factor defining the maximum timestep, adaptglok only + {tmin} value = coef + coef = factor defining the minimum timestep, adaptglok only {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} - or {leapfrog} for adaptglok - {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 + or {leapfrog} + eulerimplicit,eulerexplicit,verlet,leapfrog = integration scheme, adaptglok + {halfstepback} value = {yes} or {no}, adaptglok only :pre :ule [Examples:] @@ -80,18 +85,20 @@ 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}, {dt_grow}, {dt_shrink}, -{alpha0}, {alpha_shrink}, {tmax} and {tmin}. The default values have -been chosen to be reliable in most cases, but one could be willing -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. +to optimize the relaxation: {integrator}, {delaystep}, {dtgrow}, {dtshrink}, +{alpha0}, {alphashrink}, {tmax} and {tmin}. +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. +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. [Restrictions:] none @@ -101,7 +108,7 @@ optimum in most situations. [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, halfstepback = yes, -relaxbox = no, relaxbox_modulus = 1000000 and relaxbox_rate = 0.33. +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. diff --git a/doc/src/min_style.txt b/doc/src/min_style.txt index b93e8f99e7..9f33bac948 100644 --- a/doc/src/min_style.txt +++ b/doc/src/min_style.txt @@ -66,10 +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 +Either the {quickmin}, {fire} and {adaptglok} styles are useful in the context of nudged elastic band (NEB) calculations via the "neb"_neb.html command. NOTE: The damped dynamic minimizers use whatever timestep you have @@ -77,7 +75,7 @@ defined via the "timestep"_timestep.html command. Often they will converge more quickly if you use a timestep about 10x larger than you would normally use for dynamics simulations. -NOTE: The {quickmin}, {fire}, and {hftn} styles do not yet support the +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. diff --git a/doc/src/minimize.txt b/doc/src/minimize.txt index 6441d48693..3337fd9051 100644 --- a/doc/src/minimize.txt +++ b/doc/src/minimize.txt @@ -64,9 +64,9 @@ 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 -(Velocity-Verlet for {adaptglok}). Thus they require a -"timestep"_timestep.html be defined. +{adaptglok} perform damped dynamics using an Euler integration step. The style +{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 defined via the "timestep"_timestep.html command. Often they will diff --git a/src/min.cpp b/src/min.cpp index df92fbf153..c05b763229 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -56,17 +56,14 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp) linestyle = 1; delaystep = 20; - dt_grow = 1.1; - dt_shrink = 0.5; + dtgrow = 1.1; + dtshrink = 0.5; alpha0 = 0.25; - alpha_shrink = 0.99; + alphashrink = 0.99; 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; @@ -660,21 +657,21 @@ void Min::modify_params(int narg, char **arg) if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); delaystep = force->numeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"dt_grow") == 0) { + } else if (strcmp(arg[iarg],"dtgrow") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); - dt_grow = force->numeric(FLERR,arg[iarg+1]); + dtgrow = force->numeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"dt_shrink") == 0) { + } else if (strcmp(arg[iarg],"dtshrink") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); - dt_shrink = force->numeric(FLERR,arg[iarg+1]); + dtshrink = force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"alpha0") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); alpha0 = force->numeric(FLERR,arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"alpha_shrink") == 0) { + } else if (strcmp(arg[iarg],"alphashrink") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); - alpha_shrink = force->numeric(FLERR,arg[iarg+1]); + alphashrink = force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"tmax") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); @@ -690,28 +687,13 @@ 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; else if (strcmp(arg[iarg+1],"verlet") == 0) integrator = 1; 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"); + else error->all(FLERR,"Illegal min_modify command"); 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 1a9e4be8b5..facf41dacc 100644 --- a/src/min.h +++ b/src/min.h @@ -58,15 +58,13 @@ class Min : protected Pointers { double dmax; // max dist to move any atom in one step int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero - int delaystep; // minium steps of dynamics (adaptglok) - double dt_grow,dt_shrink; // timestep increase, decrease (adaptglok) - 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 + // only for min_adaptglok + int delaystep; // minium steps of dynamics + double dtgrow,dtshrink; // timestep increase, decrease + double alpha0,alphashrink; // mixing velocities+forces coefficient + double tmax,tmin; // timestep multiplicators max, min + int integrator; // Newton integration: euler, leapfrog, verlet... + int halfstepback_flag; // half step backward when v.f <= 0.0 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 08ea84e510..f91f4fb9a8 100644 --- a/src/min_adaptglok.cpp +++ b/src/min_adaptglok.cpp @@ -20,10 +20,6 @@ #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; @@ -41,18 +37,15 @@ void MinAdaptGlok::init() { Min::init(); + if (tmax < tmin) error->all(FLERR,"tmax cannot be smaller than tmin"); + if (dtgrow < 1.0) error->all(FLERR,"dtgrow cannot be smaller than 1.0"); + if (dtshrink > 1.0) error->all(FLERR,"dtshrink cannot be greater than 1.0"); + dt = dtinit = update->dt; dtmax = tmax * dt; 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]; - } + last_negative = ntimestep_start = update->ntimestep; } /* ---------------------------------------------------------------------- */ @@ -80,68 +73,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() -{ - 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) @@ -154,13 +85,12 @@ int MinAdaptGlok::iterate(int maxiter) alpha_final = 0.0; - double **f = atom->f; - double **v = atom->v; - // Leap Frog integration initialization if (integrator == 2) { + double **f = atom->f; + double **v = atom->v; double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; @@ -170,6 +100,7 @@ int MinAdaptGlok::iterate(int maxiter) neval++; dtf = 0.5 * dt * force->ftm2v; + if (rmass) { for (int i = 0; i < nlocal; i++) { dtfm = dtf / rmass[i]; @@ -185,7 +116,6 @@ int MinAdaptGlok::iterate(int maxiter) v[i][2] = dtfm * f[i][2]; } } - } for (int iter = 0; iter < maxiter; iter++) { @@ -196,21 +126,18 @@ int MinAdaptGlok::iterate(int maxiter) ntimestep = ++update->ntimestep; niter++; - // Relax the simulation box + // pointers - if (relaxbox_flag) relax_box(); + int nlocal = atom->nlocal; + double **v = atom->v; + double **f = atom->f; + double **x = atom->x; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; // vdotfall = v dot f - // Euler || Leap Frog integration - - if (integrator == 0 || integrator == 2) { - double **v = atom->v; - double **f = atom->f; - } - int nlocal = atom->nlocal; - double **x = atom->x; - vdotf = 0.0; for (int i = 0; i < nlocal; i++) vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2]; @@ -227,9 +154,10 @@ int MinAdaptGlok::iterate(int maxiter) // if (v dot f) > 0: // v = (1-alpha) v + alpha |v| Fhat // |v| = length of v, Fhat = unit f - // The modificatin of v is made wihtin the Verlet integration, after v update + // Only: (1-alpha) and alpha |v| Fhat is calculated here + // the modificatin of v is made wihtin the integration, after v update // if more than delaystep since v dot f was negative: - // increase timestep and decrease alpha + // increase timestep, update global timestep and decrease alpha if (vdotfall > 0.0) { vdotv = 0.0; @@ -263,30 +191,29 @@ int MinAdaptGlok::iterate(int maxiter) else scale2 = alpha * sqrt(vdotvall/fdotfall); if (ntimestep - last_negative > delaystep) { - dt = MIN(dt*dt_grow,dtmax); + dt = MIN(dt*dtgrow,dtmax); update->dt = dt; - alpha *= alpha_shrink; + alpha *= alphashrink; } // else (v dot f) <= 0 // if more than delaystep since starting the relaxation: // reset alpha - // if dt > dtmin: + // if dt*dtshrink > dtmin: // decrease timestep - // set x(t) = x(t-0.5*dt) - // set v = 0 + // update global timestep (for thermo output) + // half step back within the dynamics: x(t) = x(t-0.5*dt) + // reset velocities: v = 0 } else { last_negative = ntimestep; - // Limit decrease of timestep - if (ntimestep - ntimestep_fire > delaystep) { + if (ntimestep - ntimestep_start > delaystep) { alpha = alpha0; - if (dt > dtmin) { - dt *= dt_shrink; + if (dt*dtshrink >= dtmin) { + dt *= dtshrink; update->dt = dt; } } - double **x = atom->x; if (halfstepback_flag) { for (int i = 0; i < nlocal; i++) { x[i][0] -= 0.5 * dtv * v[i][0]; @@ -300,10 +227,6 @@ int MinAdaptGlok::iterate(int maxiter) // limit timestep so no particle moves further than dmax - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - dtvone = dt; for (int i = 0; i < nlocal; i++) { @@ -321,57 +244,20 @@ int MinAdaptGlok::iterate(int maxiter) MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld); } - // Semi-implicit Euler integration + // Dynamic integration scheme: + // 0: semi-implicit Euler + // 1: velocity Verlet + // 2: leapfrog (initial half step before the iteration loop) + // 3: explicit Euler - if (integrator == 0) { + // Semi-implicit Euler OR Leap Frog integration - dtf = dtv * force->ftm2v; + if (integrator == 0 || integrator == 2) { - if (rmass) { - for (int i = 0; i < nlocal; i++) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - if (vdotfall > 0.0) { - v[i][0] = scale1*v[i][0] + scale2*f[i][0]; - v[i][1] = scale1*v[i][1] + scale2*f[i][1]; - v[i][2] = scale1*v[i][2] + scale2*f[i][2]; - } - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } else { - for (int i = 0; i < nlocal; i++) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - if (vdotfall > 0.0) { - v[i][0] = scale1*v[i][0] + scale2*f[i][0]; - v[i][1] = scale1*v[i][1] + scale2*f[i][1]; - v[i][2] = scale1*v[i][2] + scale2*f[i][2]; - } - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } + dtf = dtv * force->ftm2v; - eprevious = ecurrent; - ecurrent = energy_force(0); - neval++; - - - // Velocity Verlet integration - - }else if (integrator == 1) { - - dtf = 0.5 * dtv * force->ftm2v; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { + if (rmass) { + for (int i = 0; i < nlocal; i++) { dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; @@ -384,9 +270,9 @@ int MinAdaptGlok::iterate(int maxiter) x[i][0] += dtv * v[i][0]; x[i][1] += dtv * v[i][1]; x[i][2] += dtv * v[i][2]; - } - } else { - for (int i = 0; i < nlocal; i++) { + } + } else { + for (int i = 0; i < nlocal; i++) { dtfm = dtf / mass[type[i]]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; @@ -399,123 +285,124 @@ int MinAdaptGlok::iterate(int maxiter) x[i][0] += dtv * v[i][0]; x[i][1] += dtv * v[i][1]; x[i][2] += dtv * v[i][2]; + } } - } - - eprevious = ecurrent; - ecurrent = energy_force(0); - neval++; - if (rmass) { - for (int i = 0; i < nlocal; i++) { + eprevious = ecurrent; + ecurrent = energy_force(0); + neval++; + + // Velocity Verlet integration + + }else if (integrator == 1) { + + dtf = 0.5 * dtv * force->ftm2v; + + if (rmass) { + for (int i = 0; i < nlocal; i++) { dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; + if (vdotfall > 0.0) { + v[i][0] = scale1*v[i][0] + scale2*f[i][0]; + v[i][1] = scale1*v[i][1] + scale2*f[i][1]; + v[i][2] = scale1*v[i][2] + scale2*f[i][2]; + } + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } else { + for (int i = 0; i < nlocal; i++) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + if (vdotfall > 0.0) { + v[i][0] = scale1*v[i][0] + scale2*f[i][0]; + v[i][1] = scale1*v[i][1] + scale2*f[i][1]; + v[i][2] = scale1*v[i][2] + scale2*f[i][2]; + } + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; } - } else { - for (int i = 0; i < nlocal; i++) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; } + + eprevious = ecurrent; + ecurrent = energy_force(0); + neval++; + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + } + } else { + for (int i = 0; i < nlocal; i++) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + } + } + + // Standard Euler integration + + }else if (integrator == 3) { + + dtf = dtv * force->ftm2v; + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + dtfm = dtf / rmass[i]; + if (vdotfall > 0.0) { + v[i][0] = scale1*v[i][0] + scale2*f[i][0]; + v[i][1] = scale1*v[i][1] + scale2*f[i][1]; + v[i][2] = scale1*v[i][2] + scale2*f[i][2]; + } + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + } + } else { + for (int i = 0; i < nlocal; i++) { + dtfm = dtf / mass[type[i]]; + if (vdotfall > 0.0) { + v[i][0] = scale1*v[i][0] + scale2*f[i][0]; + v[i][1] = scale1*v[i][1] + scale2*f[i][1]; + v[i][2] = scale1*v[i][2] + scale2*f[i][2]; + } + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + } + } + + eprevious = ecurrent; + ecurrent = energy_force(0); + neval++; + } - // Leap Frog integration - - }else if (integrator == 2) { - - dtf = dtv * force->ftm2v; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - if (vdotfall > 0.0) { - v[i][0] = scale1*v[i][0] + scale2*f[i][0]; - v[i][1] = scale1*v[i][1] + scale2*f[i][1]; - v[i][2] = scale1*v[i][2] + scale2*f[i][2]; - } - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } else { - for (int i = 0; i < nlocal; i++) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - if (vdotfall > 0.0) { - v[i][0] = scale1*v[i][0] + scale2*f[i][0]; - v[i][1] = scale1*v[i][1] + scale2*f[i][1]; - v[i][2] = scale1*v[i][2] + scale2*f[i][2]; - } - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - - eprevious = ecurrent; - ecurrent = energy_force(0); - neval++; - - // Standard Euler integration - - }else if (integrator == 3) { - - dtf = dtv * force->ftm2v; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - dtfm = dtf / rmass[i]; - if (vdotfall > 0.0) { - v[i][0] = scale1*v[i][0] + scale2*f[i][0]; - v[i][1] = scale1*v[i][1] + scale2*f[i][1]; - v[i][2] = scale1*v[i][2] + scale2*f[i][2]; - } - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - } - } else { - for (int i = 0; i < nlocal; i++) { - dtfm = dtf / mass[type[i]]; - if (vdotfall > 0.0) { - v[i][0] = scale1*v[i][0] + scale2*f[i][0]; - v[i][1] = scale1*v[i][1] + scale2*f[i][1]; - v[i][2] = scale1*v[i][2] + scale2*f[i][2]; - } - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - } - } - - eprevious = ecurrent; - ecurrent = energy_force(0); - neval++; - - } - // energy tolerance criterion // only check after delaystep elapsed since velocties reset to 0 // sync across replicas if running multi-replica minimization + // reset the timestep to the initial value 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)){ update->dt = dtinit; return ETOL; } @@ -534,6 +421,7 @@ int MinAdaptGlok::iterate(int maxiter) // force tolerance criterion // sync across replicas if running multi-replica minimization + // reset the timestep to the initial value if (update->ftol > 0.0) { fdotf = fnorm_sqr(); @@ -562,6 +450,8 @@ int MinAdaptGlok::iterate(int maxiter) } } + // reset the timestep to the initial value + update->dt = dtinit; return MAXITER; } diff --git a/src/min_adaptglok.h b/src/min_adaptglok.h index 8ae8f0a947..bd67249fc8 100644 --- a/src/min_adaptglok.h +++ b/src/min_adaptglok.h @@ -31,16 +31,12 @@ 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]; + bigint last_negative,ntimestep_start; }; }