计算化学公社

标题: 求助：Sobtop生成GAFF+UFF力场参数模拟镍基MOF，能量最小化报错 [打印本页]

作者
Author: a9471163 时间: 2023-3-30 18:57
标题: 求助：Sobtop生成GAFF+UFF力场参数模拟镍基MOF，能量最小化报错
本帖最后由 a9471163 于 2023-3-30 20:53 编辑

老师，您好，Sobtop生成GAFF+UFF力场参数模拟镍基MOF（CCDC号：924088），能量最小化报错（镍用UFF参数，其他原子是GAFF）。此外多次尝试表明，一、如用Sobtop生成UFF力场参数，也是同样的报错，
二、原子电荷如果不使用，或者用QEq电荷，都会是同样的报错，（由于计算较慢目前尚未使用REPEAT电荷），个人认为不是电荷的问题
三、能量最小化时，integrator = cg或者integrator = steep都试过了，都是同样的报错
四、给gromacs做能量最小化的初始结构，如果直接用CCDC数据库里的CIF结构，或者用UFF力场优化后的CIF结构做能量最小化，都是同样的报错（Materials Studio（MS），forcite模块选UFF力场，精度是ultra-fine做优化，MS软件输出的提示是优化成功；根据MS的计算结果，作者CCDC号：924088提供的CIF如果不做结构优化，也是一个稳定结构，个人认为该结构可以直接用于能量最小化）
我的具体操作过程如下：

1、下载CIF文件，CCDC号：924088，MOF代号TKL-107
2、GaussView载入该CIF，这是一个MOF单胞，检查成键，发现成键全部正确，如图1

(, 下载次数 Times of downloads: 5)
图1：GaussView显示的TKL-107的成键

3、GaussView将CIF扩胞为2*2*2大小，这是因为单胞太小，扩胞后晶胞参数为，a=b=30.3476, c=37.7180（单位：埃），α=β=90°，γ=120°，原子数：1296，原子种类：氟、碳、氮、氧、镍
4、GaussView将2*2*2超胞导出为mol2文件，文件已上传，文件名为924088.mol2

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

5、将该mol2文件载入Sobtop，根据Sobtop官网http://sobereva.com/soft/Sobtop/#ex10的例10，分别按键盘的1,2,4，回车，回车，生成基于GAFF+UFF力场的ITP和TOP文件。
6、把924088.mol2提交到GaussView，另存为一份924088.pdb。
7、用gmx editconf -f 924088.pdb -o system-initial.pdb -box 3.03476 3.03476 3.7718 -angles 90 90 120 生成能量最小化的初始结构文件，system-initial.pdb。该文件已经上传。

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

8、将system-initial.pdb载入VMD观察结构，确定结构无误。如图2所示

(, 下载次数 Times of downloads: 8)
图2：VMD显示的能量最小化前的初始结构，图中显示了+x和+y方向的镜像

9、基于上述生成的ITP和TOP文件，以及system-initial.pdb，做能量最小化。首先使用gmx solvate -cp system-initial.pdb -o system-water.gro -p topol.top 命令给体系加水，其中水分子用的是amber14sb中的SPCE模型。

10、这是电中性MOF，无论是使用QEq原子电荷，还是不使用原子电荷，体系都无需添加用于平衡电荷的离子。

11、做能量最小化，命令为gmx grompp -f em.mdp -c system-water.gro -p topol.top -o em.tpr -maxwarn 1 加上 gmx mdrun -v -deffnm em -ntmpi 1，其中忽略的警告为：
WARNING 1 [file topol.top, line 38]:
  1296 non-matching atom names
  atom names from topol.top will be used
  atom names from system-water.gro will be ignored，意思是，1296个MOF原子都不对应。
根据图3，可知，原子是对应的，只是后缀序号不同导致警告（Ni1 in itp文件-Ni in gro文件, C2 in itp文件-C in gro文件, C3 in itp文件-C in gro文件...）所以该警告我忽略了。

(, 下载次数 Times of downloads: 4)
图3：gro和itp文件中原子的对应关系展示

在能量最小化时有多个警告，如图4所示，因此能量最小化无法进行。我在文章开头提到的多次尝试，每次也都是报这个错误。

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

