计算化学公社

标题: 采用GMX对甲苯液体NPT模拟时盒子体积震荡 [打印本页]

作者
Author: Novice 时间: 2025-7-17 15:42
标题: 采用GMX对甲苯液体NPT模拟时盒子体积震荡
本帖最后由 Novice 于 2025-7-18 13:58 编辑

我使用packmol按照实际密度创建了一个包含1000个甲苯分子的、体积为5.62nm^3的盒子，

在使用GMX-2024.5模拟苯液体时，发现模拟过程中总是提示：
Warning: pressure scaling more than 1%, mu: 1.0154。
（注：甲苯分子top文件来自OPLS-AA的资源库https://github.com/leelasd/OPLSAA-DB。）

模拟时我是这么进行的：
1.先用cg算法进行em，收敛能量值至<100.
2. 采用如下设定进行升温NPT预平衡
dt       = 0.001  ; ps
nsteps    = 200000 ; 200ps
;
annealing          = single
annealing_npoints    = 3
annealing_time       = 0 50 100
annealing_temp       = 0 150 298.15
pbc = xyz
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb    = 1.0
vdwtype    = cut-off
rvdw       = 1.0
DispCorr    = EnerPres
;
Tcoupl  = V-rescale
tau_t = 0.2
tc_grps = system
ref_t = 298.15
;
Pcoupl    = berendsen
pcoupltype = isotropic
tau_p = 1.0
ref_p = 1.0
compressibility = 4.5e-5
;
gen_vel  = no
gen_temp = 273.15
gen_seed = -1

检查结果发现，在温度稳定区域，体系体积周期性波动，压力也会周期性偶然出现剧烈变化。

问了ChatGPT后，其建议我在以上两步之间加一个NVT升温步骤。于是，我按照如下过程进行了模拟：
1. 先用cg算法进行em，收敛能量值至<100.
2. 采用如下设定进行升温NVT模拟
dt       = 0.002  ; ps
nsteps    = 100000 ; 200ps
;
annealing          = single
annealing_npoints    = 3
annealing_time       = 0 50 100
annealing_temp       = 0 150 298.15
pbc = xyz
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb    = 1.0
vdwtype    = cut-off
rvdw       = 1.0
DispCorr    = EnerPres
;
Tcoupl  = V-rescale
tau_t = 0.2
tc_grps = system
ref_t = 298.15
;
gen_vel  = no
gen_temp = 273.15
gen_seed = -1
constraints = hbonds
3. 采用上面两步法中NPT的设定（但是去掉升温设定）进行预平衡。

但是发现无论是两步法还是3步法，平衡时NPT得到的体系体积周期性波动，压力也会周期性偶然出现剧烈变化（如下图所示）。

图1. 3步法时NVT过程中温度能达到平衡。

(, 下载次数 Times of downloads: 0)

图2. NPT时（不论是2步法还是3步法）盒子体积和压力震荡情况。
(, 下载次数 Times of downloads: 0)

请问我该如何做才能避免发生体系体积和压力的周期性波动，使体系真正平衡呢？

作者
Author: sobereva 时间: 2025-7-18 07:06
尝试c-rescale压浴，以及减小nstpcouple到10再试

作者
Author: Novice 时间: 2025-7-18 10:13
本帖最后由 Novice 于 2025-7-18 10:48 编辑

谢谢社长指教。
在增加nstpcouple=10的设定后，无论采用berendsen（tau_p = 1.0）还是采用c-rescale压浴（tau_p = 0.5）均能获得平衡的体系状态。

模拟过程中也没有发生任何类似：
Warning: pressure scaling more than 1%, mu: 1.0154
更没有类似如下错误退出的情况：
relative constraint deviation after LINCS:
……
bonds that rotated more than 30 degrees.
以及我也遇到过的：
WARNING: Listed nonbonded interaction between particles 9806 and 9807at distance 2.650 which is larger than the table limit 2.155 nm.

不过我有如下疑惑：为何两种压浴设定的能量值差别挺大的(如下图所示)?
(, 下载次数 Times of downloads: 0)

作者
Author: Novice 时间: 2025-7-18 10:46
本帖最后由 Novice 于 2025-7-18 14:09 编辑

仔细斟酌后，将我这次模拟的经历补充说明如下。

