diff --git a/src/min.cpp b/src/min.cpp index bf2149cd19..811076f37a 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -65,6 +65,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp) integrator = 0; halfstepback_flag = 1; delaystep_start_flag = 1; + max_vdotf_negatif = 2000; relaxbox_mod = 1000000; relaxbox_rate = 0.33; relaxbox_flag = 0; @@ -190,6 +191,10 @@ void Min::init() neighbor->dist_check = 1; niter = neval = 0; + + // store timestep size (important for variable timestep minimizer) + + dtinit = update->dt; } /* ---------------------------------------------------------------------- @@ -477,6 +482,10 @@ void Min::cleanup() modify->delete_fix("MINIMIZE"); domain->box_too_small_check(); + + // reset timestep size (important for variable timestep minimizer) + + update->dt = dtinit; } /* ---------------------------------------------------------------------- @@ -699,6 +708,10 @@ void Min::modify_params(int narg, char **arg) else if (strcmp(arg[iarg+1],"no") == 0) delaystep_start_flag = 0; else error->all(FLERR,"Illegal min_modify command"); iarg += 2; + } else if (strcmp(arg[iarg],"vdfmax") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); + max_vdotf_negatif = 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; @@ -924,6 +937,7 @@ char *Min::stopstrings(int n) "quadratic factors are zero", "trust region too small", "HFTN minimizer error", - "walltime limit reached"}; + "walltime limit reached", + "max iterations with v.f negative"}; return (char *) strings[n]; } diff --git a/src/min.h b/src/min.h index 6cb7fcfcd6..aaaba90957 100644 --- a/src/min.h +++ b/src/min.h @@ -48,7 +48,7 @@ class Min : protected Pointers { // possible return values of iterate() method enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE, - ZEROQUAD,TRSMALL,INTERROR,TIMEOUT}; + ZEROQUAD,TRSMALL,INTERROR,TIMEOUT,MAXVDOTF}; protected: int eflag,vflag; // flags for energy/virial computation @@ -57,6 +57,7 @@ class Min : protected Pointers { double dmax; // max dist to move any atom in one step int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero + double dtinit; // store the default timestep // only for minimize style adaptglok int delaystep; // minium steps of dynamics @@ -66,6 +67,7 @@ class Min : protected Pointers { int integrator; // Newton integration: euler, leapfrog, verlet... 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 diff --git a/src/min_adaptglok.cpp b/src/min_adaptglok.cpp index 8deaa82691..9173be9b85 100644 --- a/src/min_adaptglok.cpp +++ b/src/min_adaptglok.cpp @@ -59,11 +59,12 @@ void MinAdaptGlok::init() if (p_flag[2] && domain->zperiodic == 0) error->all(FLERR,"Cannot use boxrelax on a non-periodic dimension"); - dt = dtinit = update->dt; + dt = update->dt; dtmax = tmax * dt; dtmin = tmin * dt; alpha = alpha0; last_negative = ntimestep_start = update->ntimestep; + vdotf_negatif = 0; if (relaxbox_flag){ int icompute = modify->find_compute("thermo_press"); @@ -269,6 +270,7 @@ int MinAdaptGlok::iterate(int maxiter) if (vdotfall > 0.0) { vdotv = 0.0; + vdotf_negatif = 0; for (int i = 0; i < nlocal; i++) vdotv += v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world); @@ -325,6 +327,15 @@ int MinAdaptGlok::iterate(int maxiter) update->dt = dt; } } + + // stopping criterion while stuck in a local bassin + + vdotf_negatif++; + if (max_vdotf_negatif > 0 && vdotf_negatif > max_vdotf_negatif) + return MAXVDOTF; + + // inertia correction + if (halfstepback_flag) { for (int i = 0; i < nlocal; i++) { x[i][0] -= 0.5 * dtv * v[i][0]; @@ -529,10 +540,8 @@ int MinAdaptGlok::iterate(int maxiter) if (update->multireplica == 0) { if (fabs(ecurrent-eprevious) < update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY) - && pflag) { - update->dt = dtinit; + && pflag) return ETOL; - } } else { if (fabs(ecurrent-eprevious) < update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY) @@ -540,10 +549,8 @@ int MinAdaptGlok::iterate(int maxiter) flag = 0; else flag = 1; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld); - if (flagall == 0) { - update->dt = dtinit; + if (flagall == 0) return ETOL; - } } } @@ -554,18 +561,14 @@ int MinAdaptGlok::iterate(int maxiter) if (update->ftol > 0.0) { fdotf = fnorm_sqr(); if (update->multireplica == 0) { - if (fdotf < update->ftol*update->ftol && pflag) { - update->dt = dtinit; + if (fdotf < update->ftol*update->ftol && pflag) return FTOL; - } } else { if (fdotf < update->ftol*update->ftol && pflag) flag = 0; else flag = 1; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld); - if (flagall == 0) { - update->dt = dtinit; + if (flagall == 0) return FTOL; - } } } @@ -578,8 +581,5 @@ 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 7cc614bbfd..75d1026784 100644 --- a/src/min_adaptglok.h +++ b/src/min_adaptglok.h @@ -36,10 +36,10 @@ class MinAdaptGlok : public Min { int iterate(int); private: - double dt,dtmax,dtmin,dtinit; + double dt,dtmax,dtmin; double alpha; bigint last_negative,ntimestep_start; - int pflag; + int pflag,vdotf_negatif; class Compute *temperature,*pressure; double boxlo[3],h_inv[6]; };