模拟氮、氧掺杂石墨烯时，能量极小化报错

a9471163 · 发表于 Post on 2021-12-15 13:58:06

本帖最后由 a9471163 于 2021-12-15 14:11 编辑

各位老师，最近有个氮、氧掺杂石墨烯（氧掺杂比例高，达到30%，氮极少，0.6%左右的掺杂比）在gaff力场下用gromacs模拟，拓扑文件是acpype产生的，原子电荷是universal力场优化石墨烯结构以后（因为原子数有500多个，没有用DFT做结构优化），用orca5.0的b97-3c在隐式SMD水模型下计算了单点能，然后用multiwfn最新版计算resp电荷（即：载入文件，选择选项7-18-1，使用multiwfn自带算法进行，计算过程未调用高斯），但所得的RESP电荷个人推测是有点问题。
因为测试发现一旦使用这套RESP电荷值，能量极小化就提示：未达到收敛限，但能量极小化无法继续，自动停止（如图），此时原子上的最大受力仍然有10的5次方，强行跑后续模拟则会提示很多lincs warning最后崩溃，如果所有原子的原子电荷改成0，则可以正常模拟（但显然这么模拟不准确），所以求助各位老师，此时如何计算得到可用的电荷值，以便模拟可以正常跑呢？

继续测试发现，若氧化石墨烯中氧掺杂比例从30%降低到10%以下，则使用同样一套模拟流程是可以正常跑模拟的，正常的氧化石墨烯（简称GO）和不正常的GO我都打包上传了，里面有计算过程用到的全部文件和参数，可以清楚地看到所用GO的构型（请看压缩包的mol2文件）。有没有大佬可以指导一下我那个不正常的GO是不是建模的结构有问题呢，还是其他原因造成的？

个人已经尽力检查建模结构是否合理，是通过这几种方式判断的：
1、根据化学常识，观察并判断原子的排列和连接方式，没发现异常
2、根据文献，羧基布置在GO边缘并且在该pH下全部质子化，羟基、环氧基、N原子布置在GO的中间，这也与文献一致
3、根据acpype的Check Weird Bonds检查键连关系，没有发现有不正确的成键方式
4、计算了分子的净电荷，由于羧基去质子化带来15个负电荷
5、我在模拟前一般会通过xtb或universal力场优化分子结构，若优化后分子明显有不合理的构型，说明初始结构不合理，我会积极调整，但本次用universal优化结构，没有发现不正常现象
6、根据实验表征结果，实验证明氧掺杂比例达到30%是可能发生的并且此时氮含量较低，说明所用模型中的掺杂比例没有错

abnormal-and-normal-GO.zip (4.84 MB, 下载次数 Times of downloads: 35)

sobereva · 发表于 Post on 2021-12-16 14:57:21

我没时间具体看你的文件，但光是这么一条

计算了分子的净电荷，由于羧基去质子化带来15个负电荷

就意味着结果没意义。那么多负电荷之间库仑互斥极强，体系自然不可能稳定。如果是在水中模拟因而羧基自发质子解离，也应当考虑抗衡离子或者引入H3O+
建议你先看看已发表的文章里算氧化石墨烯都是用的什么模型。对于常见体系，自己做之前一定先多看看前人的文章免得瞎算、踩坑

a9471163 · 发表于 Post on 2021-12-16 16:45:52

本帖最后由 a9471163 于 2021-12-16 16:54 编辑

sobereva 发表于 2021-12-16 14:57
我没时间具体看你的文件，但光是这么一条
就意味着结果没意义。那么多负电荷之间库仑互斥极强，体系自然不 ...

谢谢卢老师提醒，我的模型中羧基有15个，如果分子很小的话，这么多羧基可能还是有问题的；如果分子足够大了，是可以容纳这么多羧基的，比如这篇文献：New Insight into the Aggregation of Graphene Oxide Using Molecular Dynamics Simulations and Extended Derjaguin−Landau−Verwey−Overbeek Theory，他布置的羧基比我的还多（如图）

