使用Sobtop超级方便地创建二茂铁的GROMACS的拓扑文件

sobereva

使用Sobtop超级方便地创建二茂铁的GROMACS的拓扑文件

Using Sobtop to create GROMACS topology file for ferrocene super conveniently

文/Sobereva@北京科音   2022-Feb-17

0 前言

Sobtop是笔者开发的主要用于产生各种类型体系的GROMACS拓扑文件的程序，可以在主页http://sobereva.com/soft/Sobtop免费下载。过渡金属配合物体系的拓扑文件在以往是很难产生的一类情况，因为里面牵扯到GAFF、OPLS-AA、GROMOS等主流有机分子力场不支持的元素，里面涉及到的键/键角/二面角参数在那些力场里也没有，所以acpype、Ligpargen、ATB等程序根本没法搞。而这类体系，用Sobtop则可以超级容易地产生拓扑文件。AmberTools里的MCPB.py程序虽然能产生这类体系的拓扑文件，但是过程麻烦至极，在http://bbs.keinsci.com/thread-27796-1-1.html里有人演示了怎么借助MCPB.py产生二茂铁的GROMACS的拓扑文件，估计大多数人看了后会打退堂鼓。Sobtop的详情以后会另文专门系统性地介绍，本文仅简单演示一下怎么用Sobtop快捷地产生二茂铁的拓扑文件，通过对比大家可以充分认识到用Sobtop有多么方便。

本文牵扯到的所有输入输出文件都可以在http://sobereva.com/attach/635/file.rar里找到。Sobtop用的是2022-Feb-16发布的预览版1.0(dev)。Multiwfn用的是2022-Feb-17更新的版本，Multiwfn可以在主页http://sobereva.com/multiwfn免费下载。


1 用Gaussian做优化和振动分析

首先用Gaussian对二茂铁做优化和振动分析。要注意，二茂铁的两个茂环间仅在低温晶体中是交错式的，而在真空中极小点结构是重叠式的（参看Materials, 8, 7723 (2015)和https://en.wikipedia.org/wiki/Ferrocene#Determining_the_structure中的相关信息），因此本例用重叠式（D5h点群）的初始结构。记得Gaussian计算前最好用GaussView做对称化，这样优化得又精准速度又快。

关于计算级别，本例用的TPSSh/def2-TZVP级别对于过渡金属配合物体系通常是很理想的选择。关于泛函的选择看《简谈量子化学计算中DFT泛函的选择》（http://sobereva.com/272）。如果想明显更便宜，Fe可以用SDD，配体可以用6-311G*。

此例的Gaussian输入文件如下

%chk=Ferrocene.chk
#p tpssh/def2tzvp opt freq
[空行]
Generated by Multiwfn
[空行]
0 1
C                  0.71056100   -0.97800331    1.65426528
C                 -0.00000000    1.20887857    1.65426528
C                 -0.71056100   -0.97800331    1.65426528
H                  2.17848562    0.70783289    1.62195576
H                 -2.17848562    0.70783289    1.62195576
Fe                 0.00000000    0.00000000   -0.00000000
C                  1.14971184    0.37356402    1.65426528
C                 -1.14971184    0.37356402    1.65426528
H                  1.34637816   -1.85313055    1.62195576
H                 -0.00000000    2.29059533    1.62195576
H                 -1.34637816   -1.85313055    1.62195576
C                  1.14971184    0.37356402   -1.65426528
C                 -1.14971184    0.37356402   -1.65426528
C                  0.71056100   -0.97800331   -1.65426528
C                  0.00000000    1.20887857   -1.65426528
C                 -0.71056100   -0.97800331   -1.65426528
H                  2.17848562    0.70783289   -1.62195576
H                 -2.17848562    0.70783289   -1.62195576
H                  1.34637816   -1.85313055   -1.62195576
H                  0.00000000    2.29059533   -1.62195576
H                 -1.34637816   -1.85313055   -1.62195576
[空行]
[空行]

算完了之后，用formchk把当前目录下的Ferrocene.chk转换为Ferrocene.fchk。不知道formchk是什么和怎么用的话看《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）里相关部分。

