Viscosity calculation
---------------------

The shear viscosity is a property of liquids that can be determined
easily by experiment. It is useful for parameterizing a force field
because it is a kinetic property, while most other properties which are
used for parameterization are thermodynamic. The viscosity is also an
important property, since it influences the rates of conformational
changes of molecules solvated in the liquid.

The viscosity can be calculated from an equilibrium simulation using an
Einstein relation:

.. math::  \eta = \frac{1}{2}\frac{V}{k_B T} \lim_{t \rightarrow \infty}
           \frac{\mbox{d}}{\mbox{d} t} \left\langle 
           \left( \int_{t_0}^{{t_0}+t} P_{xz}(t') \mbox{d} t' \right)^2
           \right\rangle_{t_0}
           :label: eqneinsteinrelation

This can be done with :ref:`gmx energy <gmx energy>`. This method converges
very slowly \ :ref:`149 <refHess2002a>`, and as such a nanosecond simulation might not
be long enough for an accurate determination of the viscosity. The
result is very dependent on the treatment of the electrostatics. Using a
(short) cut-off results in large noise on the off-diagonal pressure
elements, which can increase the calculated viscosity by an order of
magnitude.

|Gromacs| also has a non-equilibrium method for determining the
viscosity \ :ref:`149 <refHess2002a>`. This makes use of the fact that energy, which is
fed into system by external forces, is dissipated through viscous
friction. The generated heat is removed by coupling to a heat bath. For
a Newtonian liquid adding a small force will result in a velocity
gradient according to the following equation:

.. math:: a_x(z) + \frac{\eta}{\rho} \frac{\partial^2 v_x(z)}{\partial z^2} = 0
          :label: eqnviscositygradiant

Here we have applied an acceleration :math:`a_x(z)` in the
:math:`x`-direction, which is a function of the :math:`z`-coordinate. In
|Gromacs| the acceleration profile is:

.. math:: a_x(z) = A \cos\left(\frac{2\pi z}{l_z}\right)
          :label: eqnviscosityacceleration

where :math:`l_z` is the height of the box. The generated velocity
profile is:

.. math:: v_x(z) = V \cos\left(\frac{2\pi z}{l_z}\right)
          :label: eqnviscosityprofile1

.. math:: V = A \frac{\rho}{\eta}\left(\frac{l_z}{2\pi}\right)^2
          :label: eqnviscosityprofile2

The viscosity can be calculated from :math:`A` and :math:`V`:

.. math:: \eta = \frac{A}{V}\rho \left(\frac{l_z}{2\pi}\right)^2
          :label: eqnvisc

In the simulation :math:`V` is defined as:

.. math:: V = \frac{\displaystyle \sum_{i=1}^N m_i v_{i,x} 2 \cos\left(\frac{2\pi z}{l_z}\right)}
          {\displaystyle \sum_{i=1}^N m_i}
          :label: eqnsimulationviscosity

The generated velocity profile is not coupled to the heat bath.
Moreover, the velocity profile is excluded from the kinetic energy. One
would like :math:`V` to be as large as possible to get good statistics.
However, the shear rate should not be so high that the system gets too
far from equilibrium. The maximum shear rate occurs where the cosine is
zero, the rate being:

.. math:: \mbox{sh}_{\max} =  \max_z \left| \frac{\partial v_x(z)}{\partial z} \right|
          = A \frac{\rho}{\eta} \frac{l_z}{2\pi}
          :label: eqnshearrate

For a simulation with: :math:`\eta=10^{-3}`
[kgm\ :math:`^{-1}`\ s\ :math:`^{-1}`],
:math:`\rho=10^3`\ [kgm\ :math:`^{-3}`] and :math:`l_z=2\pi`\ [nm],
:math:`\mbox{sh}_{\max}=1`\ [psnm\ :math:`^{-1}`] :math:`A`. This shear
rate should be smaller than one over the longest correlation time in the
system. For most liquids, this will be the rotation correlation time,
which is around 10 ps. In this case, :math:`A` should be smaller than
0.1[nmps\ :math:`^{-2}`]. When the shear rate is too high, the observed
viscosity will be too low. Because :math:`V` is proportional to the
square of the box height, the optimal box is elongated in the
:math:`z`-direction. In general, a simulation length of 100 ps is enough
to obtain an accurate value for the viscosity.

The heat generated by the viscous friction is removed by coupling to a
heat bath. Because this coupling is not instantaneous the real
temperature of the liquid will be slightly lower than the observed
temperature. Berendsen derived this temperature
shift \ :ref:`31 <refBerendsen91>`, which can be written in terms of the
shear rate as:

.. math:: T_s = \frac{\eta\,\tau}{2 \rho\,C_v} \mbox{sh}_{\max}^2
          :label: eqnberendsentempshift

where :math:`\tau` is the coupling time for the Berendsen thermostat
and :math:`C_v` is the heat capacity. Using the values of the example
above, :math:`\tau=10^{-13}` [s] and :math:`C_v=2 \cdot 10^3`\ [J
kg\ :math:`^{-1}`\ K\ :math:`^{-1}`], we get:
:math:`T_s=25`\ [Kps\ :math:`^{-2}`]sh\ :math:`_{\max}^2`. When we want
the shear rate to be smaller than :math:`1/10`\ [ps\ :math:`^{-1}`],
:math:`T_s` is smaller than 0.25[K], which is negligible.

**Note** that the system has to build up the velocity profile when
starting from an equilibrium state. This build-up time is of the order
of the correlation time of the liquid.

Two quantities are written to the energy file, along with their averages
and fluctuations: :math:`V` and :math:`1/\eta`, as obtained from
(:eq:`%s <eqnvisc>`).