Hello,

I am trying to use the Muller-Plathe method to calculate viscosity in LAMMPS (version 9Jan09) using the same procedure as described in their paper. The first test is on Ar, but I can not calculate the value of the viscosity to be 3.2 (at T=0.722 and rho = 0.849; all in LJ units). Instead I get a viscosity around 6 in LJ units.

When I run LAMMPS (input script below), I get a value of 14907.6 as a final result from the viscosity fix (f_3) and implement it into the following formula j= (f_3)/(2*t*Lx*Ly), where Lx = Ly = 6 sigma, and t = 500000*6.965E-03 (number of steps times dt in LJ units). This results in j = 0.059454.

I also take the velocity profile from the end of the run reported by fix 2. I plot the x-direction velocity as a function of the z-height and take the slope of the linear fit. This is then dvx/dz, which has a value of 0.0098.

To calculate the viscosity I just divide j by dvx/dz and the result is consistently around 6. I have tried changing the frequency of momentum transfer and that does not affect the calculated viscosity. Does anything stand out as being wrong at this point? Any help is greatly appreciated!

Thanks,

Ken

#viscosity calculation of argon (T=0.722, rho =0.894)

boundary p p p

units lj

atom_style atomic

lattice fcc 0.849

region box block 0 6 0 6 0 18

create_box 1 box

create_atoms 1 box

mass 1 1.0

velocity all create 0.722 8700 loop geom

pair_style lj/cut 3.0

pair_modify tail yes

pair_coeff 1 1 1.0 1.0 3.0

neighbor 0.3 bin

neigh_modify delay 0 every 20 check no

fix 1 all nvt 0.722 0.722 0.464

thermo_style one

thermo 100

timestep 6.965e-03

run 100000

write_restart restart.create

unfix 1

neighbor 0.3 bin

neigh_modify delay 0 every 20 check no

fix 1 all nve

fix 2 all ave/spatial 100 1000 100000 z lower 0.05 vx file visc.profile units reduced

fix 3 all viscosity 100 x z 20 swap 1 vtarget 1.5

thermo_style custom step etotal ke evdwl etail temp f_3

thermo 100

timestep 6.965e-3

run 500000