还需要得到二茂铁的mol2文件，并且要求其中记录的键连关系和期望的一致，这决定了Sobtop给出的拓扑文件里都有哪些力场项。很多程序都可以实现这个，此例我们用常用的GaussView打开Ferrocene.fchk，如果看到Fe和各个C之间没显示键的话就手动连上（设几重键都无所谓，Sobtop只看有没有键）。然后保存为Ferrocene.mol2。


2 计算RESP电荷

RESP电荷很适合用于动力学模拟的目的，我们要让拓扑文件里的原子电荷对应RESP电荷。计算RESP电荷非常强大、方便、灵活的程序是Multiwfn，见《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》（http://sobereva.com/441）。要注意对二茂铁体系不宜直接按照常规方式算RESP电荷，因为拟合点的分布不会恰好满足D5h对称性，因此得到的等价原子的原子电荷会有微小的偏差。为了让所有等价的原子的电荷都严格相同，在Multiwfn中计算RESP电荷时应当设置等价性。

启动Multiwfn，然后输入
Ferrocene.fchk
7  //计算原子电荷
18  //计算RESP电荷
5  //等价性约束设置
11  //根据点群自动判断等价性
a  //根据整个体系的结构判断点群
y  //当前屏幕上显示的点群确实是D5h，列出的原子等价性没错，所以这里输入y确认。然后当前目录下有了等价性设置文件eqvcons_PG.txt
q  //返回
1  //从文本文件里读取等价性设置
eqvcons_PG.txt  //刚才产生的文件
2  //开始一步式RESP电荷拟合
[回车]  //程序提示Fe缺半径，按回车表示用推荐的半径
Multiwfn算RESP电荷效率很高，RESP电荷很快就产生了，从屏幕上给出的结果看确实原子电荷合理，而且满足等价性。按y把产生的原子电荷导出为当前目录下的Ferrocene.chg。


3 用Sobtop产生拓扑文件

启动Sobtop，然后依次输入
Ferrocene.mol2
7  //指定原子电荷
10  //从Multiwfn的chg文件中读取原子电荷
Ferrocene.chg
0  //返回
2  //产生GROMACS的.gro文件
[回车]  //此时当前目录下产生了Ferrocene.gro
1  //产生GROMACS的拓扑文件
2  //优先用GAFF原子类型，缺的自动用UFF力场的。原子类型决定了拓扑文件里的LJ参数
2  //所有成键相关参数基于振动分析产生的fchk文件里的Hessian矩阵生成
Ferrocene.fchk  //如果此文件和Ferrocene.mol2在相同目录下，此处直接按回车可以快捷载入
[回车]  //在当前目录下产生Ferrocene.top
[回车]  //在当前目录下产生Ferrocene.itp

现在拓扑文件在当前目录下出现了，Sobtop可以关了。现在可以打开itp文件看看内容，会看到没有任何问题，文件整齐干净，注释也特别清楚。


4 测试拓扑文件合理性

产生拓扑文件之后，我一般建议在真空下对这个体系跑一下动力学，根据模拟过程中结构变化情况大致判断拓扑文件是否合理。做这个模拟要用的Ferrocene.gro、Ferrocene.top、Ferrocene.itp前面我们都已经产生了。做这个模拟用的md.mdp在本文的文件包里也给了，可见是控温在298.15K下模拟100 ps，1 fs一步，每100 fs往xtc里写入一帧，没有用PBC。笔者此例用的是GROMACS 2018.8。这个mdp不兼容GROMACS较新版本，较新版本用户需要在一个足够大的空盒子里在考虑PBC的情况下跑NVT模拟。

把前述的gro、itp、top、mdp都放到当前目录下，运行以下命令
gmx grompp -f md.mdp -c Ferrocene.gro -p Ferrocene.top -o md.tpr
gmx mdrun -v -deffnm md

然后用VMD载入轨迹，用Extensions - Analysis - RMSD Trajectory Tool工具对system选区做align消除平动转动后，每10帧叠加显示一次，得到下图

可见模拟过程中二茂铁很好地维持了刚性，说明拓扑文件合理。


5 总结&其它

