Add documents and minor fix for USER-VISCOSITY package

doc/src/Commands_compute.rst

@@ -161,5 +161,6 @@ KOKKOS, o = USER-OMP, t = OPT.
    * :doc:`torque/chunk <compute_torque_chunk>`
    * :doc:`vacf <compute_vacf>`
    * :doc:`vcm/chunk <compute_vcm_chunk>`
+   * :doc:`viscosity/cos <compute_viscosity_cos>`
    * :doc:`voronoi/atom <compute_voronoi_atom>`
    * :doc:`xrd <compute_xrd>`

doc/src/Commands_fix.rst

@@ -22,6 +22,7 @@ OPT.
 .. table_from_list::
    :columns: 5
 
+   * :doc:`accelerate/cos <fix_accelerate_cos>`
    * :doc:`adapt <fix_adapt>`
    * :doc:`adapt/fep <fix_adapt_fep>`
    * :doc:`addforce <fix_addforce>`

doc/src/Packages_details.rst

@@ -104,6 +104,7 @@ page gives those details.
    * :ref:`USER-UEF <PKG-USER-UEF>`
    * :ref:`USER-VTK <PKG-USER-VTK>`
    * :ref:`USER-YAFF <PKG-USER-YAFF>`
+   * :ref:`USER-VISCOSITY <PKG-USER-VISCOSITY>`
 
 ----------
 
@@ -2365,3 +2366,29 @@ which discuss the `QuickFF <quickff_>`_ methodology.
 * :doc:`pair_style mm3/switch3/coulgauss/long <pair_mm3_switch3_coulgauss_long>`
 * :doc:`pair_style lj/switch3/coulgauss/long <pair_lj_switch3_coulgauss_long>`
 * examples/USER/yaff
+
+
+----------
+
+
+.. _PKG-USER-VISCOSITY:
+
+USER-VISCOSITY package
+----------------------
+
+**Contents:**
+
+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.
+
+**Author:** Zheng Gong (ENS de Lyon)
+
+**Supporting info:**
+
+* src/USER-VISCOSITY: filenames -> commands
+* src/USER-VISCOSITY/README
+* :doc:`fix accelerate/cos <fix_accelerate_cos>`
+* :doc:`compute viscosity/cos <compute_viscosity_cos>`
+* examples/USER/viscosity

doc/src/compute.rst

@@ -307,6 +307,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`torque/chunk <compute_torque_chunk>` - torque applied on each chunk
 * :doc:`vacf <compute_vacf>` - velocity auto-correlation function of group of atoms
 * :doc:`vcm/chunk <compute_vcm_chunk>` - velocity of center-of-mass for each chunk
+* :doc:`viscosity/cos <compute_viscosity_cos>` - velocity profile under cosine-shaped acceleration
 * :doc:`voronoi/atom <compute_voronoi_atom>` - Voronoi volume and neighbors for each atom
 * :doc:`xrd <compute_xrd>` - x-ray diffraction intensity on a mesh of reciprocal lattice nodes
 

doc/src/compute_viscosity_cos.rst