图4：能量最小化先出现警告，
(, 下载次数 Times of downloads: 4)

图4：能量最小化随后失败，自动停止

我的GROMACS 版本: 2020.3-MODIFIED，即2020.3 CUDA GPU加速版，AVX指令集，该版本下载于《GROMACS的原生Windows版的编译和安装方法》（http://sobereva.com/458）。此外，2019.6版也试了，有相同的报错
最后，附上能量最小化用的ITP、TOP、mdp文件，谢谢老师指教！ (, 下载次数 Times of downloads: 19)

作者
Author: ggdh 时间: 2023-3-31 13:00
优化一下平衡键长键角

作者
Author: a9471163 时间: 2023-3-31 16:01

ggdh 发表于 2023-3-31 13:00
优化一下平衡键长键角

好的，谢谢老师，是必须用量化计算键长键角么，我模拟的两个MOF都是一样的问题，其中有个MOF晶胞不小，怕是算不动

作者
Author: ggdh 时间: 2023-3-31 19:55

a9471163 发表于 2023-3-31 16:01
好的，谢谢老师，是必须用量化计算键长键角么，我模拟的两个MOF都是一样的问题，其中有个MOF晶胞不小，怕 ...

可以不算
你直接把晶体的键长键角当作平衡键长键角

作者
Author: a9471163 时间: 2023-3-31 20:28

ggdh 发表于 2023-3-31 19:55
可以不算
你直接把晶体的键长键角当作平衡键长键角

非常感谢，有思路计算了，我查到ztop非常适合这方面计算，看了说明后发现，“程序会在当前目录找一个名为DON.log的文件”，我没有用高斯计算，就没这个log文件，请问是否可以用mol2等其他文件代替log文件呢，如果需要手写log文件，要写成什么格式呢

作者
Author: ggdh 时间: 2023-3-31 20:47

a9471163 发表于 2023-3-31 20:28
非常感谢，有思路计算了，我查到ztop非常适合这方面计算，看了说明后发现，“程序会在当前目录找一个名为 ...

直接cp2k算，用sobtop把cp2k算出来的平衡键长键角直接写入应该也可以把？

作者
Author: a9471163 时间: 2023-3-31 21:06

ggdh 发表于 2023-3-31 20:47
直接cp2k算，用sobtop把cp2k算出来的平衡键长键角直接写入应该也可以把？

好的，谢谢老师，我试试

作者
Author: GoldenBaby 时间: 2023-3-31 23:12
我觉得你检查一下你生成的mol2文件吧，我这打开你上传mol2文件是这个样子
(, 下载次数 Times of downloads: 4)

作者
Author: a9471163 时间: 2023-3-31 23:22
本帖最后由 a9471163 于 2023-3-31 23:23 编辑

GoldenBaby 发表于 2023-3-31 23:12
我觉得你检查一下你生成的mol2文件吧，我这打开你上传mol2文件是这个样子

谢谢，这个是3D MOF。不是层状MOF，您看到成键是斜着连，应该也是对的，因为那个CIF在GaussView里成键显示是完全对的，mol2通过CIF一步转换生成，GaussView软件很智能，我还从来没见过在GaussView显示CIF成键完全正确的情况下，会生成错误的mol2，我相信GaussView不会这么令人失望的

作者
Author: GoldenBaby 时间: 2023-3-31 23:44
本帖最后由 GoldenBaby 于 2023-4-1 00:18 编辑

a9471163 发表于 2023-3-31 23:22
谢谢，这个是3D MOF。不是层状MOF，您看到成键是斜着连，应该也是对的，因为那个CIF在GaussView里成键显 ...

我理解怎么回事了。抱歉。不过我刚刚又看了一眼你这个结构，你那个镍上是不是多了个氧，那个应该是个配位水分子吧，要不你看看是不是这个地方有问题？

作者
Author: a9471163 时间: 2023-4-1 12:56
本帖最后由 a9471163 于 2023-4-1 12:58 编辑

GoldenBaby 发表于 2023-3-31 23:44
我理解怎么回事了。抱歉。不过我刚刚又看了一眼你这个结构，你那个镍上是不是多了个氧，那个应该是个配位水 ...

