add package and examples for viscosity calculation

# DFF generated Lammps input file

units          real

atom_style     full

boundary       p p p

pair_style     lj/cut/coul/long     10.0

pair_modify    mix arithmetic

pair_modify    tail yes

kspace_style   pppm 1.0e-4

dielectric     1.0

special_bonds  amber

bond_style     harmonic

angle_style    harmonic

dihedral_style none

improper_style none

read_data      data.1000SPCE.lmp

variable T equal 300

variable P equal 1.0

velocity       all create ${T} 12345 mom yes rot yes dist gaussian

timestep       1.0

# Constraint ##################################

fix com        all momentum 100 linear 1 1 1

fix rigid      all shake 1e-4 20 0 b 1 a 1

# Viscosity  ##################################

variable       A equal 0.02e-5 # angstrom/fs^2

fix            cos all accelerate/cos ${A}

compute        cos all viscosity/cos

variable       density equal density

variable       lz equal lz

variable       vMax equal c_cos[7] # velocity of atoms at z=0

variable       invVis equal v_vMax/${A}/v_density*39.4784/v_lz/v_lz*100 # reciprocal of viscosity 1/Pa/s

fix            npt all npt temp ${T} ${T} 100 iso ${P} ${P} 1000

fix_modify     npt temp cos

thermo_style   custom step cpu temp press pe density v_vMax v_invVis

thermo_modify  temp cos

thermo         100

################################################

dump           1 all custom 10000 dump.lammpstrj id mol type element q xu yu zu

dump_modify    1 sort id element O H

run            2000

This package provides fix and compute styles for calculating viscosity

with the periodic perturbation method, as described in the following paper:

    Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209–217.

There is example script for using this package in examples/USER/viscosity

The person who created this package is Zheng GONG (z.gong@outlook.com)

at ENS de Lyon.

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

   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: Zheng GONG (ENS Lyon)

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

#include <mpi.h>

#include "compute_viscosity_cos.h"

#include "atom.h"

#include "update.h"

#include "force.h"

#include "domain.h"

#include "comm.h"

#include "group.h"

#include "error.h"

#include "math_const.h"

using namespace LAMMPS_NS;

using namespace MathConst;

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

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

    Compute(lmp, narg, arg) {

  if (narg != 3) error->all(FLERR, "Illegal compute viscosity/cos command");

  scalar_flag = vector_flag = 1;

  size_vector = 7;

  extscalar = 0;

  extvector = 1;

  tempflag = 1;

  tempbias = 1;

  maxbias = 0;

  vbiasall = NULL;

  vector = new double[7];

}

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

ComputeViscosityCos::~ComputeViscosityCos() {

  if (!copymode)

    delete[] vector;

}

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

void ComputeViscosityCos::setup() {

  dynamic = 0;

  if (dynamic_user || group->dynamic[igroup]) dynamic = 1;

  dof_compute();

}

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

void ComputeViscosityCos::dof_compute() {

  adjust_dof_fix();

  natoms_temp = group->count(igroup);

  dof = domain->dimension * natoms_temp;

  dof -= extra_dof + fix_dof;

  if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);

  else tfactor = 0.0;

}

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

void ComputeViscosityCos::calc_V() {

  double zlo = domain->boxlo[2];

  double zhi = domain->boxhi[2];

  double **x = atom->x;

  double **v = atom->v;

  double *mass = atom->mass;

  double *rmass = atom->rmass;

  int *type = atom->type;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  double massone;

  double V_m[2];

  double V_m_local[2] = {0, 0};

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

    if (mask[i] & groupbit) {

      if (rmass) massone = rmass[i];

      else massone = mass[type[i]];

      V_m_local[0] +=

          2 * massone * v[i][0] * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));

      V_m_local[1] += massone;

    }

  MPI_Allreduce(V_m_local, V_m, 2, MPI_DOUBLE, MPI_SUM, world);

  V = V_m[0] / V_m[1];

}

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

double ComputeViscosityCos::compute_scalar() {

  invoked_scalar = update->ntimestep;

  double zlo = domain->boxlo[2];

  double zhi = domain->boxhi[2];

  double **x = atom->x;

  double **v = atom->v;

  double *mass = atom->mass;

  double *rmass = atom->rmass;

  int *type = atom->type;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  double t = 0.0;

  double vx_acc;

  double massone;

  calc_V();

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

    if (mask[i] & groupbit) {

      if (rmass) massone = rmass[i];

      else massone = mass[type[i]];

      vx_acc = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));

      t += ((v[i][0] - vx_acc) * (v[i][0] - vx_acc) + v[i][1] * v[i][1] +

            v[i][2] * v[i][2]) * massone;

    }

  MPI_Allreduce(&t, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);

  if (dynamic) dof_compute();

  if (dof < 0.0 && natoms_temp > 0.0)

    error->all(FLERR, "Temperature compute degrees of freedom < 0");

  scalar *= tfactor;

  return scalar;

}

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

