Merge pull request #310 from EfremBraun/master

"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,

"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).

   fix_nve_noforce

   fix_nve_sphere

   fix_nve_tri

   fix_nvk

   fix_nvt_asphere

   fix_nvt_body

   fix_nvt_manifold_rattle

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

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

   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;

}