计算化学公社

标题: 求教：使用sobtop生成MOF的两个itp文件，一个做能量最小化失败，另一个正常 [打印本页]

作者
Author: a9471163 时间: 2023-4-3 17:53
标题: 求教：使用sobtop生成MOF的两个itp文件，一个做能量最小化失败，另一个正常
老师您好，我做一个已发表MOF的分子动力学模拟，CCDC：924088，名称：TKL-107。初始结构下载于CCDC官网。
这同一个MOF的CIF载入GaussView后，经检查，成键无误。扩胞为2*2*2，随后GaussView分别保存成PDB和mol2文件（mol2已上传），

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

用 gmx editconf -f 924088.pdb -o 924088-box.pdb -box 3.03476 3.03476 3.7718 -angles 90 90 120给PDB加上盒子（该PDB结构文件已上传，晶胞参数是使用了文献值，这也是能量最小化时用的初始结构文件），

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

然后把带盒子信息的PDB提交给sobtop，依次键入1,1,0,4 (Prebuilt bonded parameters if possible, but missing ones are arbitrarily guessed），然后再输入2次回车就产生了ITP文件，另外mol2也提交给sobtop，用完全相同的操作生成另一个ITP文件，

随后发现，当用mol2生成的ITP文件时，能量最小化不成功，最大的力有10的5次方（图1）；当用PDB生成的ITP文件时，能量最小化正常进行到力小于100（图2），想请教一下为什么使用sobtop时，会有这种区别呢？

(, 下载次数 Times of downloads: 1)
图1：能量最小化不成功

(, 下载次数 Times of downloads: 2)
图2：能量最小化快速成功

个人推测：都没有用原子电荷，ITP中的力常数完全一样，MOF初始结构完全一样，力场都是UFF也一样，能量最小化都是只有MOF在盒子里，没加其他分子。这些变量都相同后，我目前发现的唯一不同是：对比两份ITP文件发现，平衡键长、键角有微小区别，但也只是小数点后第二位有极小的区别（如图3，箭头所示的是平衡键角, 其中924088-box.itp是带盒子PDB生成的，而924088.itp是mol2生成的），

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

图3：同一个化学键的平衡键角值有微小区别

所以，会不会是mol2产生的平衡键长、键角不准呢，如果每次都用量化计算纠正这些不准的平衡键长、键角，可能非常耗时、昂贵，不知还有没有其他方法。还请老师指正~

我喜欢用GAFF力场，但我试了，用mol2给sobtop生成GAFF（除金属外的其他元素）+UFF力场（MOF中的金属元素）参数是无法做能量最小化的，遇到的问题跟上面用UFF力场时一样。

最后，我上传了压缩包，包括可以正常模拟的924088-box.itp。不能正常模拟的924088.itp，我用的mdp文件，以及topol.top。以供各位感兴趣的老师测试、交流。其中UFF原子类型在top文件里，其他参数在itp文件里了。

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

作者
Author: sobereva 时间: 2023-4-4 08:02
差异不仅是那点微小的差别。itp文件的总行数都不同。而且[ angles ]部分2 3 366键角都明显不同。没具体结构文件没法说。
注意对比载入pdb和mol2文件时原子间连接关系判断的是否一致

作者
Author: a9471163 时间: 2023-4-4 20:00
本帖最后由 a9471163 于 2023-4-4 20:05 编辑

sobereva 发表于 2023-4-4 08:02
差异不仅是那点微小的差别。itp文件的总行数都不同。而且[ angles ]部分2 3 366键角都明显不同。 ...

谢谢老师提醒，具体结构文件——924088.cif来自于文献Fluorous metal-organic frameworks with enhanced stability and high H2/CO2 storage capacities，我上传了，
(, 下载次数 Times of downloads: 1)

itp文件的总行数不同，确实如此，目前我也没看出来为什么，有结果后跟您反馈。

[ angles ]部分2 3 366键角明显不同，这是一个跨盒子的角，目前已经证明，虽然两份itp中该键角不同，但二者之和为180度，也就是互为补角。

对比载入pdb和mol2文件时原子间连接关系判断的是否一致：GaussView载入观察后觉得是一致的（不过mol2有很多跨盒子的键，存在阻挡视线不方便观察的可能性），因为pdb和mol2文件都是从我这次上传的924088.cif通过同一个软件——GaussView另存为得到的，另存之前我用GaussView观察晶胞，那些成键显示都是正常的（我这里贴了图），所以理论上另存之后它们原子间连接关系应该也是一样的，GaussView总不能保存两次文件得到两个不同的连接关系。

最后，我还在用cp2k优化晶胞的键长、键角，等有结果后拿这些量化得到的键长、键角生成mol2，随后mol2提交给sobtop生成拓扑试试，有结果后跟您反馈。

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

作者
Author: sobereva 时间: 2023-4-5 06:41

a9471163 发表于 2023-4-4 20:00
谢谢老师提醒，具体结构文件——924088.cif来自于文献Fluorous metal-organic frameworks with enhanced ...

要看sobtop载入结构时显示的连接关系
sobtop对pdb文件是基于原子半径和距离判断成键。虽然gview默认也这么判断，但用的阈值也不同。所以用gview来对比没意义

作者
Author: a9471163 时间: 2023-4-9 11:37

sobereva 发表于 2023-4-5 06:41
要看sobtop载入结构时显示的连接关系
sobtop对pdb文件是基于原子半径和距离判断成键。虽然gview默认也这 ...

(, 下载次数 Times of downloads: 23)
(, 下载次数 Times of downloads: 6)
谢谢老师指点，问题已彻底解决，方法是需要手动把GaussView导出的mol2加上晶胞盒子信息，再把mol2给sobtop生成力场文件，力场用GAFF+UFF完全没问题了。成功的模拟文件已经打包上传。之前没有盒子，可能是GaussView导出的mol2在sobtop里判断成键出了些问题

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