# [lammps-users] Muller-Plate calculation of Viscosity

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)/(2tLxLy), where Lx = Ly = 6 sigma, and t = 5000006.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

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.

If you look in the MP paper (cited on the LAMMPS doc page), there is a factor
of 2 in this calculation, which you may be missing.

Steve