initial import of adapted gREM code by David Stelter and Edyta Malolepsza

/compute_pe_tally.h

/compute_plasticity_atom.cpp

/compute_plasticity_atom.h

/compute_pressure_grem.cpp

/compute_pressure_grem.h

/compute_rigid_local.cpp

/compute_rigid_local.h

/compute_spec_atom.cpp

/fix_gle.h

/fix_gpu.cpp

/fix_gpu.h

/fix_grem.cpp

/fix_grem.h

/fix_imd.cpp

/fix_imd.h

/fix_ipi.cpp

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

   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.

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

#include <mpi.h>

#include <string.h>

#include <stdlib.h>

#include "compute_pressure_grem.h"

#include "atom.h"

#include "update.h"

#include "domain.h"

#include "modify.h"

#include "fix.h"

#include "force.h"

#include "pair.h"

#include "kspace.h"

#include "error.h"

using namespace LAMMPS_NS;

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

   Last argument is the id of the gREM fix

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

ComputePressureGrem::ComputePressureGrem(LAMMPS *lmp, int narg, char **arg) :

  ComputePressure(lmp, narg-1, arg)

{

  fix_grem = arg[narg-1];

}

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

void ComputePressureGrem::init()

{

  ComputePressure::init();

  // Initialize hook to gREM fix

  int ifix = modify->find_fix(fix_grem);

  if (ifix < 0)

    error->all(FLERR,"Fix grem ID for compute pressure/grem does not exist");

  int dim;

  scale_grem = (double *)modify->fix[ifix]->extract("scale_grem",dim);

  if (scale_grem == NULL || dim != 0)

    error->all(FLERR,"Cannot extract gREM scale factor from fix grem");

}

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

   compute total pressure, averaged over Pxx, Pyy, Pzz

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

double ComputePressureGrem::compute_scalar()