@@ -0,0 +1,157 @@
+.. index:: compute viscosity/cos
+
+compute viscosity/cos command
+=============================
+
+Syntax
+""""""
+
+
+.. parsed-literal::
+
+   compute ID group-ID viscosity/cos
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* viscosity/cos = style name of this compute command
+
+
+Examples
+""""""""
+
+
+.. parsed-literal::
+
+   compute  cos all viscosity/cos
+   variable V equal c_cos[7]
+   variable A equal 0.02E-5
+   variable density equal density
+   variable lz equal lz
+   variable reciprocalViscosity equal v_V/${A}/v_density*39.4784/v_lz/v_lz*100
+
+Description
+"""""""""""
+
+Define a computation that calculates the velocity amplitude of a group of atoms
+with an cosine-shaped velocity profile and the temperature of them
+after subtracting out the velocity profile before computing the kinetic energy.
+A compute of this style can be used by any command that computes a temperature,
+e.g. :doc:`thermo_modify <thermo_modify>`, :doc:`fix npt <fix_nh>`, etc.
+
+This command together with :doc:`fix_accelerate/cos<fix_accelerate_cos>`
+enables viscosity calculation with periodic perturbation method,
+as described by :ref:`Hess<Hess>`.
+An acceleration along the x-direction is applied to the simulation system
+by using :doc:`fix_accelerate/cos<fix_accelerate_cos>` command.
+The acceleration is a periodic function along the z-direction:
+
+.. math::
+
+   a_{x}(z) = A \cos \left(\frac{2 \pi z}{l_{z}}\right)
+
+where :math:`A` is the acceleration amplitude, :math:`l_z` is the z-length
+of the simulation box. At steady state, the acceleration generates
+a velocity profile:
+
+.. math::
+
+   v_{x}(z) = V \cos \left(\frac{2 \pi z}{l_{z}}\right)
+
+The generated velocity amplitude :math:`V` is related to the
+shear viscosity :math:`\eta` by:
+
+.. math::
+
+   V = \frac{A \rho}{\eta}\left(\frac{l_{z}}{2 \pi}\right)^{2}
+
+
+and it can be obtained from ensemble average of the velocity profile:
+
+.. math::
+
+   V = \frac{\sum_i 2 m_{i} v_{i, x} \cos \left(\frac{2 \pi z_i}{l_{z}}\right)}{\sum_i m_{i}}
+
+
+where :math:`m_i`, :math:`v_{i,x}` and :math:`z_i` are the mass,
+x-component velocity and z coordinate of a particle.
+
+After the cosine-shaped collective velocity in :math:`x` direction
+has been subtracted for each atom, the temperature is calculated by the formula
+KE = dim/2 N k T, where KE = total kinetic energy of the group of
+atoms (sum of 1/2 m v\^2), dim = 2 or 3 = dimensionality of the
+simulation, N = number of atoms in the group, k = Boltzmann constant,
+and T = temperature.
+
+A kinetic energy tensor, stored as a 6-element vector, is also
+calculated by this compute for use in the computation of a pressure
+tensor. The formula for the components of the tensor is the same as
+the above formula, except that v\^2 is replaced by vx\*vy for the xy
+component, etc. The 6 components of the vector are ordered xx, yy,
+zz, xy, xz, yz.
+
+The number of atoms contributing to the temperature is assumed to be
+constant for the duration of the run; use the *dynamic* option of the
+:doc:`compute_modify <compute_modify>` command if this is not the case.
+However, in order to get meaningful result, the group ID of this compute should be all.
+
+The removal of the cosine-shaped velocity component by this command is
+essentially computing the temperature after a "bias" has been removed
+from the velocity of the atoms.  If this compute is used with a fix
+command that performs thermostatting then this bias will be subtracted
+from each atom, thermostatting of the remaining thermal velocity will
+be performed, and the bias will be added back in.  Thermostatting
+fixes that work in this way include :doc:`fix nvt <fix_nh>`, :doc:`fix temp/rescale <fix_temp_rescale>`, :doc:`fix temp/berendsen <fix_temp_berendsen>`, and :doc:`fix langevin <fix_langevin>`.
+
+This compute subtracts out degrees-of-freedom due to fixes that
+constrain molecular motion, such as :doc:`fix shake <fix_shake>` and
+:doc:`fix rigid <fix_rigid>`.  This means the temperature of groups of
+atoms that include these constraints will be computed correctly.  If
+needed, the subtracted degrees-of-freedom can be altered using the
+*extra* option of the :doc:`compute_modify <compute_modify>` command.
+
+See the :doc:`Howto thermostat <Howto_thermostat>` doc page for a
+discussion of different ways to compute temperature and perform
+thermostatting.
+
+
+----------
+
+
+**Output info:**
+
+This compute calculates a global scalar (the temperature) and a global
+vector of length 7, which can be accessed by indices 1-7.
+The first 6 elements of the vector are the KE tensor,
+and the 7-th is the cosine-shaped velocity amplitude :math:`V`,
+which can be used to calculate the reciprocal viscosity, as shown in the example.
+These values can be used by any command that uses global scalar or
+vector values from a compute as input.
+See the :doc:`Howto output <Howto_output>` doc page for an overview of LAMMPS output options.
+
+The scalar value calculated by this compute is "intensive".  The
+first 6 elements of vector values are "extensive",
+and the 7-th element of vector values is "intensive".
+
+The scalar value will be in temperature :doc:`units <units>`.  The
+first 6 elements of vector values will be in energy :doc:`units <units>`.
+The 7-th element of vector value will be in velocity :doc:`units <units>`.
+
+Restrictions
+""""""""""""
+
+This command is only available when LAMMPS was built with the USER-VISCOSITY package.
+
+Related commands
+""""""""""""""""
+
+:doc:`fix accelerate/cos <fix_accelerate_cos>`
+
+Default
+"""""""
+ none
+
+
+----------
+
+.. _Hess:
+
+**(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209–217.

doc/src/fix.rst

@@ -165,6 +165,7 @@ The individual style names on the :doc:`Commands fix <Commands_fix>` doc
 page are followed by one or more of (g,i,k,o,t) to indicate which
 accelerated styles exist.
 
+* :doc:`accelerate/cos <fix_accelerate_cos>` - apply cosine-shaped acceleration to atoms
 * :doc:`adapt <fix_adapt>` - change a simulation parameter over time
 * :doc:`adapt/fep <fix_adapt_fep>` - enhanced version of fix adapt
 * :doc:`addforce <fix_addforce>` - add a force to each atom

doc/src/fix_accelerate_cos.rst

@@ -0,0 +1,106 @@
+.. index:: fix accelerate/cos
+
+fix accelerate/cos command
+==========================
+
+Syntax
+""""""
+
+
+.. parsed-literal::
+
+   fix ID group-ID accelerate value
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* accelerate/cos = style name of this fix command
+* value = amplitude of acceleration (in unit of force/mass)
+
+
+Examples
+""""""""
+
+
+.. parsed-literal::
+
+   fix 1 all accelerate/cos 0.02E-5
+
+Description
+"""""""""""
+
+Give each atom a acceleration in x-direciton based on its z coordinate.
+The acceleration is a periodic function along the z-direction:
+
+.. math::
+
+   a_{x}(z) = A \cos \left(\frac{2 \pi z}{l_{z}}\right)
+
+where :math:`A` is the acceleration amplitude, :math:`l_z` is the z-length
+of the simulation box.
+At steady state, the acceleration generates a velocity profile:
+
+.. math::
+
+   v_{x}(z) = V \cos \left(\frac{2 \pi z}{l_{z}}\right)
+
+The generated velocity amplitude :math:`V` is related to the
+shear viscosity :math:`\eta` by:
+
+.. math::
+
+   V = \frac{A \rho}{\eta}\left(\frac{l_{z}}{2 \pi}\right)^{2}
+
+
+and it can be obtained from ensemble average of the velocity profile:
+
+.. math::
+
+   V = \frac{\sum_i 2 m_{i} v_{i, x} \cos \left(\frac{2 \pi z_i}{l_{z}}\right)}{\sum_i m_{i}}
+
+where :math:`m_i`, :math:`v_{i,x}` and :math:`z_i` are the mass,
+x-component velocity and z coordinate of a particle.
+
+The velocity amplitude :math:`V` can be calculated with :doc:`compute viscosity/cos <compute_viscosity_cos>`,
+which enables viscosity calculation with periodic perturbation method,
+as described by :ref:`Hess<Hess>`.
+Because the applied acceleration drives the system away from equilibration,
+the calculated shear viscosity is lower than the intrinsic viscosity
+due to the shear-thinning effect.
+Extrapolation to zero acceleration should generally be performed to
+predict the zero-shear viscosity.
+As the shear stress decreases, the signal-noise ratio decreases rapidly,
+the simulation time must be extended accordingly to get converged result.
+
+In order to get meaningful result, the group ID of this fix should be all.
+
+
+----------
+
+
+**Restart, fix\_modify, output, run start/stop, minimize info:**
+
+No information about this fix is written to binary restart files.
+None of the fix_modify options are relevant to this fix.
+No global or per-atom quantities are stored by this fix for access by various output commands.
+No parameter of this fix can be used with the start/stop keywords of the run command.
+This fix is not invoked during energy minimization.
+
+Restrictions
+""""""""""""
+
+This command is only available when LAMMPS was built with the USER-VISCOSITY package.
+
+Related commands
+""""""""""""""""
+
+:doc:`compute viscosity/cos <compute_viscosity_cos>`
+
+Default
+"""""""
+ none
+
+
+----------
+
+.. _Hess:
+
+**(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209–217.

examples/USER/viscosity/README

@@ -0,0 +1,5 @@
+The input script in in.1000SPCE.lmp gives an example on how to calculate
+the shear viscosity of SPCE water with the periodic perturbation method.
+For a production run, much longer simulation is required.
+Note that the acceleration amplitude of `fix accelerate/cos` should be carefully chosen,
+and an extrapolation is generally required to get the viscosity at zero-shear.

src/USER-VISCOSITY/README

@@ -2,7 +2,7 @@ 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
+There is an 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.

src/USER-VISCOSITY/compute_viscosity_cos.cpp

@@ -38,7 +38,8 @@ ComputeViscosityCos::ComputeViscosityCos(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = vector_flag = 1;
   size_vector = 7;
   extscalar = 0;
-  extvector = 1;
+  extvector = -1;
+  extlist = new int[7]{1,1,1,1,1,1,0};
   tempflag = 1;
   tempbias = 1;