1 /* ----------------------------------------------------------------------
2    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3    https://www.lammps.org/, Sandia National Laboratories
4    Steve Plimpton, sjplimp@sandia.gov
5 
6    Copyright (2003) Sandia Corporation.  Under the terms of Contract
7    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8    certain rights in this software.  This software is distributed under
9    the GNU General Public License.
10 
11    See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13 
14 #include "fix_dt_reset.h"
15 
16 #include "atom.h"
17 #include "comm.h"
18 #include "domain.h"
19 #include "dump.h"
20 #include "error.h"
21 #include "fix.h"
22 #include "force.h"
23 #include "integrate.h"
24 #include "lattice.h"
25 #include "modify.h"
26 #include "output.h"
27 #include "pair.h"
28 #include "update.h"
29 
30 #include <cmath>
31 #include <cstring>
32 
33 using namespace LAMMPS_NS;
34 using namespace FixConst;
35 
36 #define BIG 1.0e20
37 
38 /* ---------------------------------------------------------------------- */
39 
FixDtReset(LAMMPS * lmp,int narg,char ** arg)40 FixDtReset::FixDtReset(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
41 {
42   if (narg < 7) error->all(FLERR, "Illegal fix dt/reset command");
43 
44   // set time_depend, else elapsed time accumulation can be messed up
45 
46   time_depend = 1;
47   scalar_flag = 1;
48   global_freq = 1;
49   extscalar = 0;
50   extvector = 0;
51   dynamic_group_allow = 1;
52 
53   nevery = utils::inumeric(FLERR, arg[3], false, lmp);
54   if (nevery <= 0) error->all(FLERR, "Illegal fix dt/reset command");
55 
56   minbound = maxbound = 1;
57   tmin = tmax = 0.0;
58   if (strcmp(arg[4], "NULL") == 0)
59     minbound = 0;
60   else
61     tmin = utils::numeric(FLERR, arg[4], false, lmp);
62   if (strcmp(arg[5], "NULL") == 0)
63     maxbound = 0;
64   else
65     tmax = utils::numeric(FLERR, arg[5], false, lmp);
66   xmax = utils::numeric(FLERR, arg[6], false, lmp);
67 
68   if (minbound && tmin < 0.0) error->all(FLERR, "Illegal fix dt/reset command");
69   if (maxbound && tmax < 0.0) error->all(FLERR, "Illegal fix dt/reset command");
70   if (minbound && maxbound && tmin >= tmax) error->all(FLERR, "Illegal fix dt/reset command");
71   if (xmax <= 0.0) error->all(FLERR, "Illegal fix dt/reset command");
72 
73   int scaleflag = 1;
74 
75   emax = -1.0;
76   int iarg = 7;
77   while (iarg < narg) {
78     if (strcmp(arg[iarg], "units") == 0) {
79       if (iarg + 2 > narg) error->all(FLERR, "Illegal fix dt/reset command");
80       if (strcmp(arg[iarg + 1], "box") == 0)
81         scaleflag = 0;
82       else if (strcmp(arg[iarg + 1], "lattice") == 0)
83         scaleflag = 1;
84       else
85         error->all(FLERR, "Illegal fix dt/reset command");
86       iarg += 2;
87     } else if (strcmp(arg[iarg], "emax") == 0) {
88       if (iarg + 2 > narg) error->all(FLERR, "Illegal fix dt/reset command");
89       emax = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
90       if (emax <= 0.0) error->all(FLERR, "Illegal fix dt/reset command");
91       iarg += 2;
92     } else
93       error->all(FLERR, "Illegal fix dt/reset command");
94   }
95 
96   // setup scaling, based on xlattice parameter
97 
98   if (scaleflag) xmax *= domain->lattice->xlattice;
99 
100   // initializations
101 
102   t_laststep = 0.0;
103   laststep = update->ntimestep;
104 }
105 
106 /* ---------------------------------------------------------------------- */
107 
setmask()108 int FixDtReset::setmask()
109 {
110   int mask = 0;
111   mask |= END_OF_STEP;
112   return mask;
113 }
114 
115 /* ---------------------------------------------------------------------- */
116 
init()117 void FixDtReset::init()
118 {
119   // set rRESPA flag
120 
121   respaflag = 0;
122   if (utils::strmatch(update->integrate_style, "^respa")) respaflag = 1;
123 
124   // check for DCD or XTC dumps
125 
126   for (int i = 0; i < output->ndump; i++)
127     if ((strcmp(output->dump[i]->style, "dcd") == 0 ||
128          strcmp(output->dump[i]->style, "xtc") == 0) &&
129         comm->me == 0)
130       error->warning(FLERR, "Dump dcd/xtc timestamp may be wrong with fix dt/reset");
131 
132   ftm2v = force->ftm2v;
133   mvv2e = force->mvv2e;
134   dt = update->dt;
135 }
136 
137 /* ---------------------------------------------------------------------- */
138 
setup(int)139 void FixDtReset::setup(int /*vflag*/)
140 {
141   end_of_step();
142 }
143 
144 /* ---------------------------------------------------------------------- */
145 
end_of_step()146 void FixDtReset::end_of_step()
147 {
148   double dtv, dtf, dte, dtsq;
149   double vsq, fsq, massinv;
150   double delx, dely, delz, delr;
151 
152   double **v = atom->v;
153   double **f = atom->f;
154   double *mass = atom->mass;
155   double *rmass = atom->rmass;
156   int *type = atom->type;
157   int *mask = atom->mask;
158   int nlocal = atom->nlocal;
159 
160   double dtmin = BIG;
161 
162   for (int i = 0; i < nlocal; i++)
163     if (mask[i] & groupbit) {
164       if (rmass)
165         massinv = 1.0 / rmass[i];
166       else
167         massinv = 1.0 / mass[type[i]];
168       vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
169       fsq = f[i][0] * f[i][0] + f[i][1] * f[i][1] + f[i][2] * f[i][2];
170       dtv = dtf = dte = BIG;
171       if (vsq > 0.0) dtv = xmax / sqrt(vsq);
172       if (fsq > 0.0) dtf = sqrt(2.0 * xmax / (ftm2v * sqrt(fsq) * massinv));
173       dt = MIN(dtv, dtf);
174       if ((emax > 0.0) && (fsq * vsq > 0.0)) {
175         dte = emax / sqrt(fsq * vsq) / sqrt(ftm2v * mvv2e);
176         dt = MIN(dt, dte);
177       }
178       dtsq = dt * dt;
179       delx = dt * v[i][0] + 0.5 * dtsq * massinv * f[i][0] * ftm2v;
180       dely = dt * v[i][1] + 0.5 * dtsq * massinv * f[i][1] * ftm2v;
181       delz = dt * v[i][2] + 0.5 * dtsq * massinv * f[i][2] * ftm2v;
182       delr = sqrt(delx * delx + dely * dely + delz * delz);
183       if (delr > xmax) dt *= xmax / delr;
184       dtmin = MIN(dtmin, dt);
185     }
186 
187   MPI_Allreduce(&dtmin, &dt, 1, MPI_DOUBLE, MPI_MIN, world);
188 
189   if (minbound) dt = MAX(dt, tmin);
190   if (maxbound) dt = MIN(dt, tmax);
191 
192   // if timestep didn't change, just return
193   // else reset update->dt and other classes that depend on it
194   // rRESPA, pair style, fixes
195 
196   if (dt == update->dt) return;
197 
198   laststep = update->ntimestep;
199 
200   update->update_time();
201   update->dt = dt;
202   update->dt_default = 0;
203   if (respaflag) update->integrate->reset_dt();
204   if (force->pair) force->pair->reset_dt();
205   for (int i = 0; i < modify->nfix; i++) modify->fix[i]->reset_dt();
206 }
207 
208 /* ---------------------------------------------------------------------- */
209 
compute_scalar()210 double FixDtReset::compute_scalar()
211 {
212   return (double) laststep;
213 }
214