关于Green-Kubo方法模拟水的粘度

关于Green-Kubo方法模拟水的粘度

求问有没有大神用lammps做过水的模拟？
我现在改编了一下manual里面计算液氩粘度的那个算例，来模拟水的粘度，300k下得到的结果是0.2mPa.s，查得的标准值是0.8mPa.s，请路过的大神帮忙看一下是哪里出了问题，多谢！

in文件贴在下面，data文件和log文件在附件中。

# LAMMPS input script for viscosity of purewater

units real
variable T equal 300
variable V equal vol
variable dt equal 0.4
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
bond_style harmonic
angle_style harmonic

pair_style lj/cut/tip4p/long 1 2 1 1 0.125 12 15
kspace_style pppm/tip4p 1e-6
read_data sketch.data

pair_coeff 1 2 0.0 0.0 # O-H
pair_coeff 2 2 0.0 0.0 # H-H
pair_coeff 1 1 0.16275 3.16435 # O-O

timestep ${dt}
thermo $d

# equilibration and thermalization

velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 40 drag 0.2
run 100000

# viscosity calculation, switch to NVE if desired

#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}
thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
run 200000
variable v equal (v_v11+v_v22+v_v33)/3.0
variable ndens equal count(all)/vol
print "average viscosity: $v [Pa.s/ @ $T K, ${ndens} /A^3"

哥们，应该不需要把real转换成SI 这一步，即不要单位转换的因子convert
不要光看最后一个结果，要仔细分析跑动粘度，即粘度随关联时间的变化情况。
通常，跑动粘度需要一定的关联时间才能收敛。
另外，一次模拟也不够，因为算粘度误差会很大，每次得到的结果也会差别较大。

水分子选的模型会影响剪切粘度：
TIP3P TIP4P TIP5P SPC/E TIP4P/2005 Expt.a
0.321 0.494 0.699 0.729 0.855 0.896

一楼的同学说不用转化为SI ，这个我保留意见。我在自己的模拟中使用了单位转化。关于你的data文件和in文件我说一下我的看法：你使用的是tip4p模型，却使用了很多的spce模型的参数而不是tip4p模型的参数。
这是我自己的模拟结果：（300K,1.29E+5mol/m3）
spce tip4p tip4p2005
0.567 0.420 0.721
使用的参数为lammps手册“6.8 TIP4P water mode”中给出的参数，l