计算化学公社

标题: 使用Sobtop超级方便地创建二茂铁的GROMACS的拓扑文件 [打印本页]
作者
Author: sobereva    时间: 2022-2-17 12:27
标题: 使用Sobtop超级方便地创建二茂铁的GROMACS的拓扑文件
使用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帧叠加显示一次，得到下图

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

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


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的说明进行引用。引用方式在未来会有所变化，请记得发文章之前看一下主页上最新的引用说明。
作者
Author: lyj714    时间: 2022-2-17 13:35
本帖最后由 lyj714 于 2022-2-17 13:47 编辑

请问这个二面角的函数类型有什么讲究吗，我记得amber力场下一般用的是类型9，Sobtop产生的是2。这个应该是拟合方法不同而产生的差异吧，当然gromacs手册上也有说不同类型的函数形式。
作者
Author: sobereva    时间: 2022-2-17 14:02
lyj714 发表于 2022-2-17 13:35
请问这个二面角的函数类型有什么讲究吗，我记得amber力场下一般用的是类型9，acpype产生的貌似是3，Sobtop ...

关于二面角的functype
2：谐振势，一般用于实现improper项。Sobtop如果用基于Hessian产生力常数的谐振势描述某二面角项，也会用这个
4：单一普通周期势。AMBER/GAFF力场的improper项不是一般力场用的谐振势，而是普通周期势。4和1型二面角项定义完全一样，只不过为区分目的，4专门被用来记录AMBER/GAFF力场的improper项用。Sobtop基于GAFF力场产生拓扑文件时improper项用的是这个。
9：单一或多重普通周期势。一个二面角允许用多个周期性函数叠加描述。AMBER/GAFF力场的二面角项用的是这种，Sobtop基于GAFF产生拓扑文件时二面角项也用这种。
3：Ryckaert-Bellemans二面角。跟Sobtop和AMBER/GAFF力场都无关。
5：傅里叶形式二面角。跟Sobtop和AMBER/GAFF力场都无关。


早期版本（大概是4.5版以前），GROMACS不支持9形式的functype，即不允许AMBER/GAFF力场用的多个周期函数叠加描述一个扭转势的情况，因此早期acpype版本把这种扭转势等价地变换成了当时GROMACS可以用的RB形式（functype=3），这在acpype原文里也提到了。后来版本GROMACS有了functype=9了，比较新的acpype也直接用functype=9形式记录proper二面角项了。

作者
Author: lyj714    时间: 2022-2-17 14:14
sobereva 发表于 2022-2-17 14:02
关于二面角的functype
2：谐振势，一般用于实现improper项。Sobtop如果用基于Hessian产生力常数的谐振势 ...

原来如此，懂了，多谢老师指点。
作者
Author: 大村驴    时间: 2022-2-17 15:31
老师，建议像Multiwfn那样支持环境变量Sobtoppath，程序可以在当前目录、也可以在$Sobtoppath找ini、子程序和链接库文件。这样就可以将Sobtop存放在一固定目录，与用户文件互不影响。
作者
Author: sobereva    时间: 2022-2-17 15:52
大村驴 发表于 2022-2-17 15:31
老师，建议像Multiwfn那样支持环境变量Sobtoppath，程序可以在当前目录、也可以在$Sobtoppath找ini、子程序 ...

这不方实现，牵扯到零七八碎的东西
用户可以自行在.bashrc里通过alias定义一个快捷命令，一键进入Sobtop的目录
作者
Author: 大村驴    时间: 2022-2-17 16:53
老师，不是个重要的问题，当mol2的原子名的Fe写成大写的“FE”时候，会被程序认成F
作者
Author: exity    时间: 2022-2-17 17:13
佛以千万种形态在世间行走。给力！！！
作者
Author: 大村驴    时间: 2022-2-17 18:30
老师，Sobtop默认方法生成的Fc的拓扑在真空下没问题，在水溶液下会导致C3-H3严重向Fe弯折，不知是不是我的mdp有问题
作者
Author: sobereva    时间: 2022-2-17 19:48
大村驴 发表于 2022-2-17 16:53
老师，不是个重要的问题，当mol2的原子名的Fe写成大写的“FE”时候，会被程序认成F

