sobtop+RESP.sh生成碳酸氢根gaff力场，模拟崩溃

Changtx

各位老师同学好！我最近在制作一个批量生成gaff力场的脚本，使用openbabel将SMILES转换成pdb文件，通过改写RESP.sh脚本，调用Gaussian进行结构优化以及单点能的优化，将优化后的结果导出成pdb文件，使用openbabel转化成mol2文件，调用Multiwfn软件计算RESP电荷，将mol2文件导入到sobtop软件中，按照例2进行操作。上述操作均正常结束。进行了一批测试操作，有MDEA,PZ,DMF,HCO3等，在进行后续gromacs模拟过程中发现模拟崩溃，因为一直做这个体系没出现崩溃，只不过更改了力场，之前使用OPLS力场，觉得可能是力场参数有问题，对力场参数进行了检查并且一步步进行简化体系。最后发现HCO3的力场参数有问题，体系内含有HCO3就会崩溃，崩溃模式为能量最小化90步就停止，后续进行nvt模拟出现Too many LINCS warnings.HCO3的itp文件附件 体系为2000个HCO3，2000个H2O，2000个MDEAI+，em，nvt文件附件所示。体系除去HCO3后体系正常模拟。各位老师请指导下哪里出了问题。[url=]附件 Attachment[/url]

牧生

你操作应该没问题，我遇到过几乎相同的问题。

这里的原因是电荷问题，但我也没法解决。我试过很大盒子里面就只放一两个碳酸根，也照样不行。

你可以在我以前的帖子的讨论里面找找灵感，看能不能解决。

http://bbs.keinsci.com/thread-46947-1-1.html

Huschein

牧生 发表于 2025-3-3 21:40
你操作应该没问题，我遇到过几乎相同的问题。

这里的原因是电荷问题，但我也没法解决。我试过很大盒子里 ...

我看了你的帖子，是不是电荷太大容易造成体系震荡？-2的电荷分给每个原子就会很多，我看你那里三个-1一个+1的电荷，模拟过程中可能瞬时静电力变大 吸引过来然后vdw又瞬间排斥，这样的话体系就会不稳定。包括你那个steep不行但是cg可以，这很可能就是steep会一下子跳出极小值。你这个没有太好的解决方法，最好就是用专门优化的碳酸根-水 力场，他们可能对vdw和电荷专门参数化避免了这些问题，你这个自己的电荷加上GAFF的vdw效果未必好

18217265596

amber我模拟过碳酸氢根 没啥问题 应该是你参数的问题

Changtx

18217265596 发表于 2025-3-4 09:26
amber我模拟过碳酸氢根 没啥问题 应该是你参数的问题

老师能分享下你的HCO3的itp吗，我尝试一下。

18217265596

Changtx 发表于 2025-3-4 11:57
老师能分享下你的HCO3的itp吗，我尝试一下。

我是amber哦 你直接拿gaff2用amber的脚本作一套 是一样的

chuanghao

根据sob老师sobtop网站上的描述“比如曾经我用某个程序产生带磺酸基的分子的拓扑文件，看似拓扑文件没任何毛病，但跑起来就是会崩。后来一帧帧地观看模拟开始后的结构变化，发现是磺酸基上的氢在一瞬间跑到距离氧很近的地方去了。分析发现是由于这个基团上的氢带的电荷很正，氧的电荷很负，它们彼此间有1-4静电作用项因此有很强的静电吸引，再加上氢的质量很小、很容易运动，导致了上述看到的问题。像这种情况，可以把1-4作用项去了（从[ pairs ]中删除对应的），也可以把氢改成质量更大的重氢（在[ atoms ]里改原子质量）减缓其运动，或者用很小的模拟步长。显然，从拓扑文件和参数中找出原因必须对拓扑文件里的各种项、力场的函数形式都理解得很清楚才行。”

这里你可能也是氢原子的电荷和连接它的氧原子电荷相差很大，去掉pair项里和这个H相关的项可能就行了

Changtx

chuanghao 发表于 2025-3-4 12:58
根据sob老师sobtop网站上的描述“比如曾经我用某个程序产生带磺酸基的分子的拓扑文件，看似拓扑文件没任何 ...

好的我去尝试一下，感谢老师了。

牧生

HCO3.itp (2.73 KB, 下载次数 Times of downloads: 15)

2025-3-5 08:57 上传 Uploaded

点击下载
Click to download

试了下，碳酸氢根没问题。

Changtx

牧生 发表于 2025-3-5 08:57
试了下，碳酸氢根没问题。

我也试了下，您的力场文件是可以的。对比了下我俩的itp文件，区别在于电荷有点区别，可能是我生成的电荷刚好不合适。