Gromacs LINCS 键角旋转超过30度

chuxuezhe · 发表于 Post on 2018-3-31 11:37:41

各位老师好：
我需要模拟水溶液中花色苷与儿茶素形成复合物的过程。通过chemdraw建模，antechamber+acpype.py建立拓扑结构，原子类型添加到相应的ffnonbond中。能量最小化后，进行动力学模拟。力场分别试了amber99sb-ildn和gromos54a7，但是在模拟开始，均出现了：

relative constraint deviation after LINCS:
rms 0.001054, max 0.006588 (between atoms 1770 and 1772)
bonds that rotated more than 30 degrees:

曾尝试：
缩小了模拟步长，mdp文件中约束h-angles，增加盒子大小，使用不同的力场（amber和gromos），采用不同方法能量最小化，都没有解决问题。（附件是动力学输入文件，各个分子的结构和拓扑文件，以及LINCS提示的键角出问题的原子（都是儿茶素中的，与氢原子相连的，以及连接苯环的键））。

谢谢各位老师帮助。

Novice · 发表于 Post on 2025-7-16 17:26:29

本帖最后由 Novice 于 2025-7-16 18:07 编辑

我在使用packmol创建的盒子模拟苯溶液时，也遇到了类似的情况。

我具体情况是这样的：
GMX版本为2024.5，溶剂分子top文件来自www.virtualchemistry.org提供的GAFF-ESP-2018包，自己研究的有机芳香分子top用sobtop产生。

采用virtualchemistry的二氯甲烷top文件直接模拟GAFF-ESP-2012包中的二氯液体盒子时（我看了填充密度和实际接近），能顺利进行npt模拟。
如果按照实际密度推算出应填充的分子数目填充盒子时，二氯甲烷溶液和苯溶液（填充有两个其它有机分子）的em过程均没问题，但是npt模拟时几乎立刻就出现了如帖子http://bbs.keinsci.com/thread-19759-1-1.html中的Too many LINCS warnings 的错误。
二氯甲烷体系我将填充分子数目减少到实际的1/3时，发现模拟能顺利进行了。即使这样，模拟时体积还是增加了15%左右（考虑到二氯甲烷沸点低，特意将模拟温度设为为288K）。
不过我发现对于“正常”结束的体系，直接载入模拟后的.gro文件，发现分子似乎散了（有些原子未成键，测量发现未成键的C-Cl键长1.82A，实际成键的约1.79A），是virtualchemistry的二氯甲烷top不好么？