我没测试出这个问题
作者
Author: 牧生    时间: 2022-2-17 20:48
本帖最后由 牧生 于 2022-2-17 23:01 编辑

看了一阵sobtop的操作流程，有几个小疑问

第一个小问题：
Choose the Hessian-based method used for deriving rigid parameters
1 Seminario
2 Modified Seminario (mSeminario, J. Chem. Theory Comput., 14, 274)
3 The second modified Seminario by Lu (m2Seminario)
4 Diagonal elements of redundant internal coordinate Hessian (DRIH)
在http://sobereva.com/soft/Sobtop/有这样的说法，DRIH方法对于许多刚性体系比默认的modified Seminario（另一种基于Hessian产生力常数的方法）更好，而用默认的modified Seminario对此例也完全没问题。  那么是不是意味着，对于常规的有机分子，如表面活性剂等，如果想要质量更好一点的，无脑就用4。


第二个小问题：
http://sobereva.com/soft/Sobtop/的几个例子中，分别用了“参数完全基于GAFF”，“力常数全基于Hessian计算得到”，“bond和angle参数基于Hessian得到，其它部分用GAFF”，是不是“力常数全基于Hessian计算得到”是质量最高的，因为这是从量化计算的结果而得到的参数。

第三个小问题：
以后会支持OPLS么

作者
Author: tjuptz    时间: 2022-2-17 21:12
卢老师把charmm排出主流了
作者
Author: sobereva    时间: 2022-2-18 06:01
牧生 发表于 2022-2-17 20:48
看了一阵sobtop的操作流程，有几个小疑问

第一个小问题：

1 否。以后专门写帖子的时候我会系统性说明这些差别
如果追求无脑，姑且用mSeminario，但个别时候也有缺陷，比如CLi4的不同键角k是明显不等价的，但用我提出的m2Seminario没这个问题。

2 没有可比性。因为很多体系必须考虑键的可旋转性，基于Hessian只能得到谐振势的k。

3 不会。我很讨厌OPLS的原子类型命名，一堆数字恶心死了，数目还很庞大，还有很多多余的原子类型。有GAFF就足够了，而且有AMBER、GLYCAM、Lipid、Merz等兼容的力场。我意识不到有什么明显理由应当支持OPLS-AA。

很多相关问题以后在Sobtop手册和专门的介绍文章里我都会详述。
作者
Author: sobereva    时间: 2022-2-18 06:03
tjuptz 发表于 2022-2-17 21:12
卢老师把charmm排出主流了

在力场方面我是AMBER一派的
CHARMM系列我不怎么喜欢
作者
Author: sobereva    时间: 2022-2-18 06:43
大村驴 发表于 2022-2-17 18:30
老师，Sobtop默认方法生成的Fc的拓扑在真空下没问题，在水溶液下会导致C3-H3严重向Fe弯折，不知是不是我的m ...

博文里的itp在博文中的文件包里直接给了

我实际测试了在水当中常温下的模拟，那种程度的弯折是正常的。用比如acpype产生苯的拓扑文件，在水盒子中常温模拟，照样会看到C-H的明显弯折。毕竟H的质量很小，在周围分子碰撞下，很容易运动，位置波动就是会很厉害。

如果就是想让弯折程度更低的话，Sobtop里选-2 Set scale factor for rigid constants，然后选5: Set scale factor of k for all terms，输入要给k乘的倍数，比如2、3。设得越大模拟中变形程度会越低。但也不能太大，否则模拟可能崩溃。
作者
Author: xiaowu759    时间: 2022-2-18 07:00
谢谢分享！回帖标记一下。
作者
Author: a9471163    时间: 2022-2-18 18:36
本帖最后由 a9471163 于 2022-2-18 19:15 编辑
sobereva 发表于 2022-2-18 06:01
1 否。以后专门写帖子的时候我会系统性说明这些差别
如果追求无脑，姑且用mSeminario，但个别时候也有缺 ...

