计算化学公社

 找回密码 Forget password
 注册 Register
Views: 15317|回复 Reply: 22
打印 Print 上一主题 Last thread 下一主题 Next thread

[GROMACS] 使用packmol建立冰水两相,低温高压下进行NPT,冰却融化了

[复制链接 Copy URL]

1479

帖子

0

威望

4541

eV
积分
6020

Level 6 (一方通行)

本帖最后由 牧生 于 2021-2-5 21:02 编辑

对培训班中的冰的模拟有了兴趣
但自己搭建了冰水两相,低温高压下进行NPT,希望观察到水结冰过程,但结果却观察到冰融化,请问一下哪里不对

具体操作如下:
以2304.gro为基础,建立了冰水模型,如图(ice为1,水分子 400)



编写top文件

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"               
[ system ]
ice
[ molecules ]
; Compound        nmols
SOL                    1168           


再进行能量最小化
gmx grompp -f em.mdp -c mix.gro -p H2O.top -o ice.tpr -maxwarn 1
gmx mdrun -deffnm ice -v

然后直接进行npt模拟,



integrator  = md      
nsteps      = 1000000     
dt          = 0.001     
; Output control
nstxout     = 50     
nstvout     = 50      
nstenergy   = 50     
nstlog      = 50      
; Bond parameters
continuation            = no   
constraint_algorithm    = lincs   
constraints             = h-bonds
lincs_iter              = 1      
lincs_order             = 4   
; Neighborsearching
cutoff-scheme   = Verlet
ns_type         = grid     
nstlist         = 10      
rcoulomb        = 1.0      
rvdw            = 1.0     
; Electrostatics
coulombtype     = PME      
pme_order       = 4      
fourierspacing  = 0.16     
; Temperature coupling is on
tcoupl      = V-rescale         
tc-grps     =system   
tau_t       = 0.1               
ref_t       = 240             温度设定240 K
; Pressure coupling is on
pcoupl              = Parrinello-Rahman   
pcoupltype          = isotropic         
tau_p               = 2.0              
ref_p               = 200.0                 压力设定200 bar
compressibility     = 4.5e-5            
refcoord_scaling    = com
; Periodic boundary conditions
pbc     = xyz     
DispCorr    = EnerPres  
gen_vel     = no     


VMD观察,发现冰融化了,截图如下

播放轨迹动画,盒子明显收缩,但仍有真空区,表明水分子不够多,但这个应该不会影响结冰吧?

同时有另一个问题,为什么还有真空区的存在?








很快,冰就全部融化了





     
使用gmx energy -f npt.edr -o npt.xvg  看温度和压力的变化,温度倒是比较稳定,但是压力变化很大







又菜又爱玩

108

帖子

0

威望

3826

eV
积分
3934

Level 5 (御坂)

2#
发表于 Post on 2021-2-6 10:13:18 | 只看该作者 Only view this author
npt下肯定收缩

1479

帖子

0

威望

4541

eV
积分
6020

Level 6 (一方通行)

3#
 楼主 Author| 发表于 Post on 2021-2-6 10:14:43 | 只看该作者 Only view this author

npt以后还是有真空区。这是为何呢
又菜又爱玩

108

帖子

0

威望

3826

eV
积分
3934

Level 5 (御坂)

4#
发表于 Post on 2021-2-6 10:19:28 | 只看该作者 Only view this author
牧生 发表于 2021-2-6 10:14
npt以后还是有真空区。这是为何呢

收缩过程中肯定会有,你最后一张图片不就没有了嘛

1479

帖子

0

威望

4541

eV
积分
6020

Level 6 (一方通行)

5#
 楼主 Author| 发表于 Post on 2021-2-6 11:14:08 | 只看该作者 Only view this author
本帖最后由 牧生 于 2021-2-6 12:30 编辑
Kangtor 发表于 2021-2-6 10:19
收缩过程中肯定会有,你最后一张图片不就没有了嘛

在播放轨迹过程中,一开始水盒子就迅速收缩,但一直保持有真空区。过一会,冰融化了,水就到处跑了。然后真空区就看不到了

我已经设置了低温,但为何还融化了。请帮忙看看问题所在。
又菜又爱玩

553

帖子

0

威望

3324

eV
积分
3877

Level 5 (御坂)

6#
发表于 Post on 2021-2-6 12:03:00 | 只看该作者 Only view this author
那个真空就很没道理。

水的实验密度你也知道,笔算一下liquid phase要放多少水分子,加个正常的数目再跑,就好多了。

6万

帖子

99

威望

5万

eV
积分
120130

管理员

公社社长

7#
发表于 Post on 2021-2-6 20:49:43 | 只看该作者 Only view this author
SPCE水模型的冰点偏离真实冰点很大,要对冰感兴趣应当用我的gromacs培训班课上讲的TIP4P-ice来模拟
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办极高质量的各种计算化学类培训:初级量子化学培训班中级量子化学培训班高级量子化学培训班量子化学波函数分析与Multiwfn程序培训班分子动力学与GROMACS培训班CP2K第一性原理计算培训班,内容介绍以及往届资料购买请点击相应链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的高速路!培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入北京科音微信公众号获取北京科音培训的最新消息,并避免错过网上有价值的计算化学文章!
欢迎加入人气极高、专业性特别强的理论与计算化学综合交流群思想家公社QQ群(群号见此链接),合计达一万多人。北京科音培训班的学员在群中可申请VIP头衔,提问将得到群主Sobereva的最优先解答。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(十分强大、极为流行的量子化学波函数分析程序)
Google Scholar:https://scholar.google.com/citations?user=tiKE0qkAAAAJ
ResearchGate:https://www.researchgate.net/profile/Tian_Lu