结合帖子http://bbs.keinsci.com/thread-9578-1-1.html里楼主模拟正己烷遇到的问题(该贴已浏览2万余次，说明很多人遇到类似问题)，以及我模拟二氯甲烷和甲苯体系遇到的问题（该贴26楼），我觉得
对于有机溶剂体系而言，使用控压模拟（并遇到各种难以解决的奇怪问题）时，增加nstpcouple=10设置十分关键。

我尝试模拟甲苯、二氯甲烷体系好几天了，前面的具体情况见下面链接26楼：
http://bbs.keinsci.com/thread-9578-2-1.html （Gromacs LINCS 键角旋转超过30度）

c-rescale压浴是我之前的模拟的首选，但是我发现c-rescale压浴采用tau_p = 0.5的常规设置时，直接就报Too many LINCS warnings（类似帖子http://bbs.keinsci.com/thread-19759-1-1.html）. 我尝试了减小模拟步长为1fs，去掉constraints，均未能解决问题。

在上面求助无果后，我自己尝试对c-rescale压浴修改了tau_p = 0.2，发现对于甲苯稀疏填充体系（为实际填充密度1/5）的情况，能进行完NPT模拟，模拟完后分析发现甲苯盒子体积和压力出现如本贴所说的周期性震荡，实际观察发现甲苯形成了液相聚集区和气相区分离情况。

于是我就尝试问了ChatGPT，它提示tau_p较小的确容易出现震荡，建议我用berendsen压浴和长tau_p（1.0）进行预平衡。我发现这个设置的确有一些效果，之前c-rescale压浴一模拟就崩溃的甲苯按实际密度填充的盒子，也能进行了，只不过震荡问题仍然无法解决。

社长在此处回复我后，我在对控压增加了nstpcouple=10设置，顺利解决问题！

————————————————————————
PS: 友情提示，谨慎采取GPT咨询专业问题获取答案！

比如它给我的berendsen压浴和长tau_p（1.0）进行预平衡的建议，虽然表面上让我能完成模拟，但是根本上的震荡问题并没有解决，也就是并没有获得平衡的体系！

而在社长建议我将nstpcouple设为10时，我不知道这个选项的意义，就问了它，它竟然告诉我：

从 GROMACS 4.5 起就废弃了 nstpcouple；
现在压力耦合每一步都自动进行更新，nstpcouple 不再需要手动设置；
GROMACS 会在新版中提示：
nstpcouple is obsolete and has no effect.

而还建议我，控制压力耦合“缓慢”或“温和”的方法是：使用 tau_p 控制耦合强度，compressibility 控制可压缩性（防止剧烈体积变动）。
我继续问了它甲苯溶液的compressibility应设置为多少，它建议我设为1.2e-4，并给出了常见集中溶剂的compressibility，也基本上在1.xe-4作用（比水4.5e-5大了近2.7倍，我没有求证）。

我没有怀疑社长的nstpcouple建议，但是我承认我病急乱投医了，我同时调整了nstpcouple和compressibility设置后，发现运行十分平稳。
后面我在nstpcouple=10下，再次尝试了仍将compressibility设为4.5e-5，以及采用之前极易崩溃的c-rescale压浴进行模拟，发现仍然运行平稳。

可见，强如ChatGPT，在咨询具体专业问题的解决方案时，大概率也是一本正经的胡说八道。

作者
Author: student0618 时间: 2025-7-18 11:27
本帖最后由 student0618 于 2025-7-18 11:29 编辑

可见，强如ChatGPT，在咨询具体专业问题的解决方案时，大概率也是一本正经的胡说八道。

用AI前建议先看这帖的三楼 http://bbs.keinsci.com/forum.php ... 29534&fromuid=64740

尤其是有Knowledge cutoff又不给资料来源的model如chatgpt更要小心。

作者
Author: Novice 时间: 2025-7-18 13:30
本帖最后由 Novice 于 2025-7-18 16:08 编辑

再次补充总结：

对于有机溶剂体系而言，使用控压模拟（并遇到各种难以解决的奇怪问题）时，增加nstpcouple=10设置十分关键。

若还是有问题，em后先进行NVT预平衡，然后再以1fs小步长进行NPT模拟，应该就行了。