老师您好，感谢您研发了这么好的工具，想请教一下：
1、我看到您本文中的示例是使用了“优先用GAFF原子类型，缺的自动用UFF力场的”。请问后续能否提供一个选项是“优先用GAFF原子类型，缺的自动用Merz力场的，如果Merz力场里也没有，再去搜索使用UFF力场的”。据我了解，Merz力场为一些金属离子优化过参数，这些金属离子的参数可能比UFF力场的更精确一些。

2、对于不常见元素，文献有时候会自己开发、测试并使用一套参数，请问是否支持把文献所用的参数手写到一个文本文件，然后用sobtop为特定的元素导入这套参数


3、我本人之前非常喜欢用GAFF力场，所有有机体系都是先考虑GAFF力场，后来在模拟一个氮、氧掺杂石墨烯时，发现GAFF力场对包含环氧基的六元环模拟失真（我个人能力有限，不清楚是acpype的问题还是GAFF力场本身的参数不够，拓扑是acpype产生后没修改直接用,现象在环氧基密度较高的区域才能体现的比较明显，所以当时这个模拟用的石墨烯是高度氧化的），这种失真体现在http://bbs.keinsci.com/thread-26918-1-1.html中13楼的贴图，即发现C原子所连的所有原子会有可能处于该C原子的同侧，造成环的不合理形变（测试也表明，这种现象跟RESP电荷无关，1）当不使用任何电荷时也有这种问题；2）在OPLS-AA力场下用一套不合理的RESP2（0.5）做过计算，当时报lincs warning的错误，模拟一段时间后崩溃，但所有的六元环都未观察到有GAFF力场下的这种形变的）。


测试发现在使用OPLS-AA力场时与C相连的所有原子不会都处于该C原子的同侧，而是必然有位于异侧的，所以是可以维持正常的C基面的形貌，UFF力场与OPLS-AA的结果一致。而CGenFF力场与GAFF力场的结果一致。有网友在帖子http://bbs.keinsci.com/thread-26918-1-1.html中的14楼和20楼留言认为这个问题可能是GAFF力场内置参数不够造成的，目前也只发现这一例，似乎是OPLS-AA模拟得更好，但也不确定探究所得的结论是否正确，想借此机会请教您一下。我非常赞同您说的GAFF力场相比于OPLS-AA的优势，以后能用GAFF的还是优先选用



作者
Author: sobereva    时间: 2022-2-19 09:28
a9471163 发表于 2022-2-18 18:36
老师您好，感谢您研发了这么好的工具，想请教一下：
1、我看到您本文中的示例是使用了“优先用GAFF原子 ...

1 正式版里Sobtop里会有个Merz的LJ参数文件，想用Merz的代替UFF的参数就把Merz的参数放到User_LJ.dat文件里。

并不能说Merz的LJ参数就比UFF的合理。Merz的参数甚至对于重现不同属性都有不同参数，再加上配合物中金属所处的化学环境各有不同，单纯对于描述配合物与其它分子的范德华相互作用来说，金属用重现特定属性来拟合的Merz的参数很可能还不如UFF的。再加上配合物里过渡金属原子一般就一两个，过渡金属的LJ参数对于模拟实际问题并不会有太大影响。

2 支持。自己把参数写进User_LJ.dat文件就完了

3 我喜欢GAFF关键在于GAFF的原子类型定义得清楚明白。至于成键相关参数，并非GAFF就比其它的更好。根据和量子化学算的振动频率对比，对于多数情况，bond和angle项改用modified Seminario基于体系实际的Hessian算的k比用GAFF的会明显更好。
作者
Author: naoki    时间: 2022-2-20 13:04
Sob老师竟然搞出了这种好东西
偶这OPLS党要转GAFF了
作者
Author: ggdh    时间: 2022-2-20 16:33
这个工具好像完全不依赖其他的工具(multiwfn除外)?
作者
Author: sobereva    时间: 2022-2-20 17:18
ggdh 发表于 2022-2-20 16:33
这个工具好像完全不依赖其他的工具(multiwfn除外)?

对
依赖得越多，用户用起来越困难

