diff --git a/doc/src/fix_ffl.txt b/doc/src/fix_ffl.txt new file mode 100644 index 0000000000..8d6b243142 --- /dev/null +++ b/doc/src/fix_ffl.txt @@ -0,0 +1,129 @@ + + + +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +fix ffl command :h3 + +[Syntax:] + +fix ID id-group ffl tau Tstart Tstop seed \[flip-type\] :pre + +ID, group-ID are documented in "fix"_fix.html command :ulb,l +ffl = style name of this fix command :l +tau = thermostat parameter (positive real) :l +Tstart, Tstop = temperature ramp during the run :l +seed = random number seed to use for generating noise (positive integer) :l +one more value may be appended :l + flip-type = determines the flipping type, can be chosen between rescale - no_flip - hard - soft, if no flip type is given, rescale will be chosen by default :pre +:ule + +[Examples:] + +fix 3 boundary ffl 10 300 300 31415 +fix 1 all ffl 100 500 500 9265 soft :pre + +[Description:] + +Apply a Fast-Forward Langevin Equation (FFL) thermostat as described +in "(Hijazi)"_#Hijazi. Contrary to +"fix langevin"_fix_langevin.html, this fix performs both +thermostatting and evolution of the Hamiltonian equations of motion, so it +should not be used together with "fix nve"_fix_nve.html -- at least not +on the same atom groups. + +The time-evolution of a single particle undergoing Langevin dynamics is described +by the equations + +\begin\{equation\} \frac \{dq\}\{dt\} = \frac\{p\}\{m\}, \end\{equation\} + +\begin\{equation\} \frac \{dp\}\{dt\} = -\gamma p + W + F, \end\{equation\} + +where \(F\) is the physical force, \(\gamma\) is the friction coefficient, and \(W\) is a +Gaussian random force. + +The friction coefficient is the inverse of the thermostat parameter : \(\gamma = 1/\tau\), with \(\tau\) the thermostat parameter {tau}. +The thermostat parameter is given in the time units, \(\gamma\) is in inverse time units. + +Equilibrium sampling a temperature T is obtained by specifying the +target value as the {Tstart} and {Tstop} arguments, so that the internal +constants depending on the temperature are computed automatically. + +The random number {seed} must be a positive integer. A Marsaglia random +number generator is used. Each processor uses the input seed to +generate its own unique seed and its own stream of random numbers. +Thus the dynamics of the system will not be identical on two runs on +different numbers of processors. + +The flipping type {flip-type} can be chosen between 4 types described in +"(Hijazi)"_#Hijazi. The flipping operation occurs during the thermostatting +step and it flips the momenta of the atoms. If no_flip is chosen, no flip +will be executed and the integration will be the same as a standard +Langevin thermostat "(Bussi)"_#Bussi. The other flipping types are : rescale - hard - soft. + +[Restart, fix_modify, output, run start/stop, minimize info:] + +The instantaneous values of the extended variables are written to +"binary restart files"_restart.html. Because the state of the random +number generator is not saved in restart files, this means you cannot +do "exact" restarts with this fix, where the simulation continues on +the same as if no restart had taken place. However, in a statistical +sense, a restarted simulation should produce the same behavior. +Note however that you should use a different seed each time you +restart, otherwise the same sequence of random numbers will be used +each time, which might lead to stochastic synchronization and +subtle artefacts in the sampling. + +This fix can ramp its target temperature over multiple runs, using the +{start} and {stop} keywords of the "run"_run.html command. See the +"run"_run.html command for details of how to do this. + +The "fix_modify"_fix_modify.html {energy} option is supported by this +fix to add the energy change induced by Langevin thermostatting to the +system's potential energy as part of "thermodynamic +output"_thermo_style.html. + +This fix computes a global scalar which can be accessed by various +"output commands"_Section_howto.html#howto_15. The scalar is the +cumulative energy change due to this fix. The scalar value +calculated by this fix is "extensive". + +[Restrictions:] + +The GLE thermostat in its current implementation should not be used +with rigid bodies, SHAKE or RATTLE. It is expected that all the +thermostatted degrees of freedom are fully flexible, and the sampled +ensemble will not be correct otherwise. + +In order to perform constant-pressure simulations please use +"fix press/berendsen"_fix_press_berendsen.html, rather than +"fix npt"_fix_nh.html, to avoid duplicate integration of the +equations of motion. + +This fix is part of the USER-MISC package. It is only enabled if LAMMPS +was built with that package. See the "Making +LAMMPS"_Section_start.html#start_3 section for more info. + +[Related commands:] + +"fix nvt"_fix_nh.html, "fix temp/rescale"_fix_temp_rescale.html, "fix +viscous"_fix_viscous.html, "fix nvt"_fix_nh.html, "pair_style +dpd/tstat"_pair_dpd.html, "fix gld"_fix_gld.html, "fix gle"_fix_gle.html + +:line + +:link(Hijazi) +[(Hijazi)] M. Hijazi, D. M. Wilkins, M. Ceriotti, J. Chem. Phys. 148, 184109 (2018) +:link(Bussi) +[(Bussi)] G. Bussi, M. Parrinello, Phs. Rev. E 75, 056707 (2007) + diff --git a/src/USER-MISC/README b/src/USER-MISC/README index eb221c07db..b3c3e85948 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -44,6 +44,7 @@ dihedral_style table/cut, Mike Salerno, ksalerno@pha.jhu.edu, 11 May 18 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015 fix bond/react, Jacob Gissinger (CU Boulder), info at disarmmd.org, 24 Feb 2018 +fix ffl, David Wilkins (EPFL Lausanne), david.wilkins @ epfl.ch, 28 Sep 2018 fix filter/corotate, Lukas Fath (KIT), lukas.fath at kit.edu, 15 Mar 2017 fix flow/gauss, Joel Eaves (CU Boulder), Joel.Eaves@Colorado.edu, 23 Aug 2016 fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014 diff --git a/src/USER-MISC/fix_ffl.cpp b/src/USER-MISC/fix_ffl.cpp new file mode 100644 index 0000000000..54e5c03a62 --- /dev/null +++ b/src/USER-MISC/fix_ffl.cpp @@ -0,0 +1,527 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Fast-forward Langevin thermostat, see + M. Hijazi, D. M. Wilkins, M. Ceriotti, J. Chem. Phys. 148, 184109 (2018) +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Lionel Constantin (EPFL), David M. Wilkins (EPFL), + Michele Ceriotti (EPFL) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "fix_ffl.h" +#include "math_extra.h" +#include "atom.h" +#include "atom_vec_ellipsoid.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "domain.h" +#include "region.h" +#include "respa.h" +#include "comm.h" +#include "input.h" +#include "variable.h" +#include "random_mars.h" +#include "memory.h" +#include "error.h" +#include "group.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum {NOBIAS,BIAS}; +enum {CONSTANT,EQUAL,ATOM}; +//#define FFL_DEBUG 1 + +#define MAXLINE 1024 + +/* syntax for fix_ffl: + * fix nfix id-group ffl tau Tstart Tstop seed [flip_type] + * + * */ + +/* ---------------------------------------------------------------------- */ + + +FixFFL::FixFFL(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + + + if (narg < 7) + error->all(FLERR,"Illegal fix ffl command. Expecting: fix " + " ffl "); + + restart_peratom = 1; + time_integrate = 1; + scalar_flag = 1; + + //gamma = 1/ time constant(tau) + if (force->numeric(FLERR,arg[3]) <= 0) + error->all(FLERR,"Illegal fix ffl tau value, should be greater than 0"); + gamma = 1.0/force->numeric(FLERR,arg[3]); + ffl_every=1; + ffl_step=0; + + // start temperature (t ramp) + t_start = force->numeric(FLERR,arg[4]); + + // final temperature (t ramp) + t_stop = force->numeric(FLERR,arg[5]); + + // PRNG seed + int seed = force->inumeric(FLERR,arg[6]); + + // Flip type used, uses rescale if no flip is given + if (narg == 8) strcpy(flip_type,arg[7]); + else strcpy(flip_type,"rescale"); + + + // Tests if the flip type is valid + if ( strcmp(flip_type,"rescale") && strcmp(flip_type,"hard") + && strcmp(flip_type,"soft") && strcmp(flip_type,"no_flip") ) + error->all(FLERR,"Illegal fix ffl flip type, only accepts : rescale - hard - soft - no_flip"); + + if ( strcmp(flip_type,"no_flip") == 0 ) flip_int = 0; + if ( strcmp(flip_type,"rescale") == 0 ) flip_int = 1; + if ( strcmp(flip_type,"hard") == 0 ) flip_int = 2; + if ( strcmp(flip_type,"soft") == 0 ) flip_int = 3; + + t_target=t_start; + const double kT = t_target * force->boltz / force->mvv2e; + + + // initialize Marsaglia RNG with processor-unique seed + // NB: this means runs will not be the same with different numbers of processors + if (seed <= 0) error->all(FLERR,"Illegal fix ffl command"); + random = new RanMars(lmp,seed + comm->me); + + // allocate per-type arrays for mass-scaling + sqrt_m=NULL; + memory->grow(sqrt_m, atom->ntypes+1,"ffl:sqrt_m"); + + // allocates space for temporaries + ffl_tmp1=ffl_tmp2=NULL; + + grow_arrays(atom->nmax); + + // add callbacks to enable restarts + atom->add_callback(0); + atom->add_callback(1); + + energy = 0.0; +} + + +/* --- Frees up memory used by temporaries and buffers ------------------ */ + +FixFFL::~FixFFL() +{ + delete random; + + memory->destroy(sqrt_m); + memory->destroy(ffl_tmp1); + memory->destroy(ffl_tmp2); +} + +/* ---------------------------------------------------------------------- */ + +int FixFFL::setmask() +{ + int mask = 0; + + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= INITIAL_INTEGRATE_RESPA; + mask |= FINAL_INTEGRATE_RESPA; + mask |= THERMO_ENERGY; + + + return mask; +} + +/* ------- Initializes one-time quantities for FFL ---------------------- */ + +void FixFFL::init() +{ + doffl = 1; + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + + // set force prefactors + if (!atom->rmass) + { + for (int i = 1; i <= atom->ntypes; i++) + { + sqrt_m[i] = sqrt(atom->mass[i]); + } + } + + if (strstr(update->integrate_style,"respa")) + { + nlevels_respa = ((Respa *) update->integrate)->nlevels; + step_respa = ((Respa *) update->integrate)->step; + } + + init_ffl(); +} + +/* ------- Initializes constants for FFL (change with T and dt) ------- */ + +void FixFFL::init_ffl() +{ + const double kT = t_target * force->boltz / force->mvv2e; + + // compute constants for FFL + + c1 = exp ( - gamma * 0.5 * dtv ); + c2 = sqrt( (1.0 - c1*c1)* kT ); //without the mass term + + +} + + + +/* ---------------------------------------------------------------------- */ + +void FixFFL::setup(int vflag) +{ + if (strstr(update->integrate_style,"verlet")) + post_force(vflag); + else + { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } +} + +void FixFFL::ffl_integrate() +{ + double **v = atom->v; + double *rmass = atom->rmass, smi, ismi; + double factor; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // loads momentum data (mass-scaled) into the temporary vectors for the propagation + int nk=0; + double deltae=0.0; + for (int i = 0; i < nlocal; i++) + { + if (mask[i] & groupbit) + { + if (rmass) smi = sqrt(rmass[i]); + else smi = sqrt_m[type[i]]; + + for (int k = 0; k<3; k++) + { + // first loads velocities and accumulates conserved quantity + ffl_tmp2[nk] = v[i][k] * smi; + deltae += ffl_tmp2[nk] * ffl_tmp2[nk]; + nk++; + } + } + } + + //fills up a vector of random numbers + for (int i = 0; i < nk; i++) ffl_tmp1[i] = random->gaussian(); + + + // unloads momentum data (mass-scaled) from the temporary vectors + nk=0; + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) + { + if (rmass) ismi = 1.0 / sqrt(rmass[i]); + else ismi = 1.0/ sqrt_m[type[i]]; + + for (int k = 0; k<3; k++) + { + // fetches new velocities and completes computation of the conserved quantity change + v[i][k]= c1*v[i][k] + c2*ffl_tmp1[nk]*ismi; + + deltae-= v[i][k]*v[i][k] /ismi /ismi; + + //flips the sign of the momentum (HARD FLIP) + if ( flip_int == 2) + { + if (v[i][k]*ffl_tmp2[nk] < 0.0) v[i][k] = -v[i][k]; + } + + nk++; + } + } + + //rescale operation (RESCALE FLIP) + if (flip_int == 1) + { + nk=0; + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) + { + factor = sqrt ((v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) / + (ffl_tmp2[nk]*ffl_tmp2[nk] + ffl_tmp2[nk+1]*ffl_tmp2[nk+1] + + ffl_tmp2[nk+2]*ffl_tmp2[nk+2])); + + for (int k = 0; k<3; k++) + { + v[i][k]= factor * ffl_tmp2[nk]; + nk++; + } + } + } + + + //soft flip operation (SOFT FLIP) + if (flip_int == 3) + { + nk=0; + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) + { + factor = v[i][0]*ffl_tmp2[nk] + v[i][1]*ffl_tmp2[nk+1] + v[i][2]*ffl_tmp2[nk+2]; + if (factor < 0) + { + factor = factor / (ffl_tmp2[nk]*ffl_tmp2[nk] + ffl_tmp2[nk+1]*ffl_tmp2[nk+1] + + ffl_tmp2[nk+2]*ffl_tmp2[nk+2]); + + for (int k = 0; k<3; k++) + { + v[i][k] -= 2.0 * factor * ffl_tmp2[nk]; + nk++; + } + } + else + nk += 3; + } + } + + energy += deltae*0.5*force->mvv2e; + +} + +void FixFFL::initial_integrate(int vflag) +{ + double dtfm; + + // update v and x of atoms in group + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + ffl_step--; + if (doffl && ffl_step<1) ffl_integrate(); + + if (rmass) + { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + { + 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]; + 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++) + if (mask[i] & groupbit) + { + 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]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } +} + +void FixFFL::final_integrate() +{ + double dtfm; + + // update v of atoms in group + + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + if (rmass) + { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + { + 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++) + if (mask[i] & groupbit) + { + 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 (doffl && ffl_step<1) + { + ffl_integrate(); + ffl_step = ffl_every; + } + + // Change the temperature for the next step + double delta = update->ntimestep - update->beginstep; + delta /= update->endstep - update->beginstep; + t_target = t_start + delta * (t_stop - t_start); + if (t_stop != t_start) + { + // only updates if it is really necessary + init_ffl(); + } + +} +/* ---------------------------------------------------------------------- */ + +void FixFFL::initial_integrate_respa(int vflag, int ilevel, int iloop) +{ + dtv = step_respa[ilevel]; + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + + // innermost level - NVE update of v and x + // all other levels - NVE update of v + + if (ilevel==nlevels_respa-1) ffl_integrate(); + doffl=0; + if (ilevel == 0) initial_integrate(vflag); + else + { + final_integrate(); + } +} + +void FixFFL::final_integrate_respa(int ilevel, int iloop) +{ + + dtv = step_respa[ilevel]; + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + doffl=0; + final_integrate(); + if (ilevel==nlevels_respa-1) ffl_integrate(); +} + + +double FixFFL::compute_scalar() +{ + + double energy_me = energy; + double energy_all; + MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); + + return energy_all; +} + +/* ---------------------------------------------------------------------- + extract thermostat properties +------------------------------------------------------------------------- */ + +void *FixFFL::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"t_target") == 0) + { + return &t_target; + } + return NULL; +} + + +/* ---------------------------------------------------------------------- + Called when a change to the target temperature is requested mid-run +------------------------------------------------------------------------- */ + +void FixFFL::reset_target(double t_new) +{ + + t_target = t_start = t_stop = t_new; +} + +/* ---------------------------------------------------------------------- + Called when a change to the timestep is requested mid-run +------------------------------------------------------------------------- */ + +void FixFFL::reset_dt() +{ + // set the time integration constants + dtv = update->dt; + dtf = 0.5 * update->dt * (force->ftm2v); + init_ffl(); +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixFFL::memory_usage() +{ + double bytes = atom->nmax*(3*2)*sizeof(double); + return bytes; +} + + +/* ---------------------------------------------------------------------- + allocate local atom-based arrays +------------------------------------------------------------------------- */ + +void FixFFL::grow_arrays(int nmax) +{ + memory->grow(ffl_tmp1, nmax*3,"ffl:tmp1"); + memory->grow(ffl_tmp2, nmax*3,"ffl:tmp2"); + //zeroes out temporary buffers + for (int i=0; i< nmax*3; ++i) ffl_tmp1[i] = 0.0; + for (int i=0; i< nmax*3; ++i) ffl_tmp2[i] = 0.0; +} + + diff --git a/src/USER-MISC/fix_ffl.h b/src/USER-MISC/fix_ffl.h new file mode 100644 index 0000000000..1225d25989 --- /dev/null +++ b/src/USER-MISC/fix_ffl.h @@ -0,0 +1,67 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(ffl,FixFFL) + +#else + +#ifndef LMP_FIX_FFL_H +#define LMP_FIX_FFL_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixFFL : public Fix { + public: + FixFFL(class LAMMPS *, int, char **); + virtual ~FixFFL(); + int setmask(); + void init(); + void setup(int); + void ffl_integrate(); + void initial_integrate_respa(int vflag, int ilevel, int iloop); + void final_integrate_respa(int ilevel, int iloop); + void initial_integrate(int vflag); + void final_integrate(); + double compute_scalar(); + void reset_target(double); + virtual void reset_dt(); + double memory_usage(); + void grow_arrays(int); + + virtual void *extract(const char *, int &); + + void init_ffl(); + protected: + double *ffl_tmp1, *ffl_tmp2; + double t_start, t_stop, t_target; + double dtv, dtf, c1, c2, gamma; + char flip_type[10]; + + int doffl, ffl_every, ffl_step, flip_int; + class RanMars *random; + double *sqrt_m; + double *step_respa; + double energy; + int nlevels_respa; + + double **vaux; +}; + +} + +#endif +#endif