diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt index 20503e8522382f00261f2ef0d13baa3ef5abe597..51c199f27285098445d138d6e8e8b39f702bff08 100644 --- a/doc/src/Section_commands.txt +++ b/doc/src/Section_commands.txt @@ -702,6 +702,7 @@ package"_Section_start.html#start_3. "manifoldforce"_fix_manifoldforce.html, "meso/stationary"_fix_meso_stationary.html, "nve/manifold/rattle"_fix_nve_manifold_rattle.html, +"nvk"_fix_nvk.html, "nvt/manifold/rattle"_fix_nvt_manifold_rattle.html, "nph/eff"_fix_nh_eff.html, "npt/eff"_fix_nh_eff.html, diff --git a/doc/src/fix_nvk.txt b/doc/src/fix_nvk.txt new file mode 100644 index 0000000000000000000000000000000000000000..271483b44188f04a503777c20d1a0c01a6dc0c6e --- /dev/null +++ b/doc/src/fix_nvk.txt @@ -0,0 +1,71 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix nvk command :h3 + +[Syntax:] + +fix ID group-ID nvk :pre + +ID, group-ID are documented in "fix"_fix.html command +nvk = style name of this fix command :ul + +[Examples:] + +fix 1 all nvk :pre + +[Description:] + +Perform constant kinetic energy integration using the Gaussian +thermostat to update position and velocity for atoms in the group each +timestep. V is volume; K is kinetic energy. This creates a system +trajectory consistent with the isokinetic ensemble. + +The equations of motion used are those of Minary et al in +"(Minary)"_#nvk-Minary, a variant of those initially given by Zhang in +"(Zhang)"_#nvk-Zhang. + +The kinetic energy will be held constant at its value given when fix +nvk is initiated. If a different kinetic energy is desired, the +"velocity"_velocity.html command should be used to change the kinetic +energy prior to this fix. + +:line + +[Restart, fix_modify, output, run start/stop, minimize info:] + +No information about this fix is written to "binary restart +files"_restart.html. None of the "fix_modify"_fix_modify.html options +are relevant to this fix. No global or per-atom quantities are stored +by this fix for access by various "output +commands"_Section_howto.html#howto_15. No parameter of this fix can +be used with the {start/stop} keywords of the "run"_run.html command. +This fix is not invoked during "energy minimization"_minimize.html. + +[Restrictions:] + +The Gaussian thermostat only works when it is applied to all atoms in +the simulation box. Therefore, the group must be set to all. + +This fix has not yet been implemented to work with the RESPA integrator. + +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:] none + +[Default:] none + +:line + +:link(nvk-Minary) +[(Minary)] Minary, Martyna, and Tuckerman, J Chem Phys, 18, 2510 (2003). + +:link(nvk-Zhang) +[(Zhang)] Zhang, J Chem Phys, 106, 6102 (1997). diff --git a/doc/src/fixes.txt b/doc/src/fixes.txt index d17b88306ddd681e0e1e30bc8f6700cad5af6dc3..ca43a1d1e71a0407ba5ff41702074da487b79f1c 100644 --- a/doc/src/fixes.txt +++ b/doc/src/fixes.txt @@ -90,6 +90,7 @@ Fixes :h1 fix_nve_noforce fix_nve_sphere fix_nve_tri + fix_nvk fix_nvt_asphere fix_nvt_body fix_nvt_manifold_rattle diff --git a/src/USER-MISC/README b/src/USER-MISC/README index 513973b4f56db591051be1f3ec1267e4d7140930..8aed7972355e348e32d03484670a6abad28490f4 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -43,6 +43,7 @@ fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov fix grem, David Stelter, dstelter@bu.edu, 22 Nov 16 fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009 fix ipi, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014 +fix nvk, Efrem Braun (UC Berkeley), efrem.braun at gmail.com, https://github.com/lammps/lammps/pull/310 fix pimd, Yuxing Peng (U Chicago), yuxing at uchicago.edu, 24 Nov 2014 fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008 fix ti/spring, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013 diff --git a/src/USER-MISC/fix_nvk.cpp b/src/USER-MISC/fix_nvk.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5dcfe53e788c5f23539cd7f4fbfb764933112e2a --- /dev/null +++ b/src/USER-MISC/fix_nvk.cpp @@ -0,0 +1,219 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Efrem Braun (UC Berkeley) +------------------------------------------------------------------------- */ + +#include <math.h> +#include <stdio.h> +#include <string.h> +#include "fix_nvk.h" +#include "atom.h" +#include "force.h" +#include "update.h" +#include "respa.h" +#include "error.h" +#include "compute.h" +#include "math_extra.h" +#include "domain.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +FixNVK::FixNVK(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 3) + error->all(FLERR,"Illegal fix nvk command"); + if (igroup) error->all(FLERR,"Fix nvk only supports group all"); + + dynamic_group_allow = 1; + time_integrate = 1; +} + +/* ---------------------------------------------------------------------- */ + +int FixNVK::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= INITIAL_INTEGRATE_RESPA; + mask |= FINAL_INTEGRATE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixNVK::init() +{ + dtv = update->dt; + dtf = 0.5 * update->dt; + + if (strstr(update->integrate_style,"respa")) { + error->all(FLERR,"Fix nvk not yet enabled for RESPA"); + step_respa = ((Respa *) update->integrate)->step; + } + + // compute initial kinetic energy + // make better by calling compute_ke instead of copy/pasting code from compute_ke.cpp + double pfactor = 0.5 * force->mvv2e; + double **v = atom->v; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + double ke = 0.0; + if (rmass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + ke += rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + ke += mass[type[i]] * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + } + MPI_Allreduce(&ke,&K_target,1,MPI_DOUBLE,MPI_SUM,world); + K_target *= pfactor; +} + +/* ---------------------------------------------------------------------- + allow for both per-type and per-atom mass +------------------------------------------------------------------------- */ + +void FixNVK::initial_integrate(int vflag) +{ + double sm; + double a,b,sqtb,s,sdot; + + 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; + + // calculate s and sdot from Minary 2003, equations 4.12 and 4.13 + double a_local = 0.0; + double b_local = 0.0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + a_local += MathExtra::dot3(f[i], v[i]); + if (rmass) b_local += MathExtra::dot3(f[i], f[i]) / rmass[i]; + else b_local += MathExtra::dot3(f[i], f[i]) / mass[type[i]]; + } + MPI_Allreduce(&a_local,&a,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&b_local,&b,1,MPI_DOUBLE,MPI_SUM,world); + a /= (2.0*K_target); // units of inverse time + b /= (2.0*K_target * force->mvv2e); // units of inverse time squared + sqtb = sqrt(b); + s = a/b * (cosh(dtf*sqtb) - 1.0) + sinh(dtf*sqtb) / sqtb; + sdot = a/b * sqtb * sinh(dtf*sqtb) + cosh(dtf*sqtb); + + // update v and x of atoms in group per Minary 2003, equations 4.15-4.17 + // note that equation 4.15, 4.17 should read p = (p+F*s/m)/sdot + // note that equation 4.16 should read r = r + delt*p/m + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (rmass) sm = s / rmass[i]; + else sm = s / mass[type[i]]; + v[i][0] = (v[i][0] + f[i][0] * sm * force->ftm2v) / sdot; + v[i][1] = (v[i][1] + f[i][1] * sm * force->ftm2v) / sdot; + v[i][2] = (v[i][2] + f[i][2] * sm * force->ftm2v) / sdot; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNVK::final_integrate() +{ + double sm; + double a,b,sqtb,s,sdot; + + 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; + + // calculate s and sdot from Minary 2003, equations 4.12 and 4.13 + double a_local = 0.0; + double b_local = 0.0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + a_local += MathExtra::dot3(f[i], v[i]); + if (rmass) b_local += MathExtra::dot3(f[i], f[i]) / rmass[i]; + else b_local += MathExtra::dot3(f[i], f[i]) / mass[type[i]]; + } + MPI_Allreduce(&a_local,&a,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&b_local,&b,1,MPI_DOUBLE,MPI_SUM,world); + a /= (2.0*K_target); // units of inverse time + b /= (2.0*K_target * force->mvv2e); // units of inverse time squared + sqtb = sqrt(b); + s = a/b * (cosh(dtf*sqtb) - 1.0) + sinh(dtf*sqtb) / sqtb; + sdot = a/b * sqtb * sinh(dtf*sqtb) + cosh(dtf*sqtb); + + // update v and x of atoms in group per Minary 2003, equations 4.15-4.17 + // note that equation 4.15, 4.17 should read p = (p+F*s/m)/sdot + // note that equation 4.16 should read r = r + delt*p/m + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (rmass) sm = s / rmass[i]; + else sm = s / mass[type[i]]; + v[i][0] = (v[i][0] + f[i][0] * sm * force->ftm2v) / sdot; + v[i][1] = (v[i][1] + f[i][1] * sm * force->ftm2v) / sdot; + v[i][2] = (v[i][2] + f[i][2] * sm * force->ftm2v) / sdot; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNVK::initial_integrate_respa(int vflag, int ilevel, int iloop) +{ + dtv = step_respa[ilevel]; + dtf = 0.5 * step_respa[ilevel]; + + // innermost level - NVK update of v and x + // all other levels - NVK update of v + + if (ilevel == 0) initial_integrate(vflag); + else final_integrate(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNVK::final_integrate_respa(int ilevel, int iloop) +{ + dtf = 0.5 * step_respa[ilevel]; + final_integrate(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNVK::reset_dt() +{ + dtv = update->dt; + dtf = 0.5 * update->dt; +} diff --git a/src/USER-MISC/fix_nvk.h b/src/USER-MISC/fix_nvk.h new file mode 100644 index 0000000000000000000000000000000000000000..1d73e2f5dfa3bafd15043fe00ed69357c299ec91 --- /dev/null +++ b/src/USER-MISC/fix_nvk.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(nvk,FixNVK) + +#else + +#ifndef LMP_FIX_NVK_H +#define LMP_FIX_NVK_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixNVK : public Fix { + public: + FixNVK(class LAMMPS *, int, char **); + virtual ~FixNVK() {} + int setmask(); + virtual void init(); + virtual void initial_integrate(int); + virtual void final_integrate(); + virtual void initial_integrate_respa(int, int, int); + virtual void final_integrate_respa(int, int); + virtual void reset_dt(); + + protected: + double dtv,dtf; + double *step_respa; + int mass_require; + double K_target; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Fix nvk only supports group all + +Self-explanatory. + +E: Fix nvk not yet enabled for RESPA + +Self-explanatory. + +*/