本例充分体现了使用Sobtop结合Multiwfn程序产生过渡金属配合物体系拓扑文件的便利，把复杂度降低到了极限。此文只以极为简单的体系做了示例，对更复杂的情况，诸如多核金属配合物、除了刚性部分还带有可旋转键的柔性部分的配合物，Sobtop也都能很容易地产生拓扑文件。凡是体系全为刚性的情况都可以直接效仿本文的做法；对于既有刚性配位区域也有柔性基团的配合物，在Sobtop问你怎么产生成成键相关(bonded)参数时，建议选择相应选项，让bond和angle参数用特定方法基于Hessian产生，而其它的用预置力场参数，这样可以使得GAFF的参数用于描述柔性部分二面角可旋转特征。

实际中二茂铁的茂环的旋转势垒是非常低的，常温下是可以旋转的，但本文得到拓扑文件对应的是纯刚性二茂铁。本文暂不考虑茂环间可相对旋转的问题，这也不是什么重要问题。

本文的做法比起MCPB.py实在方便百倍，而且MCPB.py用的Seminario方法产生力常数也不如本文的方法合理（当前版本Sobtop默认用的是modified Seminario，产生的键角力常数远好于Seminario）。

做优化和振动分析用免费的ORCA亦可，文中用fchk的地方改成用ORCA振动分析产生的.hess文件即可。

对于体系没有对称性的情况，算RESP电荷时不需要本文这样做对称等价设置，直接在Multiwfn的RESP界面里选择1用常规的两步式拟合即可。更多细节看《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》（http://sobereva.com/441）。

本文的例子直接从chg文件中载入了RESP电荷。也可以在Sobtop里先不设置电荷，等itp产生出来之后再手动把要用的电荷写入到[ atoms ]部分。

使用Sobtop产生拓扑文件发表文章时记得按照Sobtop主页http://sobereva.com/soft/Sobtop的说明进行引用。引用方式在未来会有所变化，请记得发文章之前看一下主页上最新的引用说明。

sobereva

喵星大佬 发表于 2022-11-24 23:44
对于设计卤素的分子，需要用虚拟点描述卤键的情况，GAFF力场不涉及这样的操作，CGenFF有，是否考虑支持？

我不怎么喜欢CGenFF和CHARMM系的东西

或许目前尚没有文章提出在GAFF力场基础上考虑虚拟点描述卤键的做法，但原理上完全可以自行效仿CGenFF的做法去扩展，程序和力场方面本身不需要做什么变动。

喵星大佬

对于设计卤素的分子，需要用虚拟点描述卤键的情况，GAFF力场不涉及这样的操作，CGenFF有，是否考虑支持？

yjb

牧生 发表于 2022-11-22 18:31
看第六楼

谢谢

牧生

yjb 发表于 2022-11-22 17:57
老师，sobtop也可以像Multiwfn那样在Windows下添加环境变量，让他在任何文件夹下都可以使用吗？

看第六楼

yjb

老师，sobtop也可以像Multiwfn那样在Windows下添加环境变量，让他在任何文件夹下都可以使用吗？

sobereva

XiaoYan 发表于 2022-10-18 17:57
1.体系为带正电荷的LYS
高斯+antechamber连用计算得到frcmod文件，高斯关键词见截图。
2.根据高斯计 ...

没有叫sobtob的程序

sobtop里你选了什么都没完整体现

Gaussian使用犯了低级错误，理论方法和基组都没写，等于用的默认的HF/STO-3G，结果完全是垃圾

本来sobtop基于Hessian产生力常数的方法和antechamber也没直接可比性。sobtop用的什么在主页里都明确写明了。

XiaoYan

sobereva 发表于 2022-10-18 16:11
说清楚体系、说清楚怎么“计算”的

1.体系为带正电荷的LYS
高斯+antechamber连用计算得到frcmod文件，高斯关键词见截图。
2.根据高斯计算得到的fchk文件，通过multiwfn产生chg文件，见下图。利用chg和fchk文件通过sobtob产生所有原子的键长、键角和力常数参数。其中力常数与antechamber产生的参数差异巨大。

sobereva

