Description:
-Reset the timestep size every N steps during a run, based on current -atom velocities and forces. 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 a solid, causing a damage cascade. +
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. +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 +a solid, causing a damage cascade.
This fix overrides the timestep size setting made by the timestep command. The new timestep size dt is -computed in the following way. +computed in the following manner.
-Vmax is the maximum velocity; Amax is the maximum acceleration = -force/mass. Note that Tmin or Tmax can be specified as INF, in which -case one or both of the last 2 checks will not be performed. +
For each atom, the timestep is computed that would cause it to +displace Xmax on the next integration step, as a function of its +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. +
+Given this putative timestep for each atom, the minimum timestep value +across all atoms is computed. Then the Tmin and Tmax bounds are +applied, if specified. If one (or both) is specified as NULL, it is +not applied.
When the run style is respa, this fix resets the outer loop (largest) timestep, which is the same timestep that the diff -Naur lammps-5Jun09/doc/fix_dt_reset.txt lammps-6Jun09/doc/fix_dt_reset.txt --- lammps-5Jun09/doc/fix_dt_reset.txt 2008-02-29 18:13:20.000000000 -0700 +++ lammps-6Jun09/doc/fix_dt_reset.txt 2009-06-05 15:47:11.000000000 -0600 @@ -15,8 +15,8 @@ ID, group-ID are documented in "fix"_fix.html command dt/reset = style name of this fix command N = recompute dt every N timesteps -Tmin = minimum dt allowed (can be INF) (time units) -Tmax = maximum dt allowed (can be INF) (time units) +Tmin = minimum dt allowed (can be NULL) (time units) +Tmax = maximum dt allowed (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 @@ -31,27 +31,29 @@ [Description:] -Reset the timestep size every N steps during a run, based on current -atom velocities and forces. 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 a solid, causing a damage cascade. +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. +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 +a solid, causing a damage cascade. This fix overrides the timestep size setting made by the "timestep"_timestep.html command. The new timestep size {dt} is -computed in the following way. +computed in the following manner. -compute Vmax of any atom in group -compute Amax of any atom in group -dt1 = Xmax/Vmax -dt2 = sqrt(2 Xmax/Amax) -new dt = MIN(dt1,dt2) -if dt < Tmin, dt = Tmin -if dt > Tmax, dt = Tmax :ul - -Vmax is the maximum velocity; Amax is the maximum acceleration = -force/mass. Note that Tmin or Tmax can be specified as INF, in which -case one or both of the last 2 checks will not be performed. +For each atom, the timestep is computed that would cause it to +displace {Xmax} on the next integration step, as a function of its +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. + +Given this putative timestep for each atom, the minimum timestep value +across all atoms is computed. Then the {Tmin} and {Tmax} bounds are +applied, if specified. If one (or both) is specified as NULL, it is +not applied. When the "run style"_run_style.html is {respa}, this fix resets the outer loop (largest) timestep, which is the same timestep that the diff -Naur lammps-5Jun09/src/MANYBODY/pair_eam.cpp lammps-6Jun09/src/MANYBODY/pair_eam.cpp --- lammps-5Jun09/src/MANYBODY/pair_eam.cpp 2009-04-01 16:23:28.000000000 -0600 +++ lammps-6Jun09/src/MANYBODY/pair_eam.cpp 2009-06-05 13:14:05.000000000 -0600 @@ -488,6 +488,7 @@ { int i,j,k,m,n; int ntypes = atom->ntypes; + double sixth = 1.0/6.0; // determine max function params from all active funcfl files // active means some element is pointing at it via map @@ -540,10 +541,10 @@ k = MAX(k,2); p -= k; p = MIN(p,2.0); - cof1 = -0.166666667*p*(p-1.0)*(p-2.0); + cof1 = -sixth*p*(p-1.0)*(p-2.0); cof2 = 0.5*(p*p-1.0)*(p-2.0); cof3 = -0.5*p*(p+1.0)*(p-2.0); - cof4 = 0.166666667*p*(p*p-1.0); + cof4 = sixth*p*(p*p-1.0); frho[n][m] = cof1*file->frho[k-1] + cof2*file->frho[k] + cof3*file->frho[k+1] + cof4*file->frho[k+2]; } @@ -587,10 +588,10 @@ k = MAX(k,2); p -= k; p = MIN(p,2.0); - cof1 = -0.166666667*p*(p-1.0)*(p-2.0); + cof1 = -sixth*p*(p-1.0)*(p-2.0); cof2 = 0.5*(p*p-1.0)*(p-2.0); cof3 = -0.5*p*(p+1.0)*(p-2.0); - cof4 = 0.166666667*p*(p*p-1.0); + cof4 = sixth*p*(p*p-1.0); rhor[n][m] = cof1*file->rhor[k-1] + cof2*file->rhor[k] + cof3*file->rhor[k+1] + cof4*file->rhor[k+2]; } @@ -636,10 +637,10 @@ k = MAX(k,2); p -= k; p = MIN(p,2.0); - cof1 = -0.166666667*p*(p-1.0)*(p-2.0); + cof1 = -sixth*p*(p-1.0)*(p-2.0); cof2 = 0.5*(p*p-1.0)*(p-2.0); cof3 = -0.5*p*(p+1.0)*(p-2.0); - cof4 = 0.166666667*p*(p*p-1.0); + cof4 = sixth*p*(p*p-1.0); zri = cof1*ifile->zr[k-1] + cof2*ifile->zr[k] + cof3*ifile->zr[k+1] + cof4*ifile->zr[k+2]; @@ -649,10 +650,10 @@ k = MAX(k,2); p -= k; p = MIN(p,2.0); - cof1 = -0.166666667*p*(p-1.0)*(p-2.0); + cof1 = -sixth*p*(p-1.0)*(p-2.0); cof2 = 0.5*(p*p-1.0)*(p-2.0); cof3 = -0.5*p*(p+1.0)*(p-2.0); - cof4 = 0.166666667*p*(p*p-1.0); + cof4 = sixth*p*(p*p-1.0); zrj = cof1*jfile->zr[k-1] + cof2*jfile->zr[k] + cof3*jfile->zr[k+1] + cof4*jfile->zr[k+2]; diff -Naur lammps-5Jun09/src/fix_dt_reset.cpp lammps-6Jun09/src/fix_dt_reset.cpp --- lammps-5Jun09/src/fix_dt_reset.cpp 2009-02-12 14:36:52.000000000 -0700 +++ lammps-6Jun09/src/fix_dt_reset.cpp 2009-06-05 15:45:48.000000000 -0600 @@ -56,15 +56,15 @@ minbound = maxbound = 1; tmin = tmax = 0.0; - if (strcmp(arg[4],"INF") == 0) minbound = 0; + if (strcmp(arg[4],"NULL") == 0) minbound = 0; else tmin = atof(arg[4]); - if (strcmp(arg[5],"INF") == 0) maxbound = 0; + if (strcmp(arg[5],"NULL") == 0) maxbound = 0; else tmax = atof(arg[5]); xmax = atof(arg[6]); if (minbound && tmin < 0.0) error->all("Illegal fix dt/reset command"); if (maxbound && tmax < 0.0) error->all("Illegal fix dt/reset command"); - if (minbound && maxbound && tmin > tmax) + if (minbound && maxbound && tmin >= tmax) error->all("Illegal fix dt/reset command"); if (xmax <= 0.0) error->all("Illegal fix dt/reset command"); @@ -117,6 +117,8 @@ if ((strcmp(output->dump[i]->style,"dcd") == 0 || strcmp(output->dump[i]->style,"xtc") == 0) && comm->me == 0) error->warning("Dump dcd/xtc timestamp may be wrong with fix dt/reset"); + + ftm2v = force->ftm2v; } /* ---------------------------------------------------------------------- */ @@ -130,6 +132,10 @@ void FixDtReset::end_of_step() { + double dt,dtv,dtf,dtsq; + double vsq,fsq,massinv; + double delx,dely,delz,delr; + // accumulate total time based on previous timestep t_elapsed += (update->ntimestep - laststep) * update->dt; @@ -144,33 +150,29 @@ int *mask = atom->mask; int nlocal = atom->nlocal; - double vsq,fsq,asq,ms; - double bound[2],bound_all[2]; - bound[0] = bound[1] = 0.0; + double dtmin = BIG; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { + if (rmass) massinv = 1.0/rmass[i]; + 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]; - if (rmass) ms = rmass[i]; - else ms = mass[type[i]]; - asq = fsq/ms/ms; - bound[0] = MAX(bound[0],vsq); - bound[1] = MAX(bound[1],asq); + dtv = dtf = 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); + 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; + delz = dt*v[i][2] + 0.5*dtsq*massinv*f[i][2] * ftm2v; + delr = sqrt(delx*delx + dely*dely + delz*delz); + if (delr > xmax) dt *= xmax/delr; + dtmin = MIN(dtmin,dt); } - MPI_Allreduce(bound,bound_all,2,MPI_DOUBLE,MPI_MAX,world); - - // set new timestep - - double dt,dt1,dt2; - - if (bound_all[0] > 0.0) dt1 = xmax/sqrt(bound_all[0]); - else dt1 = BIG; - if (bound_all[1] > 0.0) dt2 = sqrt(2.0 * xmax/sqrt(bound_all[1])); - else dt2 = BIG; + MPI_Allreduce(&dtmin,&dt,1,MPI_DOUBLE,MPI_MIN,world); - dt = MIN(dt1,dt2); if (minbound) dt = MAX(dt,tmin); if (maxbound) dt = MIN(dt,tmax); diff -Naur lammps-5Jun09/src/fix_dt_reset.h lammps-6Jun09/src/fix_dt_reset.h --- lammps-5Jun09/src/fix_dt_reset.h 2008-01-09 14:56:57.000000000 -0700 +++ lammps-6Jun09/src/fix_dt_reset.h 2009-06-05 15:45:48.000000000 -0600 @@ -32,6 +32,7 @@ private: int minbound,maxbound,laststep; double tmin,tmax,xmax; + double ftm2v; double t_elapsed; int respaflag; };