gmx跑蛋白-小分子的md时是否需要将小分子用量化计算后再跑？

ezez

准备跑蛋白-小分子化合物的md，从rcsb下载的复合物文件（其中包含蛋白-化合物小分子.PDB），不确定是否需要在ligand的准备阶段用量化的方法把化合物优化一下。
之前看过一个nature chemical biology的文章，将ligand和蛋白质分别作了以下处理后再跑md：
1. 蛋白：用在线工具建模，补全了蛋白的缺失的序列。
2. ligand：用g09跑了优化，（估计是几何优化）。
3. 把ligand和蛋白进行blind dock，最后跑md

原文描述如下：
The chemical structures of MCC950 and ADP used in the BD and MD calculations were built up and optimized fully at the B3LYP/6–31+G(d,p) level within the density functional theory (DFT) framework, implemented in Gaussian09 (ref. 32). Additional vibrational calculations were conducted at the same DFT level to verify the nature of the located structures in the potential energy surface. Partial atomic charges were subsequently computed with the electrostatic potential scheme to be used during BD simulations.


ref: Nat Chem Biol 2019 Jun;15(6):560-564.


我的理解：
从rcsb拿到的复合物的结构不需要在ligand准备阶段单独对ligand跑量化，因为晶体结构是现成的，如果把ligand的结构进行优化（做量化的结合优化或者用molclus寻找能量最低的构象后），可能反而破坏了晶体结构原本的构象。
而上面NCB的文章需要对ligand跑量化进行几何结构优化，是因为作者当时并没有复合物的晶体结构，需要将ligand结构优化后再去blind dock，这一点和我已经有了复合物的晶体结构是不同的。

社长在gmx培训班上的举例子时也说用gview把配体结构打开，加氢等，未提及要跑量化来找配体的最优构象。

如果实在想优化ligand的结构，也可以只用社长的molclus寻找最低能量的构象，然后用该构象替换rcsb的复合物的晶体结构中的ligand，再跑md即可～

ps：我们单位没买g09或g16，也没法跑量化～ 听社长在培训班上说过orca不错，但是没用过orca，估计几何优化也不会跑～～

enthalpy

对于你提到的问题。做量化的计算，很多情况下不就是为了后面拟合电荷等参数。如果有配体已经有参数(来自文献或是力场包里就有)，你就顶多给小分子加氢，不用做量化计算优化结构了。

sobereva

对于算蛋白-配体复合物目的，若要优化配体，用处也仅仅是为了能够之后在一个合理结构下基于波函数算RESP、ADCH等原子电荷用于分子动力学
如果要用的原子电荷是是纯粹依赖于连接关系经验方式得到的，诸如MMFF94，完全不需要优化配体（当然，一般也不会把MMFF94电荷用于这种目的）
往届的北京科音分子动力学与GROMACS培训班上我讲用acpype的时候，优化是acpype脚本调用Ambertools程序的SQM模块在半经验级别下自动做的（除非用我课上给的noopt的版本刻意略过优化步骤）

对于要经过几何优化然后算RESP电荷的情况，基于已知的复合物中配体的初始构象，加氢后直接优化就够了，也没必要做往往很昂贵的构象搜索。本来配体在蛋白口袋里结合的构象往往就不是能量最低构象，刻意用能量最低构象算RESP电荷还不如用以复合物里的结构作为初猜直接优化的构象。

ezez