您的提醒让我想起了最关键的一点，可能不是边缘羧基的问题（我跑能量极小化之前也引入抗衡离子Na+了），而是我的模型中，碳氧个数比不合理，我看到文献计算氧化石墨烯，用的比例是C20 O2 (OH)2 (COOH)1，即文献C/O=21:6=3.5，而我的模型C/O=1.69。实验XPS表征证明，所合成的材料的C/O=2.34，我的C/O偏离了实验值和文献模拟氧化石墨烯的值，导致过多的O带来很大的静电斥力
但遗憾的是，为了验证我的想法即C/O不合理导致模拟崩溃，我又建立了一个C239 H82 N3 O77的模型，此时C/O降低到了3.10，如果按实验表征来说的话，这个比例是可能存在的，但实际模拟就遇到能量极小化可以顺利完成，正式NPT模拟后提示很多lincs warning最后崩溃的问题，我目前打算换个gromos力场试试，或许是gaff不适合模拟这种分子内静电斥力大一些的体系

sobereva · 发表于 Post on 2021-12-16 17:00:22

原理上GAFF描述这个体系没问题

最好把问题简化，是否考虑边缘的羧基、是否考虑环氧基、是否考虑羟基，全都进行测试，搞清楚问题在哪里，然后细致分析

a9471163 · 发表于 Post on 2021-12-16 17:07:18

sobereva 发表于 2021-12-16 17:00
原理上GAFF描述这个体系没问题

好的，谢谢您，有人也建议我把SMD水溶剂模型改成CPCM（water）计算单点能，并且计算RESP2代替RESP，我也正在计算，分子大了一些，在等结果，不知道能不能对模拟有改善，听说是orca中SMD模型会严重影响模型的收敛情况，换成CPCM就会好很多

a9471163 · 发表于 Post on 2021-12-17 12:34:57

sobereva 发表于 2021-12-16 14:57
我没时间具体看你的文件，但光是这么一条
就意味着结果没意义。那么多负电荷之间库仑互斥极强，体系自然不 ...

sob老师，想请教一下gaff力场下模拟使用Qeq电荷准确度能接受吗，我只知道Qeq适合跟universal力场搭配，现在我这个体系使用Qeq电荷是可以跑模拟的

另外，您说的很有道理，我领会了您的意思，后面会试一下先引入抗衡离子再去计算单点能和RESP，而不是设置成-15的负电去计算单点能，因为化学价是-15，但实际电荷值不会那么负的，我之前计算完RESP观察钠离子也就是0.8左右的电荷值，而不会到1

sobereva · 发表于 Post on 2021-12-17 12:44:53

a9471163 发表于 2021-12-17 12:34
sob老师，想请教一下gaff力场下模拟使用Qeq电荷准确度能接受吗，我只知道Qeq适合跟universal力场搭配，现 ...

QEq太烂，对静电势描述太差

显然不能带着抗衡离子算RESP电荷，RESP电荷是对氧化石墨烯自己而言计算的

a9471163 · 发表于 Post on 2021-12-17 12:54:59

sobereva 发表于 2021-12-17 12:44
QEq太烂，对静电势描述太差

显然不能带着抗衡离子算RESP电荷，RESP电荷是对氧化石墨烯自己而言计算的

好的老师，测试发现这个崩溃，确实很有可能是一个不大的分子基面去承受-15个负电，负电荷太多了，导致排斥大所以崩溃的，但根据文献，又觉得这个模型上长15个羧基也是合理的
目前有两条路：
1、有文献是让羧基全部去质子化的，这种做法适合羧基少的模型，比如3-5个羧基，否则基面上负电太多，容易崩溃；也有文献是让一部分羧基去质子化去模拟的，我可以试试后面那条路，也许适合我这个模型

2、能否把氧化石墨烯没有去质子化的状态拿去计算RESP，算完以后把所有的羧基H删了，那个H也就带0.X的正电，所以剩下的基面上，负电也就明显弱了一些，估计模拟不会崩溃，您看是否值得试一下这种方法

sobereva · 发表于 Post on 2021-12-17 21:14:49

a9471163 发表于 2021-12-17 12:54
好的老师，测试发现这个崩溃，确实很有可能是一个不大的分子基面去承受-15个负电，负电荷太多了，导致排 ...

1 原理上，genion会把抗衡阳离子加到静电势最负的位置，就算羧基多，Na离子如果在羧基附近充分，也能很大程度抵消静电互斥。一部分去质子化也不是不可以，去质子化的羧基之间的静电相互作用本身也会影响羧基的pKa使得去质子化状态稳定化程度降低

2 原理上不能。解离下来的质子带一个单位正电荷，相应地在羧基上留下一个单位负电荷。要和实际相符

a9471163 · 发表于 Post on 2021-12-18 10:51:57

sobereva 发表于 2021-12-17 21:14
1 原理上，genion会把抗衡阳离子加到静电势最负的位置，就算羧基多，Na离子如果在羧基附近充分，也能很大 ...