不过里面有个自动指认MMFF94原子电荷的功能，用户要用的话需要在ini文件里正确填入OpenBabel的可执行文件路径
作者
Author: sun877469558    时间: 2022-2-21 16:29
sob老师，此程序能否替代gmx中金属蛋白模拟时候所需要用到的MCPB.py？
作者
Author: 牧生    时间: 2022-2-21 20:02
本帖最后由 牧生 于 2022-2-21 20:03 编辑

刚发现sobtop可以在几秒内，基于GAFF，建立丙烯酰胺-丙烯酸线型共聚物的itp，但聚合度也就200而已，再提高聚合度到300就不行了。
在将来，sobtop能不能支持建立更高聚合度的线型聚合物呢？
作者
Author: sobereva    时间: 2022-2-21 20:08
牧生 发表于 2022-2-21 20:02
刚发现sobtop可以在几秒内，基于GAFF，建立丙烯酰胺-丙烯酸线型共聚物的itp，但聚合度也就200而已，再提高 ...

我不知道你具体遇到了什么情况，怎么不行
如果是调用atomtype识别原子类型的地方太慢，这没办法解决，只能用其它方法指认原子类型，比如从文件里直接读取，或者基于连接关系判断原子类型（正式版才有此功能）
作者
Author: sobereva    时间: 2022-2-21 20:09
sun877469558 发表于 2022-2-21 16:29
sob老师，此程序能否替代gmx中金属蛋白模拟时候所需要用到的MCPB.py？

MCPB.py不是直接用于gmx的
可能可以也可能不可以，目前还没专门考虑这种情况。就算直接不行，稍微扩展一下也能实现
作者
Author: 牧生    时间: 2022-2-21 20:54
本帖最后由 牧生 于 2022-2-21 20:55 编辑
sobereva 发表于 2022-2-21 20:08
我不知道你具体遇到了什么情况，怎么不行
如果是调用atomtype识别原子类型的地方太慢，这没办法解决，只 ...

已经发现问题所在。Running: atomtype -i F:\sobtop\sobtop_1.0(dev)\sobtop_1.0(dev)\polymer.mol2 -f mol2 -o MOL.ac -d ATOMTYPE_GFF.DEF
Info: The number of rings (4000) exceeded the default (1000).


因为分子链较长，在这个地方卡住了而已，我以为超过1000就不行了。。实际上，等十几秒 ~ 一分钟就好了。然后就能正常生成了。

世间竟有如此好物，建立线型聚合物不费吹灰之力，我等懒人福音。
作者
Author: sun877469558    时间: 2022-2-21 21:23
sobereva 发表于 2022-2-21 20:09
MCPB.py不是直接用于gmx的
可能可以也可能不可以，目前还没专门考虑这种情况。就算直接不行，稍微扩展一 ...

嗯嗯，对于gromacs用户来说，每次都要使用amber的MCPB.py再转换为gmx格式文件，如果sobtop能更简单更直接的gmx文件，那就很好了！期待。
作者
Author: sobereva    时间: 2022-3-3 15:58
jitou11 发表于 2022-3-3 11:27
我之前用过MCPB.py，我印象中是不可以的。MCPB.py的思路类似于QM/MM，它先是按与金属离子距离确定配体有 ...

我仔细看过MCPB.py的文章，没兴趣搞那么复杂。而且由于GROMACS没有AMBER那种leap工具，弄成MCPB.py的流程形式本来也不方便。
本身直接借助Sobtop算金属蛋白也不费事，自己挖个簇做opt freq，得到Hessian后用Sobtop自带的功能算出来需要的bond和angle的参数，然后在gromacs的拓扑文件里，在[intermolecular_interactions]里给金属和配位原子之间加入[bonds]和[angles]就完了。也没多费事，可控性还强，一个简单教程就能解释清楚。

Sobtop不会去支持AMBER，AMBER的拓扑文件乱死了。而且十几年前我从AMBER转投GROMACS后就不怎么用AMBER了。
作者
Author: sun877469558    时间: 2022-3-11 16:57
sobereva 发表于 2022-3-3 15:58
我仔细看过MCPB.py的文章，没兴趣搞那么复杂。而且由于GROMACS没有AMBER那种leap工具，弄成MCPB.py的流程 ...

