【求助】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模型用的是长程经典相互作用呢