非常感谢大佬，您这里说的对，我按您说的给氧上补了2个氢，那里应该是配位水。多亏了您的提醒，我检查了之前自己做的几个MOF，那几个MOF都正常跑MD模拟了，但是有的氧上也是没补氢，本该补氢形成配位水的，可见如果MOF能跑动力学，这里不补氢的话，也能跑，如果不能跑模拟的话，这里就算补了氢，也不能跑，因此，结论是这里补氢对解决能量最小化无法进行的问题没有帮助。我还试了不在水盒子里做能量最小化，而是只有一个MOF在真空里能量最小化，也是一样的情况，最大力只能优化到10的5次方，远远达不到收敛标准

作者
Author: a9471163 时间: 2023-4-3 19:06
本帖最后由 a9471163 于 2023-4-3 19:07 编辑

ggdh 发表于 2023-3-31 19:55
可以不算
你直接把晶体的键长键角当作平衡键长键角

老师，我有了最新进展，总结在了http://bbs.keinsci.com/thread-36296-1-1.html，您有空的话可以看看，目前我通过控制变量法，逐个排查可能的问题，觉得您说的平衡键长键角不准，是很有可能的，我cp2k算了一次结构优化，第一次跑了22小时也没跑完，收敛趋势有振荡，我看可能收敛无望，超算里充的钱也不多了就先终止计算了，我今天晚上充值后继续算第二次，微调下计算参数

我还想请教一下，您说不计算，把晶体的键长键角当作平衡键长键角，我试了，我如果拿PDB文件载入sobtop，生成的拓扑能做能量最小化，如果是mol2载入sobtop生成的拓扑没法能量最小化，这俩都是把晶体的键长键角直接用的，作者报道的值。所以说就算把实验值拿来用了，可能还是有点问题，不知道为什么。所以我现在就只能是指望cp2k的计算值能不能解决问题了，但在我心目中，计算值哪怕是量化计算值，应该也不如单晶XRD实验值吧，不知您怎么看这个问题。详情就在上面那个链接里总结了

作者
Author: 喵星大佬 时间: 2023-4-4 12:24

a9471163 发表于 2023-4-3 19:06
老师，我有了最新进展，总结在了http://bbs.keinsci.com/thread-36296-1-1.html，您有空的话可以看看，目 ...

你要搞清楚一个问题，在能量极小化的结构中，并不一定要求对于能量的每一项也要极小，换句话说，如果一种化学键平衡键长是a，一个结构里面同类化学键有b>a和c<a两种键长，这和整体能量最低并不冲突
换句话说，事实上在力场中直接指认晶体的键长/键角为该变量的极小点本身就是错误的
力场参化过程中都是通过构造一系列模型化合物来拟合键长/键角/二面角等成键参数的

作者
Author: a9471163 时间: 2023-4-4 13:25

喵星大佬发表于 2023-4-4 12:24
你要搞清楚一个问题，在能量极小化的结构中，并不一定要求对于能量的每一项也要极小，换句话说，如果一种 ...

好的，谢谢大佬，“力场参化过程中都是通过构造一系列模型化合物来拟合键长/键角/二面角等成键参数的”，我如果用cp2k做结构优化，得到平衡键长键角，是不是就是在完成您说的这步呢。我目前正在计算这个结构优化

作者
Author: WuQs 时间: 2023-4-4 16:06
你可以检查下你失败的时候产生的log文件里的具体信息提示，我之前遇到过的跟你一样，但是我在log文件里找到具体出错的位置了；提示在某个原子附近异常，然后就找到具体问题了。我浏览了一下你似乎还没排出是什么问题，可以先去log里看一下。

作者
Author: a9471163 时间: 2023-4-4 19:50

WuQs 发表于 2023-4-4 16:06
你可以检查下你失败的时候产生的log文件里的具体信息提示，我之前遇到过的跟你一样，但是我在log文件里找到 ...

谢谢老师，我后续检查进一步缩小了范围，目前最有可能的是GaussView转的mol2有问题，我把测试结果汇总在了http://bbs.keinsci.com/forum.php ... mp;page=1#pid242360

作者
Author: thor 时间: 2023-11-20 21:24
是不是二面角的原因

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