diff --git a/doc/src/Commands_fix.txt b/doc/src/Commands_fix.txt
index e002c117701110fc6dd33e83865386f52d2bae3e..2b6f4a3080f586c1d713a71db61db15b53536fee 100644
--- a/doc/src/Commands_fix.txt
+++ b/doc/src/Commands_fix.txt
@@ -65,6 +65,7 @@ OPT.
"eos/table/rx (k)"_fix_eos_table_rx.html,
"evaporate"_fix_evaporate.html,
"external"_fix_external.html,
+"ffl"_fix_ffl.html,
"filter/corotate"_fix_filter_corotate.html,
"flow/gauss"_fix_flow_gauss.html,
"freeze"_fix_freeze.html,
diff --git a/doc/src/fix_ffl.txt b/doc/src/fix_ffl.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ea33dc432f749dbb636ef245323e1390644c0d5f
--- /dev/null
+++ b/doc/src/fix_ffl.txt
@@ -0,0 +1,124 @@
+<script type="text/javascript"
+ src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
+</script>
+<script type="text/x-mathjax-config">
+ MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
+</script>
+
+"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)"_#Bussi3. 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"_Howto_output.html. The scalar is the cumulative
+energy change due to this fix. The scalar value calculated by this
+fix is "extensive".
+
+[Restrictions:]
+
+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 "Build
+package"_Build_package.html doc page 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(Bussi3)
+[(Bussi)] G. Bussi, M. Parrinello, Phs. Rev. E 75, 056707 (2007)
+
diff --git a/doc/src/fixes.txt b/doc/src/fixes.txt
index f2f257b1346240f56256dbd4107198f437d30e27..9ab299465b9b0b129ffca7b886a166e27d16b984 100644
--- a/doc/src/fixes.txt
+++ b/doc/src/fixes.txt
@@ -46,6 +46,7 @@ Fixes :h1
fix_eos_table_rx
fix_evaporate
fix_external
+ fix_ffl
fix_filter_corotate
fix_flow_gauss
fix_freeze
diff --git a/doc/src/lammps.book b/doc/src/lammps.book
index bc5df982b70c02a666c30867623ae774d9af5c01..2efa690e18cd5a20c55d841d0aa65890b87dc589 100644
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@@ -266,6 +266,7 @@ fix_eos_table.html
fix_eos_table_rx.html
fix_evaporate.html
fix_external.html
+fix_ffl.html
fix_filter_corotate.html
fix_flow_gauss.html
fix_freeze.html
diff --git a/src/.gitignore b/src/.gitignore
index 02cb90fd414963d1ced4b23e2b4ad0f6b5d07694..e3c273ebefcbebc6d23495c3fdff0616008b7e00 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -429,6 +429,8 @@
/fix_eos_table.h
/fix_evaporate.cpp
/fix_evaporate.h
+/fix_ffl.cpp
+/fix_ffl.h
/fix_filter_corotate.cpp
/fix_filter_corotate.h
/fix_viscosity.cpp
diff --git a/src/USER-MISC/README b/src/USER-MISC/README
index eb221c07db67bbd8569ac7d61f4594af29b76654..b3c3e859480b0ba3e6bfbb1a501d684f02d284b3 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 0000000000000000000000000000000000000000..e2fb0166a28473be73e33ddba4e5c746b2d6679d
--- /dev/null
+++ b/src/USER-MISC/fix_ffl.cpp
@@ -0,0 +1,483 @@
+/* ----------------------------------------------------------------------
+ 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 <mpi.h>
+#include <cmath>
+#include <cstring>
+#include <cstdlib>
+#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};
+enum {NO_FLIP, FLIP_RESCALE, FLIP_HARD, FLIP_SOFT};
+//#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 <fix-ID>"
+ " <group-ID> ffl <tau> <Tstart> <Tstop> <seed> ");
+
+ 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) {
+ if (strcmp(arg[7],"no_flip") == 0) {
+ flip_int = NO_FLIP;
+ } else if (strcmp(arg[7],"rescale") == 0) {
+ flip_int = FLIP_RESCALE;
+ } else if (strcmp(arg[7],"hard") == 0) {
+ flip_int = FLIP_HARD;
+ } else if (strcmp(arg[7],"soft") == 0) {
+ flip_int = FLIP_SOFT;
+ } else {
+ error->all(FLERR,"Illegal fix ffl flip type, only accepts : rescale - hard - soft - no_flip");
+ }
+ } else {
+ flip_int = FLIP_RESCALE;
+ }
+
+ 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 == FLIP_HARD) {
+ if (v[i][k]*ffl_tmp2[nk] < 0.0) v[i][k] = -v[i][k];
+ }
+
+ nk++;
+ }
+ }
+
+ //rescale operation (RESCALE FLIP)
+ if (flip_int == FLIP_RESCALE) {
+ 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 == FLIP_SOFT) {
+ 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 0000000000000000000000000000000000000000..1225d2598940e5a57208f73f0b4d24568abe0c92
--- /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