1479

帖子

0

威望

4541

eV
积分
6020

Level 6 (一方通行)

8#
 楼主 Author| 发表于 Post on 2021-2-6 21:35:42 | 只看该作者 Only view this author
本帖最后由 牧生 于 2021-11-10 13:07 编辑
sobereva 发表于 2021-2-6 20:49
SPCE水模型的冰点偏离真实冰点很大,要对冰感兴趣应当用我的gromacs培训班课上讲的TIP4P-ice来模拟

谢谢。目前就卡在三点水转四点水的这一步。在centos下运行的命令,理应转化后,有n个水分子,就会多出来n个虚原子,就要改gro文件中的数目,加上n个。但似乎我运行命令后,并没有变化。正在找原因。
2021.11.10 给出运动命令后没变化的原因,
cat em-sol.gro |awk '{print $0;if($2=="OW"){a=$1;b=$3;c=$4;d=$5;e=$6;}if($2=="HW2")printf("%8s     MW %4d%8.3f%8.3f%8.3f\n",a,b,c,d,e)}' > newicebox.gro

这个命令,需要em-sol中的水是HW, OW等,如果是H, O,就不行了。
又菜又爱玩

236

帖子

0

威望

5064

eV
积分
5300

Level 6 (一方通行)

9#
发表于 Post on 2021-2-6 23:35:03 | 只看该作者 Only view this author
液体部分体积太大了,固液共存的胞应该是这样:

屏幕截图 2021-02-06 233458.png (237.79 KB, 下载次数 Times of downloads: 56)

屏幕截图 2021-02-06 233458.png

6

帖子

0

威望

77

eV
积分
83

Level 2 能力者

10#
发表于 Post on 2023-2-10 12:47:13 | 只看该作者 Only view this author
牧生 发表于 2021-2-6 21:35
谢谢。目前就卡在三点水转四点水的这一步。在centos下运行的命令,理应转化后,有n个水分子,就会多出来n ...

我也遇到了这个问题 请问现在解决了吗 老师

365

帖子

5

威望

3866

eV
积分
4331

Level 6 (一方通行)

Nerv

11#
发表于 Post on 2023-2-10 13:02:55 | 只看该作者 Only view this author
hzh12138 发表于 2023-2-10 12:47
我也遇到了这个问题 请问现在解决了吗 老师

有必要做三点水转四点水吗?tip4p水模型下冰的结构可以直接通过genice程序获得,先对纯冰体系跑个npt平衡一下,然后把盒子沿某一方向扩大一点,跑个nvt自然就得到冰水混合的体系了。
God's in his heaven,all is right with the world

1479

帖子

0

威望

4541

eV
积分
6020

Level 6 (一方通行)

12#
 楼主 Author| 发表于 Post on 2023-2-10 13:07:07 | 只看该作者 Only view this author
hzh12138 发表于 2023-2-10 12:47
我也遇到了这个问题 请问现在解决了吗 老师

解决了的。

http://bbs.keinsci.com/thread-9958-1-1.html     第七楼
又菜又爱玩

6

帖子

0

威望

77

eV
积分
83

Level 2 能力者

13#
发表于 Post on 2023-2-11 20:28:42 | 只看该作者 Only view this author
Lacrimosa 发表于 2023-2-10 13:02
有必要做三点水转四点水吗?tip4p水模型下冰的结构可以直接通过genice程序获得,先对纯冰体系跑个npt平衡 ...

谢谢您的提醒,这个想法对我这个初学者来说很有意思,有时间我也把这部分做做看,genice暂时还看不太懂

6

帖子

0

威望

77

eV
积分
83

Level 2 能力者

14#
发表于 Post on 2023-2-11 20:33:29 | 只看该作者 Only view this author
牧生 发表于 2023-2-10 13:07
解决了的。

http://bbs.keinsci.com/thread-9958-1-1.html     第七楼

谢谢。那个脚本很管用,我转化好后又重新做了一下,200 bar冰还是融化了,后来我调到450 bar就好了很多,可能初始体系不太理想吧,初始有3500个水分子和3072个冰分子.

1479

帖子

0

威望

4541

eV
积分
6020

Level 6 (一方通行)

15#
 楼主 Author| 发表于 Post on 2023-2-11 21:34:10 | 只看该作者 Only view this author
本帖最后由 牧生 于 2023-2-11 21:36 编辑
hzh12138 发表于 2023-2-11 20:33
谢谢。那个脚本很管用,我转化好后又重新做了一下,200 bar冰还是融化了,后来我调到450 bar就好了很多, ...

用tip4pice水模型,只要温度低于0度,压力即使是1bar,也会结冰的。如果冰化了,考虑自己用的水模型问题
genice可以直接通过pip安装,一行命令就解决。
又菜又爱玩

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2025-8-14 23:05 , Processed in 0.276907 second(s), 23 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list