Merge pull request #933 from rtoijala/fix_dt_reset_energy

Limit atom energy change in fix dt/reset

doc/src/fix_dt_reset.txt

@@ -19,7 +19,9 @@ Tmin = minimum dt allowed which can be NULL (time units)
 Tmax = maximum dt allowed which can be NULL (time units)
 Xmax = maximum distance for an atom to move in one timestep (distance units)
 zero or more keyword/value pairs may be appended
-keyword = {units} :ul
+keyword = {emax} or {units} :ul
+  {emax} value = Emax
+    Emax = maximum kinetic energy change for an atom in one timestep (energy units)
   {units} value = {lattice} or {box}
     lattice = Xmax is defined in lattice units
     box = Xmax is defined in simulation box units :pre
@@ -27,12 +29,14 @@ keyword = {units} :ul
 [Examples:]
 
 fix 5 all dt/reset 10 1.0e-5 0.01 0.1
-fix 5 all dt/reset 10 0.01 2.0 0.2 units box :pre
+fix 5 all dt/reset 10 0.01 2.0 0.2 units box
+fix 5 all dt/reset 5 NULL 0.001 0.5 emax 30 units box :pre
 
 [Description:]
 
 Reset the timestep size every N steps during a run, so that no atom
-moves further than Xmax, based on current atom velocities and forces.
+moves further than Xmax, based on current atom velocities and forces,
+and (optionally) no atom's kinetic energy changes by more than Emax.
 This can be useful when starting from a configuration with overlapping
 atoms, where forces will be large.  Or it can be useful when running
 an impact simulation where one or more high-energy atoms collide with
@@ -48,7 +52,9 @@ current velocity and force.  Since performing this calculation exactly
 would require the solution to a quartic equation, a cheaper estimate
 is generated.  The estimate is conservative in that the atom's
 displacement is guaranteed not to exceed {Xmax}, though it may be
-smaller.
+smaller. Also, if the {Emax} value is given, for each atom, the
+timestep is limited to a value that allows the atom's kinetic energy
+to change by at most {Emax}.
 
 Given this putative timestep for each atom, the minimum timestep value
 across all atoms is computed.  Then the {Tmin} and {Tmax} bounds are
@@ -87,4 +93,4 @@ minimization"_minimize.html.
 
 [Default:]
 
-The option defaults is units = lattice.
+The option defaults is units = lattice, and no kinetic energy change limit.

src/fix_dt_reset.cpp

@@ -69,6 +69,7 @@ FixDtReset::FixDtReset(LAMMPS *lmp, int narg, char **arg) :
 
   int scaleflag = 1;
 
+  emax = -1.0;
   int iarg = 7;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"units") == 0) {
@@ -77,6 +78,11 @@ FixDtReset::FixDtReset(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
       else error->all(FLERR,"Illegal fix dt/reset command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"emax") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix dt/reset command");
+      emax = force->numeric(FLERR,arg[iarg+1]);
+      if (emax <= 0.0) error->all(FLERR,"Illegal fix dt/reset command");
+      iarg += 2;
     } else error->all(FLERR,"Illegal fix dt/reset command");
   }
 
@@ -117,6 +123,7 @@ void FixDtReset::init()
                      "Dump dcd/xtc timestamp may be wrong with fix dt/reset");
 
   ftm2v = force->ftm2v;
+  mvv2e = force->mvv2e;
   dt = update->dt;
 }
 
@@ -131,12 +138,10 @@ void FixDtReset::setup(int vflag)
 
 void FixDtReset::end_of_step()
 {
-  double dtv,dtf,dtsq;
+  double dtv,dtf,dte,dtsq;
   double vsq,fsq,massinv;
   double delx,dely,delz,delr;
 
-  // compute vmax and amax of any atom in group
-
   double **v = atom->v;
   double **f = atom->f;
   double *mass = atom->mass;
@@ -153,10 +158,14 @@ void FixDtReset::end_of_step()
       else massinv = 1.0/mass[type[i]];
       vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
       fsq = f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
-      dtv = dtf = BIG;
+      dtv = dtf = dte = BIG;
       if (vsq > 0.0) dtv = xmax/sqrt(vsq);
       if (fsq > 0.0) dtf = sqrt(2.0*xmax/(ftm2v*sqrt(fsq)*massinv));
       dt = MIN(dtv,dtf);
+      if (emax > 0.0 && vsq > 0.0 && fsq > 0.0) {
+        dte = emax/sqrt(fsq*vsq)/sqrt(ftm2v*mvv2e);
+        dt = MIN(dt, dte);
+      }
       dtsq = dt*dt;
       delx = dt*v[i][0] + 0.5*dtsq*massinv*f[i][0] * ftm2v;
       dely = dt*v[i][1] + 0.5*dtsq*massinv*f[i][1] * ftm2v;

src/fix_dt_reset.h

@@ -37,8 +37,8 @@ class FixDtReset : public Fix {
  private:
   bigint laststep;
   int minbound,maxbound;
-  double tmin,tmax,xmax;
-  double ftm2v;
+  double tmin,tmax,xmax,emax;
+  double ftm2v,mvv2e;
   double dt,t_laststep;
   int respaflag;
 };