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:

(409)\[\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}\]

This can be done with gmx energy. This method converges rather slowly 149, and usually hundreds of nanoseconds are needed 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. It is most convenient to use the Einstein relation. Because gmx mdrun stores averages of quantities computed every nstcalcenergy steps, the Einstein relation can use these averages and thus writing to energy file (nstenergy) can be done infrequently. This avoids the overhead of very large energy files that are needed with the autocorrelation function approach.

GROMACS also has two non-equilibrium methods for determining the viscosity. The recommended method is to apply a shear to the system, see the next section. The viscosity can the be measured as the stress of the corresponding off-diagonal element of the pressure tensor divided by the shear rate. This is a straightforward procedure. The only disadvantage is that one needs to balance the cost of long simulatiosn at low shear rate due to low signal to noise ratio to the risk of shear thinning appearing at higher shear rates. Running at multiple shear rates might be necessary to ensure that one is in the linear regime.

The second non-equilibrium method is called “cosine acceleration”. 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:

(410)\[a_x(z) + \frac{\eta}{\rho} \frac{\partial^2 v_x(z)}{\partial z^2} = 0\]

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

(411)\[a_x(z) = A \cos\left(\frac{2\pi z}{l_z}\right)\]

where \(l_z\) is the height of the box. The generated velocity profile is:

(412)\[v_x(z) = V \cos\left(\frac{2\pi z}{l_z}\right)\]

(413)\[V = A \frac{\rho}{\eta}\left(\frac{l_z}{2\pi}\right)^2\]

The viscosity can be calculated from \(A\) and \(V\):

(414)\[\eta = \frac{A}{V}\rho \left(\frac{l_z}{2\pi}\right)^2\]

In the simulation \(V\) is defined as:

(415)\[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}\]

The generated velocity profile is not coupled to the heat bath. Moreover, the velocity profile is excluded from the kinetic energy. One would like \(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:

(416)\[\mbox{sh}_{\max} = \max_z \left| \frac{\partial v_x(z)}{\partial z} \right| = A \frac{\rho}{\eta} \frac{l_z}{2\pi}\]

For a simulation with: \(\eta=10^{-3}\) [kgm\(^{-1}\)s\(^{-1}\)], \(\rho=10^3\)[kgm\(^{-3}\)] and \(l_z=2\pi\)[nm], \(\mbox{sh}_{\max}=1\)[psnm\(^{-1}\)] \(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, \(A\) should be smaller than 0.1[nmps\(^{-2}\)]. When the shear rate is too high, the observed viscosity will be too low. Because \(V\) is proportional to the square of the box height, the optimal box is elongated in the \(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 31, which can be written in terms of the shear rate as:

(417)\[T_s = \frac{\eta\,\tau}{2 \rho\,C_v} \mbox{sh}_{\max}^2\]

where \(\tau\) is the coupling time for the Berendsen thermostat and \(C_v\) is the heat capacity. Using the values of the example above, \(\tau=10^{-13}\) [s] and \(C_v=2 \cdot 10^3\)[J kg\(^{-1}\)K\(^{-1}\)], we get: \(T_s=25\)[Kps\(^{-2}\)]sh\(_{\max}^2\). When we want the shear rate to be smaller than \(1/10\)[ps\(^{-1}\)], \(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: \(V\) and \(1/\eta\), as obtained from ((414)).