Fast-forward Langevin functionality included

<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)"_#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)

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

/* ----------------------------------------------------------------------

   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 <math.h>

#include <string.h>

#include <stdlib.h>

#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 <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) 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;

}

/* -*- 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;

};

}