docs/reference-manual/special/viscosity-calculation.rst

   1 Viscosity calculation
   2 ---------------------
   3
   4 The shear viscosity is a property of liquids that can be determined
   5 easily by experiment. It is useful for parameterizing a force field
   6 because it is a kinetic property, while most other properties which are
   7 used for parameterization are thermodynamic. The viscosity is also an
   8 important property, since it influences the rates of conformational
   9 changes of molecules solvated in the liquid.
  10
  11 The viscosity can be calculated from an equilibrium simulation using an
  12 Einstein relation:
  13
  14 .. math::  \eta = \frac{1}{2}\frac{V}{k_B T} \lim_{t \rightarrow \infty}
  15            \frac{\mbox{d}}{\mbox{d} t} \left\langle
  16            \left( \int_{t_0}^{{t_0}+t} P_{xz}(t') \mbox{d} t' \right)^2
  17            \right\rangle_{t_0}
  18            :label: eqneinsteinrelation
  19
  20 This can be done with :ref:`gmx energy <gmx energy>`. This method converges
  21 very slowly \ :ref:`149 <refHess2002a>`, and as such a nanosecond simulation might not
  22 be long enough for an accurate determination of the viscosity. The
  23 result is very dependent on the treatment of the electrostatics. Using a
  24 (short) cut-off results in large noise on the off-diagonal pressure
  25 elements, which can increase the calculated viscosity by an order of
  26 magnitude.
  27
  28 |Gromacs| also has a non-equilibrium method for determining the
  29 viscosity \ :ref:`149 <refHess2002a>`. This makes use of the fact that energy, which is
  30 fed into system by external forces, is dissipated through viscous
  31 friction. The generated heat is removed by coupling to a heat bath. For
  32 a Newtonian liquid adding a small force will result in a velocity
  33 gradient according to the following equation:
  34
  35 .. math:: a_x(z) + \frac{\eta}{\rho} \frac{\partial^2 v_x(z)}{\partial z^2} = 0
  36           :label: eqnviscositygradiant
  37
  38 Here we have applied an acceleration :math:`a_x(z)` in the
  39 :math:`x`-direction, which is a function of the :math:`z`-coordinate. In
  40 |Gromacs| the acceleration profile is:
  41
  42 .. math:: a_x(z) = A \cos\left(\frac{2\pi z}{l_z}\right)
  43           :label: eqnviscosityacceleration
  44
  45 where :math:`l_z` is the height of the box. The generated velocity
  46 profile is:
  47
  48 .. math:: v_x(z) = V \cos\left(\frac{2\pi z}{l_z}\right)
  49           :label: eqnviscosityprofile1
  50
  51 .. math:: V = A \frac{\rho}{\eta}\left(\frac{l_z}{2\pi}\right)^2
  52           :label: eqnviscosityprofile2
  53
  54 The viscosity can be calculated from :math:`A` and :math:`V`:
  55
  56 .. math:: \eta = \frac{A}{V}\rho \left(\frac{l_z}{2\pi}\right)^2
  57           :label: eqnvisc
  58
  59 In the simulation :math:`V` is defined as:
  60
  61 .. math:: V = \frac{\displaystyle \sum_{i=1}^N m_i v_{i,x} 2 \cos\left(\frac{2\pi z}{l_z}\right)}
  62           {\displaystyle \sum_{i=1}^N m_i}
  63           :label: eqnsimulationviscosity
  64
  65 The generated velocity profile is not coupled to the heat bath.
  66 Moreover, the velocity profile is excluded from the kinetic energy. One
  67 would like :math:`V` to be as large as possible to get good statistics.
  68 However, the shear rate should not be so high that the system gets too
  69 far from equilibrium. The maximum shear rate occurs where the cosine is
  70 zero, the rate being:
  71
  72 .. math:: \mbox{sh}_{\max} =  \max_z \left| \frac{\partial v_x(z)}{\partial z} \right|
  73           = A \frac{\rho}{\eta} \frac{l_z}{2\pi}
  74           :label: eqnshearrate
  75
  76 For a simulation with: :math:`\eta=10^{-3}`
  77 [kgm:math:`^{-1}`\ s\ :math:`^{-1}`],
  78 :math:`\rho=10^3`\ [kgm:math:`^{-3}`] and :math:`l_z=2\pi`\ [nm],
  79 :math:`\mbox{sh}_{\max}=1`\ [psnm:math:`^{-1}`] :math:`A`. This shear
  80 rate should be smaller than one over the longest correlation time in the
  81 system. For most liquids, this will be the rotation correlation time,
  82 which is around 10 ps. In this case, :math:`A` should be smaller than
  83 0.1[nmps\ :math:`^{-2}`]. When the shear rate is too high, the observed
  84 viscosity will be too low. Because :math:`V` is proportional to the
  85 square of the box height, the optimal box is elongated in the
  86 :math:`z`-direction. In general, a simulation length of 100 ps is enough
  87 to obtain an accurate value for the viscosity.
  88
  89 The heat generated by the viscous friction is removed by coupling to a
  90 heat bath. Because this coupling is not instantaneous the real
  91 temperature of the liquid will be slightly lower than the observed
  92 temperature. Berendsen derived this temperature
  93 shift \ :ref:`31 <refBerendsen91>`, which can be written in terms of the
  94 shear rate as:
  95
  96 .. math:: T_s = \frac{\eta\,\tau}{2 \rho\,C_v} \mbox{sh}_{\max}^2
  97           :label: eqnberendsentempshift
  98
  99 where :math:`\tau` is the coupling time for the Berendsen thermostat
 100 and :math:`C_v` is the heat capacity. Using the values of the example
 101 above, :math:`\tau=10^{-13}` [s] and :math:`C_v=2 \cdot 10^3`\ [J
 102 kg\ :math:`^{-1}`\ K\ :math:`^{-1}`], we get:
 103 :math:`T_s=25`\ [Kps:math:`^{-2}`]sh\ :math:`_{\max}^2`. When we want
 104 the shear rate to be smaller than :math:`1/10`\ [ps:math:`^{-1}`],
 105 :math:`T_s` is smaller than 0.25[K], which is negligible.
 106
 107 **Note** that the system has to build up the velocity profile when
 108 starting from an equilibrium state. This build-up time is of the order
 109 of the correlation time of the liquid.
 110
 111 Two quantities are written to the energy file, along with their averages
 112 and fluctuations: :math:`V` and :math:`1/\eta`, as obtained from
 113 (:eq:`%s <eqnvisc>`).
 114