diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 13336fa89b..8230a88986 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -161,5 +161,6 @@ KOKKOS, o = USER-OMP, t = OPT. * :doc:`torque/chunk ` * :doc:`vacf ` * :doc:`vcm/chunk ` + * :doc:`viscosity/cos ` * :doc:`voronoi/atom ` * :doc:`xrd ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index ce7cade0f6..67696912a8 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -22,6 +22,7 @@ OPT. .. table_from_list:: :columns: 5 + * :doc:`accelerate/cos ` * :doc:`adapt ` * :doc:`adapt/fep ` * :doc:`addforce ` diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index 8251c5301e..7c8a849f65 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -104,6 +104,7 @@ page gives those details. * :ref:`USER-UEF ` * :ref:`USER-VTK ` * :ref:`USER-YAFF ` + * :ref:`USER-VISCOSITY ` ---------- @@ -2365,3 +2366,29 @@ which discuss the `QuickFF `_ methodology. * :doc:`pair_style mm3/switch3/coulgauss/long ` * :doc:`pair_style 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 ` +* :doc:`compute viscosity/cos ` +* examples/USER/viscosity diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 0726a502bc..c09ba7fbeb 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -307,6 +307,7 @@ The individual style names on the :doc:`Commands compute ` doc * :doc:`torque/chunk ` - torque applied on each chunk * :doc:`vacf ` - velocity auto-correlation function of group of atoms * :doc:`vcm/chunk ` - velocity of center-of-mass for each chunk +* :doc:`viscosity/cos ` - velocity profile under cosine-shaped acceleration * :doc:`voronoi/atom ` - Voronoi volume and neighbors for each atom * :doc:`xrd ` - x-ray diffraction intensity on a mesh of reciprocal lattice nodes diff --git a/doc/src/compute_viscosity_cos.rst b/doc/src/compute_viscosity_cos.rst new file mode 100644 index 0000000000..8e7bb2b70e --- /dev/null +++ b/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 ` 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 `, :doc:`fix npt `, etc. + +This command together with :doc:`fix_accelerate/cos` +enables viscosity calculation with periodic perturbation method, +as described by :ref:`Hess`. +An acceleration along the x-direction is applied to the simulation system +by using :doc:`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 ` 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 `, :doc:`fix temp/rescale `, :doc:`fix temp/berendsen `, and :doc:`fix langevin `. + +This compute subtracts out degrees-of-freedom due to fixes that +constrain molecular motion, such as :doc:`fix shake ` and +:doc:`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 ` command. + +See the :doc:`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 ` 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 `. The +first 6 elements of vector values will be in energy :doc:`units `. +The 7-th element of vector value will be in velocity :doc:`units `. + +Restrictions +"""""""""""" + +This command is only available when LAMMPS was built with the USER-VISCOSITY package. + +Related commands +"""""""""""""""" + +:doc:`fix accelerate/cos ` + +Default +""""""" + none + + +---------- + +.. _Hess: + +**(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209–217. diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 2b5ed48ac9..cb2e4bb5b1 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -165,6 +165,7 @@ The individual style names on the :doc:`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 ` - apply cosine-shaped acceleration to atoms * :doc:`adapt ` - change a simulation parameter over time * :doc:`adapt/fep ` - enhanced version of fix adapt * :doc:`addforce ` - add a force to each atom diff --git a/doc/src/fix_accelerate_cos.rst b/doc/src/fix_accelerate_cos.rst new file mode 100644 index 0000000000..0a780ae59c --- /dev/null +++ b/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 ` 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 `, +which enables viscosity calculation with periodic perturbation method, +as described by :ref:`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 ` + +Default +""""""" + none + + +---------- + +.. _Hess: + +**(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209–217. diff --git a/examples/USER/viscosity/README b/examples/USER/viscosity/README new file mode 100644 index 0000000000..1221509ea7 --- /dev/null +++ b/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. diff --git a/src/USER-VISCOSITY/README b/src/USER-VISCOSITY/README index f24bf6c2b5..13c50ec286 100644 --- a/src/USER-VISCOSITY/README +++ b/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. diff --git a/src/USER-VISCOSITY/compute_viscosity_cos.cpp b/src/USER-VISCOSITY/compute_viscosity_cos.cpp index 44103b501f..6f69969759 100644 --- a/src/USER-VISCOSITY/compute_viscosity_cos.cpp +++ b/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;