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模型用的是长程经典相互作用呢？这个文献上面是tip4p模型，大概和我的in文件也不同吧。因为我用的应该是tip4p/ew模型。

TIP4P/Ice (2005) force field for simulation of ice 1h

From: Ankit Mishra <ankitmis@...114...> Date: Fri, 12 May 2017 01:47:49 -0700

Dear lammps Users,

I am trying to simulate ice 1h structure which is hexagonal in nature using TIP4P/Ice (2005) force field. My objective is to obtain phonon density of states of water at T = 100 K (though temperature does not effect DOS, but still I am doing at low temperature to avoid diffusion effects) using a force field which properly describes this phase. And compare it with the behaviour of water at similar T & P using a reaxFF implementation for a system, which also has water in it. This way I can compare how good the reaxFF implementation is.

Since ice 1h is the most stable phase at my required condition of temperature and pressure, hence I am trying to obtain its phonon DOS.

I have constructed the input script based on the information related to charges, and LJ parameters on the lammps web site (attached with this email).

But after a trying a lot of things, I am still unable to obtain a stable structure at even 10K.

So here I list all the problems I encountered and solutions I tried:

1. I encountered the first problem in kspace_style, since TIP4P/Ice require kspace_style pppm/tip4p for coloumb interaction calculations, and this style only works with orthogonal strcutures.