计算化学公社

标题: 求助：氨基化合物在水相中MD模拟奔溃 [打印本页]

作者
Author: nku_linz 时间: 2023-9-28 16:48
标题: 求助：氨基化合物在水相中MD模拟奔溃
本帖最后由 nku_linz 于 2023-9-28 21:37 编辑

采用gromacs模拟氨基化合物与皂苷的相互作用

两个化合物通过ORCA的B97-3c进行优化和振动分析，参考sobtop程序主页http://sobereva.com/soft/Sobtop/中的例1生成拓扑文件

接着于ORCA进行单点计算，并参考博文http://sobereva.com/441，通过Multiwfn计算RESP电荷（施加了等价性约束），把电荷信息换为RESP电荷

构建5*5*5的盒子放入一个CS（氨基化合物）和一个GA（皂苷），加入水分子并用Cl离子平衡盒子电荷

md进行能量最小化 emtol  = 10，最小化过程顺利，提示为
Polak-Ribiere Conjugate Gradients converged to machine precision in 2125 steps, but did not reach the requested Fmax < 10.
Potential Energy  = -2.4436003e+05
Maximum force    =  5.8923340e+02 on atom 176
Norm of force    =  6.9977082e+00

接着跑平衡相模拟，步长为2fs，步数为50000，使用Berendsen压浴，tau_p=0.5, 控温使用V-rescale，tau_t=0.2, 升温时长为30 ps，从0升到298.15K

平衡跑几步就报错，提示：
step 14580: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates

排查崩溃原因，把平衡相步长缩短至0.5fs，步数为200000，去掉约束，升温速率减慢，升温时长为80 ps，从0升到298.15K

平衡相可以跑更长时间，但仍然提示错误：
Step 178540  Warning: pressure scaling more than 1%, mu: 0.98725 0.98725 0.98725
step 178546: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.

看到压力有变化，调整tau_p=0.8，步数增加到250000，升温时长改为100ps

模拟到后面还是崩溃：
step 198716: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates

单独对两个分子做模拟，发现GA分子单独在水中跑模拟可以顺利进行，怀疑是氨基化合物CS的问题，gmx产生了模拟轨迹，调出崩溃前的一帧到VMD中查看

发现CS中61和13号原子，25号和80号原子大幅度偏移
(, 下载次数 Times of downloads: 8)

检查结构文件和拓扑文件，没看出哪里不妥，请教各位老师，这是哪里操作不当导致？[atomtypes]项统一放到ffnonbonded.itp文件中

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

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

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

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

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

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

作者
Author: sobereva 时间: 2023-9-29 05:52
CS直接用sobtop创建基于GAFF力场的拓扑文件就完了，根本不需要也不应该基于Hessian产生谐振势形式的力常数
仔细看http://sobereva.com/soft/Sobtop#FAQ11
而且网页里还都粗体字特意强调了

注：绝对不要以为什么时候都得像本例这样提供含有Hessian的文件并让Sobtop计算力常数。如果你要对普通有机体系产生拓扑文件，应该效仿的是例2直接用GAFF的参数，那是最省事的。老有人只看了本例，居然连例2都不看就去胡搞。也记得把FAQ10和FAQ11看了，搞清楚什么时候才有必要基于Hessian算力常数。

作者
Author: zdworld 时间: 2023-10-3 12:02

sobereva 发表于 2023-9-29 05:52
CS直接用sobtop创建基于GAFF力场的拓扑文件就完了，根本不需要也不应该基于Hessian产生谐振势形式的力常数
...

您好，我直接用sobtop的gaff力场参数化丙酮但是密度到了810 理论应该792，这种程度的偏差合理吗？如果基于DFT拟合键参数也差不多应该是非键参数的问题。

作者
Author: 喵星大佬 时间: 2023-10-3 12:13

zdworld 发表于 2023-10-3 12:02
您好，我直接用sobtop的gaff力场参数化丙酮但是密度到了810 理论应该792，这种程度的偏差合理吗？如果基 ...

合理，这误差已经不大了，就2.2%你还想咋样
gaff的原子类型考虑的化合物类型很多，vdW参数这点误差已经很小了
成键参数对密度影响很小，就是非键的问题，最多你再调调电荷

作者
Author: sobereva 时间: 2023-10-3 14:17

zdworld 发表于 2023-10-3 12:02
您好，我直接用sobtop的gaff力场参数化丙酮但是密度到了810 理论应该792，这种程度的偏差合理吗？如果基 ...

GAFF是形式简单的通用力场，能达到这样的准确度就可以了。如果你想要更高精度，那需要你去优化分子专一性的LJ参数/原子电荷，类似于优化水模型一样。届时可以以当前的参数作为初猜。

作者
Author: zdworld 时间: 2023-10-3 16:51

喵星大佬发表于 2023-10-3 12:13
合理，这误差已经不大了，就2.2%你还想咋样
gaff的原子类型考虑的化合物类型很多，vdW参数这点误差已经 ...

好的，非常感谢

作者
Author: zdworld 时间: 2023-10-3 16:52

sobereva 发表于 2023-10-3 14:17
GAFF是形式简单的通用力场，能达到这样的准确度就可以了。如果你想要更高精度，那需要你去优化分子专一性 ...

谢谢，我去了解一下

作者
Author: nku_linz 时间: 2023-10-14 12:59

sobereva 发表于 2023-9-29 05:52
CS直接用sobtop创建基于GAFF力场的拓扑文件就完了，根本不需要也不应该基于Hessian产生谐振势形式的力常数
...

谢谢老师指出错误，我没有仔细看FAQ10和FAQ11, 折腾不讨好！已经改用例2的方法产生.itp文件，MD已可以正常跑起来。

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