XiaoYan 发表于 2022-10-18 09:43
Sob老师您好，我在使用Sobtop补充参数时发现用Sobtop计算出来的参数值与antechamber功能补充的参数值差异很 ...

说清楚体系、说清楚怎么“计算”的

XiaoYan

Sob老师您好，我在使用Sobtop补充参数时发现用Sobtop计算出来的参数值与antechamber功能补充的参数值差异很大，请问这是什么原因导致的呢？第一张图红色方框内是Sobtop计算的结果，第二张图是antechamber和parmchk连用得到的数值。

sobereva

naoki 发表于 2022-8-1 14:29
多谢您的指导，我正在根据您的帖子编译cp2k。还想请教一下您，这个COF晶体是个ABC错层形式，一个晶胞里有 ...

可以

naoki

sobereva 发表于 2022-8-1 02:24
可以
记得CP2K算COF的REPEAT电荷的时候只用一层，使得每个原子都能和上下的真空区接触。

多谢您的指导，我正在根据您的帖子编译cp2k。还想请教一下您，这个COF晶体是个ABC错层形式，一个晶胞里有三层COF（如下图所示），cif文件中每一层内的原子编号没有明显规律，但是整个晶胞的原子编号顺序是先第一层再第二层再第三层，这种情况下我是不是应该分别把每次层COF提出来加真空层算一遍REPEAT电荷，再按层顺序合并到完整晶胞的文件中呢，感谢~

sobereva

naoki 发表于 2022-8-1 00:28
谢谢您的回复！我去试一下您说的方法。我是想让把COF放在溶剂盒子中，在盒子xy平面上跨周期，z轴方向的两 ...

可以
记得CP2K算COF的REPEAT电荷的时候只用一层，使得每个原子都能和上下的真空区接触。

naoki

sobereva 发表于 2022-8-1 00:02
Multiwfn不支持QE产生的波函数。
对于周期性体系做动力学的目的（而且是原子和真空区有接触的情况，如CO ...

谢谢您的回复！我去试一下您说的方法。我是想让把COF放在溶剂盒子中，在盒子xy平面上跨周期，z轴方向的两侧与有机溶剂和离子接触，用gromacs跑动力学，这样的体系用REPEAT电荷合适吗，谢谢您~

sobereva

naoki 发表于 2022-7-31 22:48
sob老师好，谢谢您的回复。
还想请教一下您，我想算的是周期性COF晶体的电荷用于md，这个COF臂上带有吡 ...

Multiwfn不支持QE产生的波函数。
对于周期性体系做动力学的目的（而且是原子和真空区有接触的情况，如COF），原子电荷可以用CP2K算REPEAT电荷，下文提到了
使用Multiwfn非常便利地创建CP2K程序的输入文件
http://sobereva.com/587（http://bbs.keinsci.com/thread-21668-1-1.html）

只要mol2文件里记录的连接关系体现出了周期性，就能照常根据化学环境指认GAFF原子类型。

naoki

sobereva 发表于 2022-7-31 22:26
MMFF94原子电荷的计算从原理上就没法设总电荷。换句话说只能算中性体系的原子电荷。

sob老师好，谢谢您的回复。
还想请教一下您，我想算的是周期性COF晶体的电荷用于md，这个COF臂上带有吡啶阳离子结构，整体COF骨架荷正电（如下图a所示）。虽然我的mol2文件里已经定义好了键型，但一些基于成键关系获取电荷的程序都不能正确识别图a结构，而实识别为图b这种还原态结构（实际上mol2文件里吡啶环的成键方式为ar共振式，不是双键）。正确结构应该是左边a，比右边b少两个电子。



这么看来貌似只能通过量子化学的方法获取原子电荷了。对于这个COF，我目前打算用multiwfn载入QE的输出文件来计算它晶胞的MK电荷（跑md的时候我打算给COF增加位置限制势），然后再带着电荷扩胞用于md计算，不知这个做法是否正确可行，您觉得怎样的操作更为合适呢。还有一个问题就是对于这种骨架上有吡啶阳离子结构的周期性COF，sobtop能否正确识别gaff原子类型生成拓扑文件呢，谢谢您~