\(\renewcommand{\AA}{\text{Å}}\)

# 10.3.7. Calculate viscosity

The shear viscosity \(\eta\) of a fluid can be measured in at least 6 ways using various options in LAMMPS. See the `examples/VISCOSITY` directory for scripts that implement the 5 methods discussed here for a simple Lennard-Jones fluid model and 1 method for SPC/E water model. Also, see the page on calculating thermal conductivity for an analogous discussion for thermal conductivity.

\(\eta\) is a measure of the propensity of a fluid to transmit momentum in a direction perpendicular to the direction of velocity or momentum flow. Alternatively it is the resistance the fluid has to being sheared. It is given by

\[J = -\eta \cdot \text{grad}(V_{\text{stream}})\]

where \(J\) is the momentum flux in units of momentum per area per time. and \(\text{grad}(V_{\text{stream}})\) is the spatial gradient of the velocity of the fluid moving in another direction, normal to the area through which the momentum flows. Viscosity thus has units of pressure-time.

The first method is to perform a non-equilibrium MD (NEMD) simulation by shearing the simulation box via the fix deform command, and using the fix nvt/sllod command to thermostat the fluid via the SLLOD equations of motion. Alternatively, as a second method, one or more moving walls can be used to shear the fluid in between them, again with some kind of thermostat that modifies only the thermal (non-shearing) components of velocity to prevent the fluid from heating up.

Note

A recent (2017) book by (Daivis and Todd) discusses use of the SLLOD method and non-equilibrium MD (NEMD) thermostatting generally, for both simple and complex fluids, e.g. molecular systems. The latter can be tricky to do correctly.

In both cases, the velocity profile setup in the fluid by this procedure can be monitored by the fix ave/chunk command, which determines \(\text{grad}(V_{\text{stream}})\) in the equation above. E.g. the derivative in the y-direction of the \(V_x\) component of fluid motion or \(\text{grad}(V_{\text{stream}}) = \frac{\text{d} V_x}{\text{d} y}\). The \(P_{xy}\) off-diagonal component of the pressure or stress tensor, as calculated by the compute pressure command, can also be monitored, which is the \(J\) term in the equation above. See the Howto nemd page for details on NEMD simulations.

The third method is to perform a reverse non-equilibrium MD simulation using the fix viscosity command which implements the rNEMD algorithm of Muller-Plathe. Momentum in one dimension is swapped between atoms in two different layers of the simulation box in a different dimension. This induces a velocity gradient which can be monitored with the fix ave/chunk command. The fix tallies the cumulative momentum transfer that it performs. See the fix viscosity command for details.

The fourth method is based on the Green-Kubo (GK) formula which relates the ensemble average of the auto-correlation of the stress/pressure tensor to \(\eta\). This can be done in a fully equilibrated simulation which is in contrast to the two preceding non-equilibrium methods, where momentum flows continuously through the simulation box.

Here is an example input script that calculates the viscosity of liquid Ar via the GK formalism:

    # Sample LAMMPS input script for viscosity of liquid Ar

    units       real
    variable    T equal 200.0       # run temperature
    variable    Tinit equal 250.0   # equilibration temperature
    variable    V equal vol
    variable    dt equal 4.0
    variable    p equal 400     # correlation length
    variable    s equal 5       # sample interval
    variable    d equal $p*$s   # dump interval

    # convert from LAMMPS real units to SI

    variable    kB equal 1.3806504e-23    # [J/K] Boltzmann
    variable    atm2Pa equal 101325.0
    variable    A2m equal 1.0e-10
    variable    fs2s equal 1.0e-15
    variable    convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}

    # setup problem

    dimension    3
    boundary     p p p
    lattice      fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
    region       box block 0 4 0 4 0 4
    create_box   1 box
    create_atoms 1 box
    mass         1 39.948
    pair_style   lj/cut 13.0
    pair_coeff   * * 0.2381 3.405
    timestep     ${dt}
    thermo       $d

    # equilibration and thermalization

    velocity     all create ${Tinit} 102486 mom yes rot yes dist gaussian
    fix          NVT all nvt temp ${Tinit} ${Tinit} 10 drag 0.2
    run          8000

    # viscosity calculation, switch to NVE if desired

    velocity     all create $T 102486 mom yes rot yes dist gaussian
    fix          NVT all nvt temp $T $T 10 drag 0.2
    #unfix       NVT
    #fix         NVE all nve

    reset_timestep 0
    variable     pxy equal pxy
    variable     pxz equal pxz
    variable     pyz equal pyz
    fix          SS all ave/correlate $s $p $d &
                 v_pxy v_pxz v_pyz type auto file S0St.dat ave running
    variable     scale equal ${convert}/(${kB}*$T)*$V*$s*${dt}
    variable     v11 equal trap(f_SS[3])*${scale}
    variable     v22 equal trap(f_SS[4])*${scale}
    variable     v33 equal trap(f_SS[5])*${scale}
    thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
    run          100000
    variable     v equal (v_v11+v_v22+v_v33)/3.0
    variable     ndens equal count(all)/vol
    print        "average viscosity: $v [Pa.s] @ $T K, ${ndens} atoms/A^3"

The fifth method is related to the above Green-Kubo method, but uses the Einstein formulation, analogous to the Einstein mean-square-displacement formulation for self-diffusivity. The time-integrated momentum fluxes play the role of Cartesian coordinates, whose mean-square displacement increases linearly with time at sufficiently long times.

The sixth is the periodic perturbation method, which is also a non-equilibrium MD method. However, instead of measuring the momentum flux in response to an applied velocity gradient, it measures the velocity profile in response to applied stress. A cosine-shaped periodic acceleration is added to the system via the fix accelerate/cos command, and the compute viscosity/cos command is used to monitor the generated velocity profile and remove the velocity bias before thermostatting.

Note

An article by (Hess) discussed the accuracy and efficiency of these methods.

(Daivis and Todd) Daivis and Todd, Nonequilibrium Molecular Dynamics (book), Cambridge University Press, https://doi.org/10.1017/9781139017848, (2017).

(Hess) Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209-217.