如果sob老师有精力的话，还是请老师做一个算金属蛋白教程，因为对于gmx用户来说，用MCPB.py做实在是复杂，还要再转成gmx，增加了很多工作量
作者
Author: sobereva    时间: 2022-3-12 05:14
sun877469558 发表于 2022-3-11 16:57
如果sob老师有精力的话，还是请老师做一个算金属蛋白教程，因为对于gmx用户来说，用MCPB.py做实在是复杂 ...

最近没时间，有时间会写一个
作者
Author: sun666    时间: 2022-3-16 09:52
sobtop有一个重复的选项 如图4和6
作者
Author: Entropy.S.I    时间: 2022-3-16 22:50
sun666 发表于 2022-3-16 09:52
sobtop有一个重复的选项 如图4和6

两个“above”指代的内容不一样
作者
Author: nianbin    时间: 2022-3-20 10:17
这个简直太棒了
作者
Author: 竹会飞    时间: 2022-3-29 11:35
老师 您好我想问一下 ，我用sobtop 做环丁砜 的gromacs 的 itp文件 ，但是我发现 非键相互作用只给出一种 C 的sigma (nm)    epsilon (kJ/mol)   在原子电荷上是给出的两种C的电荷，我觉得从C的环境上 应该是两种C，这种情况应该怎么办？
作者
Author: heiheihaha    时间: 2022-3-30 12:02
Sob老师您好，请问我用Multiwfn处理金属配合物的RESP电荷是，提示Too many symmetry operations. Try a lower tolerance，但我设置更小的tolerance（0.09,0.08……）以后提示“Detected point group: C1 No symmetry-equivalence atoms were found”。请问这种情况怎么处理呢？谢谢~
作者
Author: sobereva    时间: 2022-4-1 01:29
竹会飞 发表于 2022-3-29 11:35
老师 您好我想问一下 ，我用sobtop 做环丁砜 的gromacs 的 itp文件 ，但是我发现 非键相互作用只给出一种 C ...

指认原子类型那一步仔细看屏幕上已指认的原子类型的解释，凭化学常识判断原子类型指认得对不对，不对就自己手动修改（根据LJ_param.dat里的原子类型的注释手动设为更合适的原子类型）
作者
Author: sobereva    时间: 2022-4-1 01:31
heiheihaha 发表于 2022-3-30 12:02
Sob老师您好，请问我用Multiwfn处理金属配合物的RESP电荷是，提示Too many symmetry operations. Try a low ...

如果目前官网上的最新版本Multiwfn也这样，把你的输入文件压缩后上传，说明具体操作步骤
作者
Author: youknowdcf    时间: 2022-5-11 09:50
感谢，有了这个再也没用过网站上的其他力场工具了。顺便请问社长，基于sobtop得到的配体是力场是GAFF的，与AMBER力场配套，在生物组分选择力场时选择amber99sb-ildn可以吗
作者
Author: sobereva    时间: 2022-5-11 20:13
youknowdcf 发表于 2022-5-11 09:50
感谢，有了这个再也没用过网站上的其他力场工具了。顺便请问社长，基于sobtop得到的配体是力场是GAFF的，与 ...

可以
不过这年头最好用AMBER14SB或AMER19SB
作者
Author: youknowdcf    时间: 2022-5-12 09:05
sobereva 发表于 2022-5-11 20:13
可以
不过这年头最好用AMBER14SB或AMER19SB

感谢社长，刚看了http://bbs.keinsci.com/thread-15094-1-1.html里的链接地址，但发现gmx的页面改版了，找不到user-contributions页面了，也找不到amber14sb和19sb了，想厚颜向社长讨要一下
作者
Author: 阿di达思    时间: 2022-7-30 19:39
感谢社长，Sobtop后期版本会将INTERFACEFF参数写进User_LJ.dat文件吗？研究材料界面作用的用户应该很多。
作者
Author: sobereva    时间: 2022-7-31 03:40
阿di达思 发表于 2022-7-30 19:39
感谢社长，Sobtop后期版本会将INTERFACEFF参数写进User_LJ.dat文件吗？研究材料界面作用的用户应该很多。

