calculating viscosity for LJ system via pressure tensor in 3D using Green-Kubo

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I am trying to calculate the viscosity for a LJ system in 3D using Green-Kubo but I am having trouble replicating a known example. I understand the example Green-Kubo calculation provided by LAMMPS is just an example, and I have been trying to modify it so as to give a known literature value for the viscosity. This involves the autocorrelation function of the pressure tensor component, p_{xy}. That is, I am trying to find a well behaved quantity

< p_{xy}(t_0) * p_{xy}(t_0 + t) >

I have been reading the documentation and the couple threads that appear on this topic, for example:

http://lammps.sandia.gov/tutorials/italy14/italy_kappa_viscosity_Mar14.pdf

I should be able to just dump the pressure tensor components and build this autocorrelation function post processing. But, I think something is fundamentally wrong with my simulations. No matter how frequently I dump the pressure tensor component “p_{xy}”, the autocorrelation function does not show a smooth decay like I would expect for an autocorrelation function. I have attached a pic of my autocorrelation function showing the independence of my thermo dump frequency.

The autocorrelation function simply decays immediately after the first step (I have normalized the autocorrelation function). The literature data I am comparing to is Erpenbecks “Shear viscosity of the LJ fluid near the triple point: Green-Kubo results.” I have included his fig5 data for the 4000 particle case. Erpenbeck uses a slightly different simulation, but I would expect similar results.

I have included my input below. I am doing NVE after an NVT equilibration (I pick up a previous restart file with a considerable equilibration period), which should be the proper ensemble/method for Green-Kubo. I make sure to skip the NVT thermo output when building my autocorrelation function.

Any thoughts or comments? Thank you in advanced,
–Joe

#sample LAMMPS input script for viscosity of liquid Ar

units lj
variable T equal 0.722
variable V equal vol
variable dt equal 0.001
variable rc equal 2.5
variable ThermoStep equal 2

variable x equal 10
variable y equal $x
variable z equal $x

setup problem

dimension 3
boundary p p p

read_restart EQUIL2.restart

pair_style lj/cut {rc} pair_coeff * * 1.0 1.0 pair_modify shift yes timestep {dt}
thermo ${ThermoStep}

variable theV equal vol
print “vol = ${theV}”

variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz

thermo_style custom step temp press v_pxy v_pxz v_pyz

equilibration and thermalization

velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 0.1
run 10000 # just to pick up

reset_timestep 0

unfix NVT

fix NVE all nve

run 10000

AutoCorr.eps (184 KB)

2

The example in examples/VISCOSITY/in.gk.2d is reproducing

known values accurately (see the README in the same dir).
As is examples/KAPPA/in.heatflux
for a 3d problem for thermal conductivity (which is conceptually
very similar).

Steve

3

I revisited the example and README for the 2d case as you pointed out. Do you recall what reference gives these known values? I do not see mention of one. The 2d example mentions a paper by Daivis and Evans, but this seems to be only pertinent for the theory of choosing components of the pressure tensor, as discussed in the appendix of that paper.

I have read the paper, Woodcock, AIChE Journal, 52, 438 (2006), you mention in your slides (http://lammps.sandia.gov/tutorials/italy14/italy_kappa_viscosity_Mar14.pdf) when comparing the various methods for calculating the viscosity (including Green-Kubo). But I don’t see the state point ρ* = 0.6, T* = 1.0, Rc = 2.5 σ anywhere in that paper.

My apologies in advanced if I missed a reference somewhere, but I have looked through the manual and the examples multiple times, and have followed up on some of the mentioned papers and I am not sure where to find the published reference for the state point mentioned.

Thanks very much,
–Joe

4

I did some digging but don’t have an ideal answer for you.

You’re right, the Woodcock paper doesn’t have that exact state point,
not sure why I cited that paper. It also seems to indicate that
the 3d state point rho=0.6, T=1.0 is 2-phase, but I don’t think
that is the case.

Two variables to worry about are the dimension (2d vs 3d) and the
cutoff. The Woodcock paper is 3d (as are all others I could find),
whereas the examples are 2d, mainly to enable simple pictures/movies
and quick calculations.

The examples are also using cut = 2.5 sigma, and I couldn’t find the value
in the Woodcock paper. A lot of papers use 3.0. I don’t think either
of these variables make a big difference when at a stable state point
(e.g. away from a triple point), but I’m not certain.

The safest thing for you to do is probably modify the scripts in
examples/VISCOSITY to run 3d problems for which there are literature
values.

Gary Grest assures me that in 3d for rho = 0.6, T = 1.0, cut = 3.0,
the correct viscosity = ~1.0. See the start of section B in the
attached paper. Note that it says density = 0.66 (for a 0 pressure
run). Also it is quoting kinematic viscosity = 1.53, but you have to
multiply by density to get shear viscosity.

Also note that while the 2d scripts in examples are quick/dirty runs,
they are all giving roughly the same answer for the 4 methods, so that
gives some confidence they are giving the right value for the 2d state
point.

Steve

JChemPhys_132_174106.pdf (767 KB)