老师您好，我那个氧化石墨烯体系，想请教一下有没有办法通过调整参数，继续让模拟更稳定，昨天看您的PPT把步长从2飞秒改到1飞秒，效果明显，可以延长模拟时间，最后模拟5纳秒，最后仍然崩溃，每次崩溃前都有共性报错是：One or more water molecules can not be settled（用或不用lincs基本都是这个时间段模拟崩溃，都是这个报错）。如果能模拟10纳秒以上我也就满足了，下面这是我的mdp文件（请看附件）

我那个GO又对比了几篇文献，觉得自己的模型合适，文献如果是把一部分羧基去质子化，那也是只有10%的羧基不做去质子化（如下图），我15个羧基，就相当于只有1个羧基不做去质子化，那基面就带-14的负电，这跟我之前计算用的-15相差无几，而且我RESP电荷计算方法也无误，模拟时抗衡离子也加的非常规范。所以就想看看能不能继续调整NPT模拟的参数，让这个模拟能更稳定，非常感谢老师

mdout.mdp (11.53 KB, 下载次数 Times of downloads: 4)

sobereva · 发表于 Post on 2021-12-18 22:20:38

a9471163 发表于 2021-12-18 10:51
老师您好，我那个氧化石墨烯体系，想请教一下有没有办法通过调整参数，继续让模拟更稳定，昨天看您的PPT ...

用柔性水。此时不牵扯到用settle约束

a9471163 · 发表于 Post on 2021-12-19 11:56:06

sobereva 发表于 2021-12-18 22:20
用柔性水。此时不牵扯到用settle约束

非常感谢sob老师，我今天试试

a9471163 · 发表于 Post on 2021-12-19 23:47:34

sobereva 发表于 2021-12-18 22:20
用柔性水。此时不牵扯到用settle约束

老师您好，已经尝试使用柔性水进行正式的NPT模拟，结果如下：1、模拟前做能量极小化的方法是：使用cg算法，手动修改能量极小化的收敛限，逐渐降低收敛限，然后用上一次优化完的结构继续、反复做能量极小化，最后最大力能优化到小于700
2、进行限制性分子动力学模拟，模拟步长1飞秒报错，使用0.1飞秒跑完了该模拟，但苯环已经扭曲变形（如图，自知继续NPT正式模拟已无意义，但还是稍微尝试了一下，果不其然NPT模拟报错，报错内容如下）
3、开始进行正式NPT模拟，使用柔性水，报错内容如下（模拟步长0.5飞秒也不行，中途崩溃，想请问一下没必要尝试强行模拟了吧，测试到这一步，是不是可以理解成就是模型不合理，导致RESP电荷计算的不合理）

喵星大佬 · 发表于 Post on 2021-12-20 07:38:37

本帖最后由喵星大佬于 2021-12-20 07:59 编辑

仔细检查你的成键项参数，尤其是涉及环氧的

扭曲到这个程度不可能是静电作用计算的问题，一定是成键参数有问题，盲猜是环氧的键角和二面角项

另外，一般的分力场(尤其是键角仅用谐振势的)对于这种高张力小环是描述很差的，结果普遍不可信。

你仔细看你的图，有问题的地方基本都是因为环氧的C试图变成正四面体构型而导致的扭曲，或者是由于环氧本身键角太小而模拟时另外三个C与O都在同一个半球(没有imprepor项限制)导致的。描述环氧这种需要专门的原子类型和成键项参数，包括限制另外三个原子进入同一半球的improper二面角项，不记得gaff有没有这个了，印象中是没有，基本可以确定这块一定是有问题的。最后就是看你要模拟的是什么性质了，如果实在搞不定参数可以用gfn-ff搞搞看

a9471163 · 发表于 Post on 2021-12-20 11:39:22

喵星大佬发表于 2021-12-20 07:38
仔细检查你的成键项参数，尤其是涉及环氧的

扭曲到这个程度不可能是静电作用计算的问题，一定是成键参数 ...

非常感谢您的提醒，gfn-ff计算速度很慢，想请教一下OPLS-AA能不能做这种体系，有文献是OPLS-AA做的氧化石墨烯，但环氧的比例比我的模型少，我这个是生物炭，按实验表征的30%的氧比例建模的，这时候O掺杂是到了极限的，几乎没法再加上去更多的氧

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

[GROMACS] 模拟氮、氧掺杂石墨烯时，能量极小化报错

评分 Rate

浏览过的版块