现在没明确打算，也没时间弄
如果有谁有兴趣写一份我也欢迎
作者
Author: naoki    时间: 2022-7-31 16:03
本帖最后由 naoki 于 2022-7-31 17:08 编辑

sob老师好，我想请教您一下sobtop能否在调用openbabel生成MMFF94电荷的时候设置体系总电荷为某一定值，谢谢~
好像openbabel的MMFF94插件本身就没有设定net charge的地方

作者
Author: sobereva    时间: 2022-7-31 22:26
naoki 发表于 2022-7-31 16:03
sob老师好，我想请教您一下sobtop能否在调用openbabel生成MMFF94电荷的时候设置体系总电荷为某一定值，谢谢 ...

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

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

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

这么看来貌似只能通过量子化学的方法获取原子电荷了。对于这个COF，我目前打算用multiwfn载入QE的输出文件来计算它晶胞的MK电荷（跑md的时候我打算给COF增加位置限制势），然后再带着电荷扩胞用于md计算，不知这个做法是否正确可行，您觉得怎样的操作更为合适呢。还有一个问题就是对于这种骨架上有吡啶阳离子结构的周期性COF，sobtop能否正确识别gaff原子类型生成拓扑文件呢，谢谢您~
作者
Author: sobereva    时间: 2022-8-1 00:02
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原子类型。
作者
Author: naoki    时间: 2022-8-1 00:28
sobereva 发表于 2022-8-1 00:02
Multiwfn不支持QE产生的波函数。
对于周期性体系做动力学的目的（而且是原子和真空区有接触的情况，如CO ...

谢谢您的回复！我去试一下您说的方法。我是想让把COF放在溶剂盒子中，在盒子xy平面上跨周期，z轴方向的两侧与有机溶剂和离子接触，用gromacs跑动力学，这样的体系用REPEAT电荷合适吗，谢谢您~
作者
Author: sobereva    时间: 2022-8-1 02:24
naoki 发表于 2022-8-1 00:28
谢谢您的回复！我去试一下您说的方法。我是想让把COF放在溶剂盒子中，在盒子xy平面上跨周期，z轴方向的两 ...

可以
记得CP2K算COF的REPEAT电荷的时候只用一层，使得每个原子都能和上下的真空区接触。
作者
Author: naoki    时间: 2022-8-1 14:29
本帖最后由 naoki 于 2022-8-1 14:31 编辑
sobereva 发表于 2022-8-1 02:24
可以
记得CP2K算COF的REPEAT电荷的时候只用一层，使得每个原子都能和上下的真空区接触。

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

(, 下载次数 Times of downloads: 47)
作者
Author: sobereva    时间: 2022-8-1 22:06
naoki 发表于 2022-8-1 14:29
多谢您的指导，我正在根据您的帖子编译cp2k。还想请教一下您，这个COF晶体是个ABC错层形式，一个晶胞里有 ...

可以
作者
Author: XiaoYan    时间: 2022-10-18 09:43
Sob老师您好，我在使用Sobtop补充参数时发现用Sobtop计算出来的参数值与antechamber功能补充的参数值差异很大，请问这是什么原因导致的呢？第一张图红色方框内是Sobtop计算的结果，第二张图是antechamber和parmchk连用得到的数值。
作者
Author: sobereva    时间: 2022-10-18 16:11
XiaoYan 发表于 2022-10-18 09:43
Sob老师您好，我在使用Sobtop补充参数时发现用Sobtop计算出来的参数值与antechamber功能补充的参数值差异很 ...

说清楚体系、说清楚怎么“计算”的
作者
Author: XiaoYan    时间: 2022-10-18 17:57
sobereva 发表于 2022-10-18 16:11
说清楚体系、说清楚怎么“计算”的

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

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

没有叫sobtob的程序

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

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

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

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

看第六楼
作者
Author: yjb    时间: 2022-11-23 09:27
牧生 发表于 2022-11-22 18:31
看第六楼

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

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

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



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3