void ComputeViscosityCos::compute_vector() {

  int i;

  invoked_vector = update->ntimestep;

  double zlo = domain->boxlo[2];

  double zhi = domain->boxhi[2];

  double **x = atom->x;

  double **v = atom->v;

  double *mass = atom->mass;

  double *rmass = atom->rmass;

  int *type = atom->type;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  double vx_acc;

  double massone, t[6];

  for (i = 0; i < 6; i++) t[i] = 0.0;

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

    if (mask[i] & groupbit) {

      if (rmass) massone = rmass[i];

      else massone = mass[type[i]];

      vx_acc = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));

      t[0] += massone * (v[i][0] - vx_acc) * (v[i][0] - vx_acc);

      t[1] += massone * v[i][1] * v[i][1];

      t[2] += massone * v[i][2] * v[i][2];

      t[3] += massone * (v[i][0] - vx_acc) * v[i][1];

      t[4] += massone * (v[i][0] - vx_acc) * v[i][2];

      t[5] += massone * v[i][1] * v[i][2];

    }

  MPI_Allreduce(t, vector, 6, MPI_DOUBLE, MPI_SUM, world);

  for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;

  vector[6] = V;

}

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

   remove velocity bias from atom I to leave thermal velocity

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

void ComputeViscosityCos::remove_bias(int i, double *v) {

  double zlo = domain->boxlo[2];

  double zhi = domain->boxhi[2];

  double **x = atom->x;

  vbias[0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));

  vbias[1] = 0;

  vbias[2] = 0;

  v[0] -= vbias[0];

//  v[1] -= vbias[1];

//  v[2] -= vbias[2];

}

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

   remove velocity bias from atom I to leave thermal velocity

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

void ComputeViscosityCos::remove_bias_thr(int i, double *v, double *b) {

  double zlo = domain->boxlo[2];

  double zhi = domain->boxhi[2];

  double **x = atom->x;

  b[0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));

  b[1] = 0;

  b[2] = 0;

  v[0] -= b[0];

//  v[1] -= b[1];

//  v[2] -= b[2];

}

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

   remove velocity bias from all atoms to leave thermal velocity

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

void ComputeViscosityCos::remove_bias_all() {

  double zlo = domain->boxlo[2];

  double zhi = domain->boxhi[2];

  double **x = atom->x;

  double **v = atom->v;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

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

    if (mask[i] & groupbit) {

      vbiasall[i][0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));

      vbiasall[i][1] = 0;

      vbiasall[i][2] = 0;

      v[i][0] -= vbiasall[i][0];

//      v[i][1] -= vbiasall[i][1];

//      v[i][2] -= vbiasall[i][2];

    }

}

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

   add back in velocity bias to atom I removed by remove_bias()

   assume remove_bias() was previously called

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

void ComputeViscosityCos::restore_bias(int i, double *v) {

  v[0] += vbias[0];

  v[1] += vbias[1];

  v[2] += vbias[2];

}

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

   add back in velocity bias to atom I removed by remove_bias_thr()

   assume remove_bias_thr() was previously called with the same buffer b

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

void ComputeViscosityCos::restore_bias_thr(int i, double *v, double *b) {

  v[0] += b[0];

  v[1] += b[1];

  v[2] += b[2];

}

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

   add back in velocity bias to all atoms removed by remove_bias_all()

   assume remove_bias_all() was previously called

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

void ComputeViscosityCos::restore_bias_all() {

  double **v = atom->v;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

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

    if (mask[i] & groupbit) {

      v[i][0] += vbiasall[i][0];

      v[i][1] += vbiasall[i][1];

      v[i][2] += vbiasall[i][2];

    }

}

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

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

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

   Contributing author: Zheng GONG (ENS Lyon)

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

#ifdef COMPUTE_CLASS

ComputeStyle(viscosity/cos,ComputeViscosityCos)

#else

#ifndef LMP_COMPUTE_VISCOSITY_COS_H

#define LMP_COMPUTE_VISCOSITY_COS_H

#include "compute.h"

namespace LAMMPS_NS {

class ComputeViscosityCos : public Compute {

 public:

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

  virtual ~ComputeViscosityCos();

  void init() {}

  void setup();

  virtual double compute_scalar();

  virtual void compute_vector();

  void remove_bias(int, double *);

  void remove_bias_thr(int, double *, double *);

  void remove_bias_all();

  void restore_bias(int, double *);

  void restore_bias_thr(int, double *, double *);

  void restore_bias_all();

 protected:

  double tfactor;

  double V;

  void dof_compute();

  void calc_V();

};

}

#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: Temperature compute degrees of freedom < 0

This should not happen if you are calculating the temperature

on a valid set of atoms.

*/