计算化学公社

标题: oplsaa力场搭配RESP2(0.5)电荷模拟多孔生物炭报错 [打印本页]

作者
Author: a9471163 时间: 2022-2-1 21:59
标题: oplsaa力场搭配RESP2(0.5)电荷模拟多孔生物炭报错
本帖最后由 a9471163 于 2022-2-1 22:18 编辑

各位老师春节快乐、虎年大吉，我使用gromacs2018.4和2019.6版本，使用oplsaa力场搭配RESP2(0.5)电荷模拟多孔生物炭时都报错误（报错内容相同，总结在下面，模拟参数如图，已上传图片。所以我觉得不是gromacs版本的问题。多孔生物炭类似于氧化石墨烯，但基面上的缺陷比较多）
随后我使用materials studio，利用forcite模块中的compass力场分配的电荷，替换了RESP2(0.5)电荷，这套电荷搭配oplsaa力场，在gromacs软件中运行动力学模拟是可以正常进行的，因此报错的原因被锁定在RESP2（0.5）电荷的问题（RESP电荷也是报错的）

oplsaa力场模拟这类碳材料的稳健性已有文献证明，如New Insight into the Aggregation of Graphene Oxide Using Molecular Dynamics Simulations and Extended Derjaguin-Landau-V erwey-Overbeek Theory，因此我个人认为力场参数的问题不大（生物炭拓扑文件17.itp上传了），从能量极小化后的输出结构来看，各种C、H、O、N原子的描述较好，其最终构型与UFF力场或xtb优化得到的构型类似，没有观察到明显的不合理变形现象发生

分子模型文件是final6.mol2，我已上传。还请大佬把把关看下模型有无问题，我个人认为分子模型合理，设置模型时的元素比例参考了实验数据，官能团位置和种类参考了文献，

观察compass力场的电荷和RESP2发现有一定的区别，RESP2那套电荷的绝对值明显更大，很多原子的电荷值大于1，而compass力场的电荷的绝对值都比较小，两套电荷值我都上传了，文件是charge.txt。生物炭拓扑文件是17.itp

gromacs中模拟的报错内容总结如下：
正式模拟NPT时，出现很多lincs warning，警告涉及的原子也不固定，相当于是遍地开花，碳基面上到处都有警告，最终因为too many lincs warning而模拟失败。为了解决报错，我尝试了这些办法

1、不再使用lincs ，即设置constraints=none 并且模拟步长从2飞秒缩短到1飞秒，此时出现新的报错，内容如图，意为water cannot be settled
2、在上面一步的基础上使用柔性水试图解决报错，此时water cannot be settled就没有了，但有新的报错，如图，意为nonbonded interaction larger than table limit，对此我没有什么办法
3、尝试仍然使用lincs，但lincs-warnangle从默认的30度增大到90度，这也是软件允许设置的最大值，此时就没有lincs warning，但不一会儿模拟就仍然崩溃，主要报错是cannot fix pbc，如图

上述报错都是在模拟很短时间时就出现的

值得一提的是，能量极小化和限制性分子动力学（NPT）时，也有错误（总结如下），我是自己勉强解决了报错，才走到了正式模拟这步，但正式模拟目前没什么办法能进行下去，如果大佬有办法还请支支招吧，感谢~

能量极小化是遇到力越优化越大，根本不能收敛的问题，碳基面上的羟基H，由于静电吸引作用，与碳基面上碳原子可能出现吸引，最后羟基H跟C原子位置重叠，受力非常大，能量极小化失败（对应的是黑色背景图，图上标示了原子类型和RESP2电荷值）

我用XTB优化调整了生物炭的初始结构，拉远了这个羟基和碳原子之间的距离，避免它们吸引在一起，最后能量极小化能够在Fmax小于100的精度下收敛了

随后进行限制性分子动力学把水分子弛豫开，此时position_restraints中fcx=fcy=fcz=1000，此时跟正式模拟报错一样，多个原子上都会出现lincs warning，解决办法是lincs-warnangle设置为90，此时只有一对原子上会一直警告说lincs会超过90度，无奈之下把这对原子的fcx=fcy=fcz=1000调整为10000，此后限制性分子动力学可以完成了

计算RESP2（0.5）时我觉得操作无误，简述操作步骤，大佬们把把关吧，感谢~：

使用的分子模型就是同一个final6.mol2，羧基去质子化，电荷为-3
考虑到原子数较多，没有用DFT优化几何结构，而是在XTB中的extreme精度下添加H2O溶剂模型进行了几何优化，
随后用优化了的几何结构在ORCA中使用 B3LYP D3 def2-SVP def2/J RIJCOSX noautostart miniprint nopop计算真空和水里的两个单点能，溶剂模型是默认的SMD
产生的输出文件提交给Multiwfn计算真空、水里两个RESP电荷并取平均值得到RESP2，计算RESP的方法是参考sob老师《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》一文中天冬氨酸的示例进行，除了读取构象列表没有进行外其他步骤相同，电荷约束施加在羧基的COO三个原子上，约束为-1，等价性约束施加在羧基的两个O上

