lammps 模拟tip4p冰

【求助】lammps 模拟tip4p冰，密度总是调不对 （1）

想请各位帮忙看看我的输入文件哪里有问题。
压强等于0,得到的密度是0.933g/cm3左右。文献上是0.933g/cm3。总是不一样。输入的结构和周期性边界条件也反复检查过。xyz三个方向分别重复了5周期。

Radial distribution functions and densities for the SPC/E, TIP4P and TIP5P models for liquid water and ices Ih, Ic, II, III, IV, V, VI, VII, VIII, IX, XI and XIIw

种类1是O,2是H；L-J参数选用的是那篇文献的参数epsilong=0.006722  sigma=3.154000 

输入的in文件：
log                 03.ThermoMD.log
units               metal
atom_style       full

read_data        large.water

mass             1 15.9994
mass             2 1.008

group            g1 type 1
group            g2 type 2

pair_style       lj/cut/coul/long/tip4p 1 2 1 1 0.15 8.5 8.5
pair_coeff       1 1   0.006722  3.154000 8.5
pair_coeff       * 2    0.000000  0.000000 8.5

bond_style       harmonic
bond_coeff       1 527.20    0.9572
angle_style      harmonic
angle_coeff      1 37.95     104.52

kspace_style     pppm/tip4p 0.0001

neighbor        2 bin
neigh_modify    delay 0

timestep         0.002
thermo           100

thermo_style     custom step temp pe ke  etotal press vol v_varDensitySI pxx pyy pzz enthalpy epair evdwl ecoul  xlo xhi ylo yhi  zlo zhi
thermo_modify   flush yes line one

dump            water all custom 10000 04.Traj.part/04.Trajectory.*.xyi id type x y z
restart         100000 X.restart/X.restart

velocity        all create 250 143573
fix               fixShake all shake 0.0001 20 0 b 1 a 1
fix               T      all npt 250 250 0.2 xyz 0.0 0.0 2.0  drag  0.2
run             2000000
unfix           T
输入的坐标文件：

    3000    atoms
    2000    bonds
    1000    angles

    2       atom types
    1       bond types
    1       angle types

     -11.2250000      11.2250000       xlo xhi
     -19.4417000      19.4417000       ylo yhi
     -18.2950000      18.2950000       zlo zhi

Masses

       1    15.99940
       2     1.00800

Atoms

      1      1   1   -1.04   -8.209 -18.493 -17.992
      2      1   2    0.52   -8.209 -18.502 -17.035
      3      1   2    0.52   -8.209 -19.495 -18.320
      4      2   1   -1.04   -8.209 -18.417 -15.253
      5      2   2    0.52   -7.414 -18.022 -14.924
      6      2   2    0.52   -9.004 -18.022 -14.924
      7      3   1   -1.04   -8.209 -13.265 -14.328
      8      3   2    0.52   -8.209 -13.256 -13.371
      9      3   2    0.52   -8.209 -12.263 -14.656
     10      4   1   -1.04   -8.209 -13.341 -11.589

Response：

即使你使用完全相同的输入文件，使用同一个版本的程序来计算，但机器不一样，编译程序不一样，都会导致数值计算结果的一些很小的差别，这样的差别在MD计算中，随着运行步数的增加会累积和放大。假如输入文件不完全相同，比如说初始速度的分布采用不同的参数、原子位置输入时小数点位数不一样，温压控制参数不同，计算步数和统计平均的步数不一样，计算结果有差别就是更自然的事了。 另外我刚才看了一下，你的文献是用Monte Carlo方法模拟的，你是用分子动力学模拟的，方法就不同嘛，结果有差别，这再也正常不过了。文献列出的实验值是0.920，你得到0.933，他是0.937，你的结果还好过他，你该把这次模拟看成是成功的模拟才对。你纠结什么啊。

ASk：

非常感谢你细致的建议！另外还有个问题就是：  我对于tip4p/ew模型和tip4p模型的差别不是很清楚。tip4p/ew是不是用的截断的静电相互作用，而tip4p模型用的是长程经典相互作用呢࿱