我记得参加培训时社长说过，NVT预平衡是不必要的，经过我的这个例子我对此观点持保留意见。我发现即使经过比较严格的em（emtol = 60.0），如果直接用em结果跑升温NPT预平衡的话，很容易就出现Warning: pressure scaling more than 1%，甚至bonds that rotated more than 30 degrees退出的情况。而先进行NVT预平衡，然后再进行NPT模拟，模拟时就不会有任何问题。

-------------------------------------
此外说下我的另一个发现：
遇到问题时，切换一下力场可能也有帮助。
比如我进行甲苯模拟时，发现采用GAFF力场的甲苯盒子，即使能量最小化，以及NVT预平衡后，以2fs步长进行NPT模拟还是会频繁报
Warning: pressure scaling more than 1%，
最后甚至出现了LINCS WARNING：bonds that rotated more than 30 degrees而退出的的情况。
此时我将步长设为1fs，就没有出现任何问题了。

同样是甲苯体系，采用OPLS-AA力场，能量最小化，以及NVT预平衡后，以2fs步长进行NPT模拟虽然也频繁报
Warning: pressure scaling more than 1%，
但是体系能正常运行结束，且检查盒子压力、温度、体积也都是平衡的。
相信OPLS-AA力场若采用1fs步长模拟时，体系应该更稳定。

PS: 两次采用的盒子是分别采用packmol根据OPLS-AA DB和virtualchemistry.org提供的甲苯力场文件中的gro转化为pdb填充而来（我没仔细对比原子顺序是否一致，怕麻烦因此生成了两次盒子），初始构象并不完全相同，因此可能对以上结论有影响。如果不合理也欢迎大佬们指正。

作者
Author: sobereva 时间: 2025-7-18 18:58

如果直接用em结果跑升温NPT预平衡的话，很容易就出现Warning: pressure scaling more than 1%，甚至bonds that rotated more than 30 degrees退出的情况

关键不在于是否非得做NVT。NPT过程若用很大的tau_p，哪怕初始状态偏离平衡显著，也可以让盒子不容易产生较快的变化而减少出现此情况。只不过这样的tau_p对于一般目的的模拟往往就不适合了。另外，模拟初期的Warning: pressure scaling more than 1%很多情况本身就是无害的，等体系跑得较为平衡后自然就不容易再出现了。

此外，我的培训里也并非无条件地抵制做NVT，专门说了适合使用NVT的条件：
(, 下载次数 Times of downloads: 0)

关于nstpcouple，gromacs的开发者（尤其是Berk Hess）一味地追求计算效率，在较新版本引入了过激的默认设置，所以我建议你改回老版本的默认值
(, 下载次数 Times of downloads: 0)

作者
Author: Novice 时间: 2025-7-18 21:29
本帖最后由 Novice 于 2025-7-18 22:12 编辑

sobereva 发表于 2025-7-18 18:58
关键不在于是否非得做NVT。NPT过程若用很大的tau_p，哪怕初始状态偏离平衡显著，也可以让盒子不容易产生 ...

6楼关于NVT不是必须的说法我后来看您讲义上的内容，翻到了这一页。讲义真的是面面俱到，毫无瑕疵，太强了！想改说法时发现在后台审核了，所以就没改。

我刚尝试了在em后直接进行升温NPT并对c-rescale用较大的tau_p（1.5）和2fs步长，能不发生任何错误或警告而平稳完成预平衡。
倒是之前用较小的tau_p，虽然也能运行结束，但是导致体系盒子体积震荡。
我现在很好奇为啥常规值的tau_p会导致体系出现问题？

另外，能不能麻烦您解释下3楼我关于不同压浴设定为什么对体系能量影响这么大的疑问。

作者
Author: sobereva 时间: 2025-7-25 15:09

Novice 发表于 2025-7-18 21:29
6楼关于NVT不是必须的说法我后来看您讲义上的内容，翻到了这一页。讲义真的是面面俱到，毫无瑕疵，太强了 ...

较新版本默认的nstpcouple偏大，导致默认的tau_p下也容易有那个问题，gromacs论坛上也早有人遇到过

先看看温度曲线有什么差异。温度正比于原子平均动能

欢迎光临计算化学公社 (http://bbs.keinsci.com/)