I want to get viscosity using NEMD water shear viscosity,below is my script .The water model is SPC/E.I am getting 0.0013 Pas.s .Do any one has any idea why the value I am getting is not the same as literature (0.0004Pa.s) at shear rate 0.5ps-1.Any help would be highly appreciated.!Is there anything wrong in my script or any wrong with conversion?
variable T equal 300.0 # run temperature
variable V equal vol
variable dt equal 1.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}
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
remove hybrid if not necessary
pair_style hybrid lj/cut/coul/long 12.0 12.0
pair_modify mix geometric tail yes
read_data eq.lmp
pair_coeff 1 1 lj/cut/coul/long 0.000000 0.000000 # Hw Hw
pair_coeff 2 2 lj/cut/coul/long 0.155425 3.165500 # Ow Ow
fix SHAKE all shake 0.0001 20 0 b 1
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
change_box all triclinic
kspace_style pppm 1.0e-5
fix 11 all nve
fix 21 all langevin $T $T 10.0 768012
restart 10000 nve.restart
timestep 1.0
thermo 1000
#thermo_style custom step temp pe ke etotal vol press density
run 100000
unfix 11
reset_timestep 0
velocity all scale T
variable srate0 equal 5.0e13 # shear rate [1/s]
variable srate equal {srate0}/1e15 # convert shear rate unit from s^-1 to fs^-1
variable xyrate equal ${srate}/ly # Shear rate per unit length (fs^-1/Å)
fix 1 all nvt/sllod temp $T T 10.0 fix 2 all deform 1 xy erate {xyrate} remap v # Apply deformation
compute layers all chunk/atom bin/1d y center 0.05 units reduced
fix 3 all ave/chunk 20 250 5000 layers vx file profile.nemd.3d-${srate}
compute temp after deforming box
compute usual all temp
compute tilt all temp/deform
restart 100000 shear.restart
thermo 1000
thermo_style custom step temp c_usual pe etotal press vol density pxy
thermo_modify temp tilt
dump TRAJ all custom 10000 dump.lammpstrj id mol type element x y z
dump_modify TRAJ element H O sort id
run 100000
Convert viscosity to Pa·s
variable scale equal {convert}/({kB}$T)$V*s*{dt}
variable visc equal -pxy/v_xyrate*{scale}
fix vave all ave/time 10 100 1000 v_visc ave running file visc-avg-{srate}
thermo 5000
thermo_style custom step temp c_usual pe etotal press vol density pxy v_visc f_vave
thermo_modify temp tilt
run 500000
print “Average Viscosity in Pa·s: ${visc}”