Lab: viscosity by hand

The setup

A shear flow (layers sliding past each other) is resisted by a shear stress, the off-diagonal pressure P_xy. At rest its average is zero, but it fluctuates. The time these fluctuations take to die away is the viscosity:

System: the 3D LJ liquid of Lab 6.9 (256 atoms, ρ = 0.8, T = 1.0, cutoff 1.6, V = 320), 100 000 steps with a gentle Bussi thermostat. In an isotropic liquid five combinations of the pressure tensor carry the same information: P_xy, P_xz, P_yz, (P_xx − P_yy)/2 and (P_yy − P_zz)/2.

Step 1

The shear stress by hand

Why: P_xy measures how much x-momentum crosses a plane perpendicular to y, per unit area and time. Atoms carry momentum across by moving (kinetic part) and push across by forces (virial part). Both are needed, and both come from quantities the simulation already has.

What the numbers say: One configuration gives a small number of either sign. In the LJ liquid P_xy averages to zero (−0.02 ± 0.27 per step, the mean is pure noise) and fluctuates by ±0.27, larger than any shear stress a flow would create in a simulation. The viscosity is not in the average but in how long a fluctuation lasts. Note that the diagonal pressure of Lab 3.8 is the same formula with x and x instead of x and y.

Step 2

Green–Kubo by hand

Why: The stress autocorrelation C(t) = ⟨P_xy(0)P_xy(t)⟩ starts at the size of the fluctuations and decays as the liquid relaxes. The longer it takes, the more viscous the liquid. Its integral, times V/k_BT, is η: the same recipe as D from the VACF in Lab 6.7.

What the numbers say: Both routes are a few multiplications. The factor V is the important difference from diffusion: the viscosity is a property of the whole box, one number per time origin, like the conductivity in Lab 6.7, so it is a collective quantity and as noisy as one. With only ten numbers the two routes disagree (0.367 vs 0.300); the question of where to stop the integral is already visible, since C turns negative at lag 4.

Step 3

The LJ liquid: where to stop the integral

Why: Now the real data: 100 000 steps, all five stress components averaged. The running integral η(t) should rise and level off once C(t) has decayed. In practice the noise in C(t) keeps adding to the integral, so the plateau has to be chosen.

What the numbers say: Most of the viscosity comes from the first 0.2 time units, where the stress correlation decays; by t = 1–2 the integral sits at 2.40–2.46, and both routes agree on about 2.42. Beyond that the integral drifts (2.86 at t = 4, 3.66 at t = 8), not because the liquid has a long memory but because the noisy tail of C(t) is integrated too: the random walk of the noise wins. Practical rules: plot η(t), choose the plateau where C(t) has decayed to the noise, and take the error bar from independent blocks or runs (here ±0.12). In viscous liquids such as ionic liquids the plateau comes much later, and fitting C(t) or using many independent runs becomes necessary.

Step 4

Cross-checks: Yeh–Hummer and Stokes–Einstein

Why: A viscosity should be consistent with diffusion. Two checks: the box-size dependence of D in Lab 6.9 must follow the Yeh–Hummer slope ξk_BT/(6πη), and Stokes–Einstein, D = k_BT/(cπηR), should give a sensible particle radius R.

What the numbers say: Two independent routes agree: the viscosity from stress fluctuations (2.42) predicts exactly how D changes with box size (the three larger boxes collapse onto 0.0688), and the slope of D vs 1/L gives back 2.43. Stokes–Einstein, a formula for a sphere in a continuum, works surprisingly well for atoms if the "slip" boundary condition is used: it gives a diameter of 0.96σ, almost exactly the LJ size. With "stick" the particle comes out too small. Limitation: Stokes–Einstein is a continuum law; it fails for small solutes in large solvents, in supercooled liquids, and for ions with strong local structure, which is exactly where deviations from it are interesting.

Step 5

Convergence: why viscosity needs long runs

Why: Viscosity, like conductivity, is a property of the whole box, so it does not get more samples from more atoms. How long must the run be for a given precision, and how much do the five equivalent stress components help?

What the numbers say: A run of 20 time units, plenty for D of this liquid, gives η only to ±30 %; it takes about 400 time units (100 000 steps) to reach ±5 %, and the error falls only as one over the square root of the run length. Using all five shear components is free and is worth five times more simulation. Larger boxes do not help much: the fluctuations shrink as 1/V while the V in front of the integral grows, so the relative noise is about the same. For viscous liquids, where the stress correlation lasts 100–1000 times longer, people use many independent runs, the "time-decomposition" method, or nonequilibrium shear (NEMD) instead.

Step 6

Units: from LJ to mPa s

Why: Like every reduced quantity, the viscosity becomes physical once we pick a substance. For argon the LJ unit of viscosity is √(mε)/σ².

What the numbers say: A liquid of LJ argon at these conditions comes out at about 0.22 mPa s, a quarter of water's viscosity, a sensible value for a simple liquid of small atoms (with our shortened cutoff it is a model, not a precise prediction for argon). In real-unit codes the Green–Kubo prefactor V/k_BT and the pressure units are the usual sources of factor-of-10⁵ mistakes; checking a code's viscosity on the LJ liquid against a published value is a good first test. LAMMPS: fix ave/correlate on the pressure tensor; GROMACS: gmx energy -vis.

Summary