{

  invoked_scalar = update->ntimestep;

  if (update->vflag_global != invoked_scalar)

    error->all(FLERR,"Virial was not tallied on needed timestep");

  // invoke temperature if it hasn't been already

  double t;

  if (keflag) {

    if (temperature->invoked_scalar != update->ntimestep)

      t = temperature->compute_scalar() * (*scale_grem);

    else t = temperature->scalar * (*scale_grem);

  }

  if (dimension == 3) {

    inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);

    virial_compute(3,3);

    if (keflag)

      scalar = (temperature->dof * boltz * t +

                virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;

    else

      scalar = (virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;

  } else {

    inv_volume = 1.0 / (domain->xprd * domain->yprd);

    virial_compute(2,2);

    if (keflag)

      scalar = (temperature->dof * boltz * t +

                virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;

    else

      scalar = (virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;

  }

  return scalar;

}

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

   compute pressure tensor

   assume KE tensor has already been computed

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

void ComputePressureGrem::compute_vector()

{

  invoked_vector = update->ntimestep;

  if (update->vflag_global != invoked_vector)

    error->all(FLERR,"Virial was not tallied on needed timestep");

  if (force->kspace && kspace_virial && force->kspace->scalar_pressure_flag)

    error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' for "

	       "tensor components with kspace_style msm");

  // invoke temperature if it hasn't been already

  double *ke_tensor;

  if (keflag) {

    if (temperature->invoked_vector != update->ntimestep)

      temperature->compute_vector();

    ke_tensor = temperature->vector;

  }

  for (int i = 0; i < 6; i++)

    ke_tensor[i] *= *scale_grem;

  if (dimension == 3) {

    inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);

    virial_compute(6,3);

    if (keflag) {

      for (int i = 0; i < 6; i++)

        vector[i] = (ke_tensor[i] + virial[i]) * inv_volume * nktv2p;

    } else

      for (int i = 0; i < 6; i++)

        vector[i] = virial[i] * inv_volume * nktv2p;

  } else {

    inv_volume = 1.0 / (domain->xprd * domain->yprd);

    virial_compute(4,2);

    if (keflag) {

      vector[0] = (ke_tensor[0] + virial[0]) * inv_volume * nktv2p;

      vector[1] = (ke_tensor[1] + virial[1]) * inv_volume * nktv2p;

      vector[3] = (ke_tensor[3] + virial[3]) * inv_volume * nktv2p;

      vector[2] = vector[4] = vector[5] = 0.0;

    } else {

      vector[0] = virial[0] * inv_volume * nktv2p;

      vector[1] = virial[1] * inv_volume * nktv2p;

      vector[3] = virial[3] * inv_volume * nktv2p;

      vector[2] = vector[4] = vector[5] = 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 COMPUTE_CLASS

ComputeStyle(pressure/grem,ComputePressureGrem)

#else

#ifndef LMP_COMPUTE_PRESSURE_GREM_H

#define LMP_COMPUTE_PRESSURE_GREM_H

#include "compute_pressure.h"

namespace LAMMPS_NS {

class ComputePressureGrem : public ComputePressure {

 public:

  ComputePressureGrem(class LAMMPS *, int, char **);

  virtual ~ComputePressureGrem() {};

  virtual void init();

  virtual double compute_scalar();

  virtual void compute_vector();

 protected:

  // Access to gREM fix scale factor

  char   *fix_grem;

  double *scale_grem;

};

}

#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: Compute pressure must use group all

Virial contributions computed by potentials (pair, bond, etc) are

computed on all atoms.

E: Could not find compute pressure temperature ID

The compute ID for calculating temperature does not exist.

E: Compute pressure temperature ID does not compute temperature

The compute ID assigned to a pressure computation must compute

temperature.

E: Compute pressure requires temperature ID to include kinetic energy

The keflag cannot be used unless a temperature compute is provided.

E: Virial was not tallied on needed timestep

You are using a thermo keyword that requires potentials to

have tallied the virial, but they didn't on this timestep.  See the

variable doc page for ideas on how to make this work.

E: Must use 'kspace_modify pressure/scalar no' for tensor components with kspace_style msm

Otherwise MSM will compute only a scalar pressure.  See the kspace_modify

command for details on this setting.

E: Fix grem ID for compute pressure/grem does not exist

Compute pressure/grem was passed an invalid fix id

E: Cannot extract gREM scale factor from fix grem

The fix id passed to compute pressure/grem refers to an incompatible fix

*/

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

   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.

   Force scaling fix for gREM.

   Cite: http://dx.doi.org/10.1063/1.3432176

   Cite: http://dx.doi.org/10.1021/acs.jpcb.5b07614

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

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

   Contributing author: Edyta Malolepsza (Broad Institute)

   Contributing author: David Stelter (Boston University)

   Contributing author: Tom Keyes (Boston University)

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

#include <string.h>

#include <stdlib.h>

#include <math.h>

#include "comm.h"

#include "fix_grem.h"

#include "atom.h"

#include "force.h"

#include "update.h"

#include "modify.h"

#include "domain.h"

#include "input.h"

#include "compute.h"

#include "memory.h"

#include "error.h"

using namespace LAMMPS_NS;

using namespace FixConst;

enum{NONE,CONSTANT,EQUAL,ATOM};

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

FixGrem::FixGrem(LAMMPS *lmp, int narg, char **arg) :

  Fix(lmp, narg, arg)

{

  if (narg < 8) error->all(FLERR,"Illegal fix grem command");

  scalar_flag = 1;

  vector_flag = 1;

  size_vector = 3;

  global_freq = 1;

  extscalar = 1;

  extvector = 1;

  // tbath - temp of bath, the same as defined in thermostat

  lambda = force->numeric(FLERR,arg[3]);

  eta = force->numeric(FLERR,arg[4]);

  h0 = force->numeric(FLERR,arg[5]);

  tbath = force->numeric(FLERR,arg[6]);

  pressref = force->numeric(FLERR,arg[7]);

  // create a new compute temp style

  // id = fix-ID + temp

  // compute group = all since pressure is always global (group all)

  //   and thus its KE/temperature contribution should use group all

  int n = strlen(id) + 6;

  id_temp = new char[n];

  strcpy(id_temp,id);

  strcat(id_temp,"_temp");

  char **newarg = new char*[3];

  newarg[0] = id_temp;

  newarg[1] = (char *) "all";

  newarg[2] = (char *) "temp";

  modify->add_compute(3,newarg);

  delete [] newarg;

  tflag = 1;

  // create a new compute pressure style

  // id = fix-ID + press, compute group = all

  // pass id_temp as 4th arg to pressure constructor

  n = strlen(id) + 7;

  id_press = new char[n];

  strcpy(id_press,id);

  strcat(id_press,"_press");

  newarg = new char*[5];

  newarg[0] = id_press;

  newarg[1] = (char *) "all";

  newarg[2] = (char *) "pressure/grem";

  newarg[3] = id_temp;

  newarg[4] = id;

  modify->add_compute(5,newarg);

  delete [] newarg;

  // create a new compute ke style

  // id = fix-ID + ke

  n = strlen(id) + 8;

  id_ke = new char[n];

  strcpy(id_ke,id);

  strcat(id_ke,"_ke");

  newarg = new char*[3];

  newarg[0] = id_ke;

  newarg[1] = (char *) "all";

  newarg[2] = (char *) "ke";

  modify->add_compute(3,newarg);

  delete [] newarg;

  keflag = 1;

  // create a new compute pe style

  // id = fix-ID + pe

  n = strlen(id) + 9;

  id_pe = new char[n];

  strcpy(id_pe,id);

  strcat(id_pe,"_pe");

  newarg = new char*[3];

  newarg[0] = id_pe;

  newarg[1] = (char *) "all";

  newarg[2] = (char *) "pe";

  modify->add_compute(3,newarg);

  delete [] newarg;

  peflag = 1;

}

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

FixGrem::~FixGrem()

{

  // delete temperature, pressure and energies if fix created them

  if (tflag) modify->delete_compute(id_temp);

  if (pflag) modify->delete_compute(id_press);

  if (keflag) modify->delete_compute(id_ke);

  if (peflag) modify->delete_compute(id_pe);

  delete [] id_temp;

  delete [] id_press;

  delete [] id_ke;

  delete [] id_pe;

}

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

int FixGrem::setmask()

{

  int mask = 0;

  mask |= POST_FORCE;

  return mask;

}

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

void FixGrem::init()

{

  // add check of sign of eta

  // set temperature and pressure ptrs

  int icompute = modify->find_compute(id_temp);

  if (icompute < 0)

    error->all(FLERR,"Temperature ID for fix scaledforce does not exist");

  temperature = modify->compute[icompute];

  icompute = modify->find_compute(id_ke);

  if (icompute < 0)

    error->all(FLERR,"Ke ID for fix scaledforce does not exist");

  ke = modify->compute[icompute];

  icompute = modify->find_compute(id_pe);

  if (icompute < 0)

    error->all(FLERR,"Pe ID for fix scaledforce does not exist");

  pe = modify->compute[icompute];  

}

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

void FixGrem::setup(int vflag)

{

  if (strstr(update->integrate_style,"verlet"))

    post_force(vflag);

}

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

void FixGrem::min_setup(int vflag)

{

  post_force(vflag);

}

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

void FixGrem::post_force(int vflag)

{

  double **f = atom->f;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  double tmpvolume = domain->xprd * domain->yprd * domain->zprd;

  double tmppe = pe->compute_scalar();

  // potential energy

  double tmpenthalpy = tmppe+pressref*tmpvolume/(force->nktv2p);

  double teffective = lambda+eta*(tmpenthalpy-h0);

  scale_grem = tbath/teffective;

  for (int i = 0; i < nlocal; i++)

    if (mask[i] & groupbit) {

      f[i][0] *= scale_grem;

      f[i][1] *= scale_grem;

      f[i][2] *= scale_grem;

    }

  pe->addstep(update->ntimestep+1);

}

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

   extract scale factor

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

void *FixGrem::extract(const char *str, int &dim)

{

  dim=0;

  if (strcmp(str,"scale_grem") == 0) {

    return &scale_grem;

  }

  return NULL;

}

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

   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(grem,FixGrem)

#else

#ifndef LMP_FIX_GREM_H

#define LMP_FIX_GREM_H

#include "fix.h"

namespace LAMMPS_NS {

class FixGrem : public Fix {

 public:

  FixGrem(class LAMMPS *, int, char **);

  ~FixGrem();

  int setmask();

  void init();

  void setup(int);

  void min_setup(int);

  void post_force(int);

  void *extract(const char *, int &);

  double scale_grem;

 private:

  double lambda,eta,h0,tbath,pressref;

 protected:

  char *id_temp,*id_press,*id_ke,*id_pe;

  class Compute *temperature,*pressure,*ke,*pe;

  int pflag,tflag,keflag,peflag;

};

}

#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: Region ID for fix grem does not exist

Self-explanatory.

E: Variable name for fix grem does not exist

Self-explanatory.

E: Variable for fix grem is invalid style

Self-explanatory.

E: Cannot use variable energy with constant force in fix grem

This is because for constant force, LAMMPS can compute the change

in energy directly.

E: Must use variable energy with fix grem

Must define an energy vartiable when applyting a dynamic

force during minimization.

*/