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)