我目前最关心的是，RESP2（0.5）的稳健性大部分情况下都是没问题的，我这里模拟不下去，是继续使用RESP2但进行一些必要的调整，还是更换一套电荷用于gromacs模拟呢？谢谢各位老师

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

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

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

(, 下载次数 Times of downloads: 26)
(, 下载次数 Times of downloads: 34)
(, 下载次数 Times of downloads: 20)
(, 下载次数 Times of downloads: 26) (, 下载次数 Times of downloads: 19)

作者
Author: sobereva 时间: 2022-2-2 01:18
如我群里说的，这俩原子电荷一个很正一个很负，再加上氢的质量小因此移动快，如果考虑它俩的静电相互作用，可能导致有严重吸到一起的倾向导致出现问题。去掉[pairs]里相应项不计算1-4作用就可以避免。另外，羧基部分没有绝对必要非得总电荷约束为-1。去掉这个约束的话说不定氧上的电荷会更负一些，而那个悬着的碳的电荷没那么负了。
另外，GFN-xTB优化这样的体系所得结构精度还是比较牵强。用B97-3c在CPCM下做个优化也不至于花太多时间（几何优化收敛限可以设得比较宽松）。

作者
Author: a9471163 时间: 2022-2-2 21:28

sobereva 发表于 2022-2-2 01:18
如我群里说的，这俩原子电荷一个很正一个很负，再加上氢的质量小因此移动快，如果考虑它俩的静电相互作用， ...

谢谢sob老师，您说的方法在尝试了，后面看一下效果

作者
Author: 喵星大佬 时间: 2022-2-3 06:15
我觉得一般来说不会出现那么大的原子电荷

作者
Author: a9471163 时间: 2022-2-5 15:01

sobereva 发表于 2022-2-2 01:18
如我群里说的，这俩原子电荷一个很正一个很负，再加上氢的质量小因此移动快，如果考虑它俩的静电相互作用， ...

非常感谢sob老师，按您说的做模拟可以正常完成

作者
Author: a9471163 时间: 2022-2-5 15:01

喵星大佬发表于 2022-2-3 06:15
我觉得一般来说不会出现那么大的原子电荷

您说的对，按sob老师的办法重算RESP2之后，原子电荷的绝对值变小了

作者
Author: bxc 时间: 2022-6-28 22:47
各位老师好，请教一下，您这个石墨烯就总共就316原子吗，如果我是拿氧化石墨作为基板，需要20*3nm，2000多原子，请问我该如何计算resp2（0.5）电荷，orca算不动呀@sobereva

作者
Author: sobereva 时间: 2022-6-29 01:23

bxc 发表于 2022-6-28 22:47
各位老师好，请教一下，您这个石墨烯就总共就316原子吗，如果我是拿氧化石墨作为基板，需要20*3nm，2000多 ...

整个体系本来就算不动
氧化石墨烯上不同位点适合的原子电荷在文献里绝对能搜到

作者
Author: bxc 时间: 2022-6-29 10:40

sobereva 发表于 2022-6-29 01:23
整个体系本来就算不动
氧化石墨烯上不同位点适合的原子电荷在文献里绝对能搜到

那是不是用sobtop定义好每种原子类型，同时定义好原子电荷，然后直接生成itp就可以啦，想问下老师这是最快捷的方法了吗@sobereva

作者
Author: sobereva 时间: 2022-6-30 00:50

bxc 发表于 2022-6-29 10:40
那是不是用sobtop定义好每种原子类型，同时定义好原子电荷，然后直接生成itp就可以啦，想问下老师这是最 ...

这种体系应当在assign_AT.dat里定义原子类型判断规则，并且各种原子类型用的原子电荷直接写文献里的值，到时候sobtop就能根据你的体系结构通过assign_AT.dat里的规则自动指认原子类型和赋予原子电荷了。

作者
Author: bxc 时间: 2022-6-30 21:25

sobereva 发表于 2022-6-30 00:50
这种体系应当在assign_AT.dat里定义原子类型判断规则，并且各种原子类型用的原子电荷直接写文献里的值， ...

谢谢老师

作者
Author: jay6013 时间: 2023-6-9 17:08

a9471163 发表于 2022-2-5 15:01
非常感谢sob老师，按您说的做模拟可以正常完成

LZ可以私聊我一下吗？因为新号没有权限发私信，希望可以向你学习一下

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