到了苯体系，我像二氯甲烷体系一样先尝试了按实际密度填充，发现npt模拟时出现Too many LINCS warnings的错误。
采用苯填充分子数目减少到实际的1/6体系模拟时，进行大概500ps后，提示体系偏离平衡太多而退出的错误（我记得是这样，具体情况忘了），有意思的是我发现此时盒子内部出现了苯的团聚液滴(其它为气象）的现象。

我试了调整不同的压浴设置和模拟步长，发现均不能进行NPT模拟；我尝试了关闭压浴（NVT），发现能顺利进行模拟。

实在找不到思路了，我就单独模拟纯苯的溶液体系，由于GAFF-ESP-2012包中没有苯盒子文件，所以我使用packmol自行创建了一个10nm^3，填充4200个苯分子（按密度计算实际能填>6000个）的盒子。em顺利结束后进行npt模拟时，发现仍然不能通过，报错信息就如楼主说的Gromacs LINCS 键角旋转超过30度。

我就想是不是virtualchemistry的苯top文件有问题，于是用sobtop基于Hessian产生了苯的top文件，对比后发现电荷、键长参数基本相似；含氢的键角参数virtualchemistry的要比sobtop产生的要大2倍；dihedrals项参数由于使用的是不同函数类型所以我没法比较。但是我发现使用sobtop基于Hessian产生的top也不能进行苯液体的npt模拟，而且我觉得virtualchemistry的苯top应该也没问题。

找到这篇帖子后，我试了采用1fs小步长,切换berendsen和C-rescale压浴，发现均不能顺利进行模拟（Gromacs LINCS 键角旋转超过30度报错）。而关闭约束后，更是直接报错
step 400, will finish Wed Jul 16 16:46:30 2025 Segmentation fault (core dumped)
out文件也没有更多信息。

模拟中使用的相关文件如下

md.zip (1.35 MB, 下载次数 Times of downloads: 0)

1.为什么溶解有两个分子的二氯甲烷体系需要设置的比实际密度低很多才能进行NPT模拟?
2.我应该如何进行苯溶液体系的NPT模拟呢？

请有经验的老师给指个明路。

liusy401 · 发表于 Post on 2024-1-28 19:20:14

请问去掉约束是constraint_algorithm = lincs， constraints = h-bonds 这两个约束都去掉吗

社会主义小战士 · 发表于 Post on 2023-4-28 14:59:19

sobereva 发表于 2023-4-28 12:29
1fs跑一阵后，如果2fs续跑时能正常跑，很长的模拟都用2fs

谢谢社长！

sobereva · 发表于 Post on 2023-4-28 12:29:11

社会主义小战士发表于 2023-4-27 13:12
社长好！
这里去掉约束，用1fs步长的话，只是用来做预平衡么？后期跑很长的模拟也要用1FS么？

1fs跑一阵后，如果2fs续跑时能正常跑，很长的模拟都用2fs

社会主义小战士 · 发表于 Post on 2023-4-27 13:12:31

sobereva 发表于 2018-3-31 13:45
把约束去掉再试，1fs步长时候不需要做约束
并且通过退火设定，升温时从0K慢慢升温
注意grompp时候的任何 ...

社长好！
这里去掉约束，用1fs步长的话，只是用来做预平衡么？后期跑很长的模拟也要用1FS么？

chuxuezhe · 发表于 Post on 2022-11-4 16:56:40

ZHIMING123 发表于 2022-11-4 11:14
你好，您是如何解决的，我也是运行nvt上来就报错这个，电脑还有弹窗提示出来。

你好，去掉约束，时间设置为1fs就行。不知道体系搭建情况，必要时可以设置一个较低的温度先平衡下体系。

ZHIMING123 · 发表于 Post on 2022-11-4 11:14:36

chuxuezhe 发表于 2018-4-2 06:22
谢谢Sob老师，问题已经解决了。

你好，您是如何解决的，我也是运行nvt上来就报错这个，电脑还有弹窗提示出来。

sobereva · 发表于 Post on 2022-9-21 22:39:16

JCenter 发表于 2022-9-20 12:01
好的，谢谢sob老师。那我的产生相控压方法就用berendsen了。使用此方法对于后续的MSD和数量密度分析，RDF ...

没问题

JCenter · 发表于 Post on 2022-9-20 12:01:23

sobereva 发表于 2022-9-20 05:12
用berendsen就完了
或者尝试gmx >=2021版本里的C-rescale压浴

好的，谢谢sob老师。那我的产生相控压方法就用berendsen了。使用此方法对于后续的MSD和数量密度分析，RDF分析有影响吗

sobereva · 发表于 Post on 2022-9-20 05:12:04

JCenter 发表于 2022-9-20 02:00
好的，谢谢老师。在您回复前，我今天尝试过了。MD_NPT中使用berendsen是可以正常跑完的，盒子也没问题。 ...

用berendsen就完了
或者尝试gmx >=2021版本里的C-rescale压浴

JCenter · 发表于 Post on 2022-9-20 02:00:14

本帖最后由 JCenter 于 2022-9-20 02:17 编辑

sobereva 发表于 2022-9-19 23:43
继续使用heat_NPT时用的Berendsen压浴看情况

好的，谢谢老师。在您回复前，我今天尝试过了。MD_NPT中使用berendsen是可以正常跑完的，盒子也没问题。嗯~，目前的情况自己表达的可能有点乱，所以我重新叙述一遍存在的问题和我尝试解决的过程，也方便其他和我一样的新手看下我尝试的方法，问题是：对于1个5*5*5nm且含有580个正己烷分子的盒子进行MD模拟，在em.mdp(emtol为50)和heat_NPT1ns(0-200ps升温0-298k,后800ps保持温度，berendsenNPT) 正常跑完，检查能量最小化后体系的势能收敛，检查heat_NPT后体系的温度，压力和密度均达到平衡状态。在heat_NPT的最后一帧的基础上进行MD_NPT2ns(控压方法pr)续跑，一段时间后，出现崩溃，系统提示我多条bonds that rotated more than 30 degrees。后面我按照sob老师的http://sobereva.com/soft/Sobtop/#FAQ8一文与老师发过的模拟崩溃的原因ppt图片进行错误的排除。
1.检查结构文件，拓扑文件与坐标文件的原子顺序一致。因为是用sob老师的sobtop程序创建的top文件，另外，其创建的gro文件转成pdb用packmol创建盒子。另外，分子也是在b3lyp-d3(bj)/6-311g**级别下优化过的。
2.检查拓扑文件，文件内的[ bonds ]、[ angles ]、[ dihedrals ]、[ pairs ]字段对照结构看了一下，没问题，另一方面也是sobtop程序产生的，gaffl力场，比较放心。而且我也对照过antechamber+acpype.py产生的拓扑文件，相一致。再就是原子电荷，我是按照老师的resp2(0.5)电荷文章进行求算代入的。这方面我按照老师的文章一步一步做的，没理由相信原子电荷是错误的。但也有可能错误，为此我上传了文件。3.并行计算方式，在使用-ntmpi 1后依旧失败，这方面的原因排除。
4.最后检查mdp文件。我已经使用了能量最小化消除了不合理的接触。heat_NPT1ns的结果很平衡。推测问题出MD_NPT.mdp上，我先将之前的2fs步长+h-bonds约束改为1fs步长+none，结果出现lincs报警和停止。接着压控部分，我在heat_NPT.mdp中将berendsen更换为PR,尽管这样做不合理，这样是为了验证方法的效果，最后在跑heat_NPT后系统崩溃，也是出现了bonds that rotated more than 30 degrees。反过来，我在MD_NPT.mdp中将PR更换为berendsen，结果正常跑完。我也尝试过在pr控压方法的基础上对tau_p,从2.0改为3.0.结果虽然正常跑完。但是分子间变得极其稀疏，盒子边长从5nm 变成150nm。这个现象目前还没有找到答案。就目前来说，在MD_NPT中压控采用berendsen方法可以正常结束。但是，不确定在产生相中使用berendsen会不会被审稿人质疑(http://bbs.keinsci.com/thread-12837-1-1.html)。另外，对于tau_p=3.0这方面，不知道是够有突破口。
以上是我的问题排除的过程，不知道下一步该怎么办了。老师，您有什么建议吗。

sobereva · 发表于 Post on 2022-9-19 23:43:01

JCenter 发表于 2022-9-19 19:22
老师，我用了gmx grompp -f MD_NPT.mdp -c hexane_heat_NPT.gro -p hexane.top -t hexane_heat_NPT.cpt - ...

继续使用heat_NPT时用的Berendsen压浴看情况

JCenter · 发表于 Post on 2022-9-19 19:22:05

JCenter 发表于 2022-9-19 12:03
好的，谢谢sob老师。

老师，我用了gmx grompp -f MD_NPT.mdp -c hexane_heat_NPT.gro -p hexane.top -t hexane_heat_NPT.cpt -o hexane_MD.tpr 在heat_NPT(1ns,没问题)基础上续跑MD_NPT(2ns)。但跑了一段时间结果还是出现bonds that rotated more than 30 degrees的问题。

JCenter · 发表于 Post on 2022-9-19 12:03:34

sobereva 发表于 2022-9-19 03:10
heat_NPT能正常跑就没理由之后的动力学没法正常跑。续跑时候记得延续之前动力学最后的状态

好的，谢谢sob老师。

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[GROMACS] Gromacs LINCS 键角旋转超过30度

浏览过的版块