enthalpy 发表于 2022-8-2 21:00
对于你提到的问题。做量化的计算，很多情况下不就是为了后面拟合电荷等参数。如果有配体已经有参数(来自文 ...

感谢～

我先试试，后面有体会了再来讨论～～

ezez

sobereva 发表于 2022-8-3 04:50
对于算蛋白-配体复合物目的，若要优化配体，用处也仅仅是为了能够之后在一个合理结构下基于波函数算RESP、A ...

“本来配体在蛋白口袋里结合的构象往往就不是能量最低构象”
谢谢社长提醒，我之前一直以为配体是最低能量构象，看来并不一定如此啊～～

不过想问一下，这个是有什么依据吗？rcsb里的晶体应该不是最低能量构象，这一点好理解（因为蛋白复合物在结晶时也是需要特殊条件的），但是小分子的构象这么多，真正和蛋白结合的究竟是哪个？只能跑完md得到最低势能的那一帧后dump才能获得最低最稳定的构象吗～～

sobereva

ezez 发表于 2022-8-3 16:33
“本来配体在蛋白口袋里结合的构象往往就不是能量最低构象”
谢谢社长提醒，我之前一直以为配体是最低能 ...

本来很多配体分子就是高度柔性的，配体和蛋白质残基存在相互作用，显然会影响配体的构象
复合物处于能量最低状态，显然也不意味着其中某个部分，诸如配体，处于能量最低状态

ezez

sobereva 发表于 2022-8-3 04:50
对于算蛋白-配体复合物目的，若要优化配体，用处也仅仅是为了能够之后在一个合理结构下基于波函数算RESP、A ...

社长，看了您的回复，我的理解是：

1. 对于本例中从rcsb中获得蛋白+化合物晶体结构的情况，小分子的结构文件还是直接用rcsb里获得的化合物结构-加氢-去准备md的文件即可。
2. 可以直接用已经从rcsb里获得的配体结构加氢，再用orca进行opt+frq，然后算resp2电荷，用于跑md---该步骤是为了获得比较严谨的resp2电荷。

3. 对于没有获得蛋白+化合物的复合物晶体的情况，则用：
1）xtb在GFN0-xTB下做动力学模拟
2）Molclus调用xtb在GFN0-xTB下做批量优化
3）用Molclus调用xtb在GFN2-xTB结合隐式水模型下做批量优化
4）Molclus调用ORCA得到水环境下的高精度自由能（orca 算opt+frq，再用orca算高精度单点能，再用isostat获得最佳化合物构象）
5）以上一步中的最佳构象和高精度单点能来准备md文件，以及用orca的frq任务产生的.hess文件通过Sobtop生成.gro及.top文件，跑md。
对吗？
ps：感觉如果想让orca产生hess文件，不能通过molclus调用orca，应该单独使用orca来计算，否则molclus目录下没找到hess文件～（是否可以通过设置setting.ini文件中的ibkhess = 1来实现？）

sobereva

ezez 发表于 2022-10-1 15:52
社长，看了您的回复，我的理解是：

1. 对于本例中从rcsb中获得蛋白+化合物晶体结构的情况，小分子的结 ...

是freq不是frq

关于算配体的RESP电荷用的构象问题，看下文文末的说明
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441（http://bbs.keinsci.com/thread-10880-1-1.html）

构象搜索没必要非得用xtb跑动力学的方式，可旋转的键不是很多的情况下用gentor更方便，而且确保采样全面
gentor：扫描方式做分子构象搜索的便捷工具
http://bbs.keinsci.com/thread-2388-1-1.html

可以用ibkhess

一般有机分子直接让sobtop产生纯GAFF参数的拓扑文件就行了，没必要基于Hessian产生谐振势的参数，看“FAQ 11：什么时候才需要Sobtop基于Hessian算力常数？”（http://sobereva.com/soft/Sobtop/#FAQ11）

jim

sobereva 发表于 2022-10-2 06:32
是freq不是frq

关于算配体的RESP电荷用的构象问题，看下文文末的说明

请问社长，同一个配体，利用autodock对接可能产生不同坐标的pdb文件，请问：不同坐标的配体应该只算一次RESP2电荷就可以了吧？不用每个不同坐标的pdb文件都算一次电荷吧？谢谢
---------------------------------分割线补充
联系到计算RESP2中用到的mol2文件（包含坐标信息），以及后期MD所需的配体.gro文件里的坐标信息，还是每个不同坐标的pdb文档算一次RESP2电荷为好吧（除非乐于更改坐标信息）？

sobereva

jim 发表于 2024-1-3 22:37
请问社长，同一个配体，利用autodock对接可能产生不同坐标的pdb文件，请问：不同坐标的配体应该只算一次R ...

通常来说只对最重要的构象算（或者Multiwfn算RESP电荷时考虑构象权重平均）即可。但假定有不同比较稳定的结合方式、MD过程中小分子构象也不怎么变，那么分别对不同构象的小分子拟合RESP电荷做相应的MD，原理上还会更好一些。

更多讨论看此文文末
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441（http://bbs.keinsci.com/thread-10880-1-1.html）