I am currently utilizing the GK method to compute the viscosity of a polymer system comprising primarily benzene rings and hydrocarbon chains. The input script used is as follows:
units real
variable T equal 408.15 # run temperature
variable Tinit equal 408.15 # equilibration temperature
variable V equal vol
variable dt equal 1
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}
# setup problem
dimension 3
boundary p p p
atom_style full
neighbor 2.0 bin
neigh_modify delay 5 check yes
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style cvff
pair_style lj/cut/coul/long 10
kspace_style pppm 1.0e-4
read_data polymer.data
timestep ${dt}
thermo $d
# equilibration and thermalization
velocity all create ${Tinit} 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp ${Tinit} ${Tinit} 10 drag 0.2
run 8000
# viscosity calculation, switch to NVE if desired
velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 10 drag 0.2
unfix NVT
fix NVE all nve
reset_timestep 0
variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz
fix SS all ave/correlate $s $p $d &
v_pxy v_pxz v_pyz type auto file S0St.dat ave running
variable scale equal ${convert}/(${kB}*$T)*$V*$s*${dt}
variable v11 equal trap(f_SS[3])*${scale}
variable v22 equal trap(f_SS[4])*${scale}
variable v33 equal trap(f_SS[5])*${scale}
variable v equal (v_v11+v_v22+v_v33)/3.0
thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33 v_v
run 1000000
variable ndens equal count(all)/vol
print "average viscosity: $v [Pa.s] @ $T K, ${ndens} atoms/A^3"
However, the outcome I received is:
average viscosity: 0.00180735846760119 [Pa.s] @ 408.15 K, 0.104814548328098 atoms/A^3
approximately 1000 times lower than the actual value ranging from 1-3 pa.s at the same temperature. I am uncertain whether the reason behind this discrepancy is due to the size of the system or simulation settings or whether the GK method is unsuitable for such a polymer system. For reference, my model size is 4040 40 Angstrom. NVT and NVE ensembles give similar results.