受体蛋白-配体小分子分子动力学模拟基本流程--根据第14届培训内容整理

lijilo

写在前面：
2025年7月2日更新，补充部分常见分析方法数据输出及其相应命令

我参加了第14届北京科音分子动力学与GROMACS培训班，那可真是让我大开眼界，收获满满啊！卢天老师的讲解极为细致，真牛啊！我在这里真心建议有分子动力学需求的哥们儿去线下参加这个培训，绝对值，绝对不会亏！

备注：在内容之中的双斜杠//后的内容为注解，无双斜杠//的内容为运行程序，部分内容可能有刊误，请大家批评指正。

1.预处理受体蛋白、配体分子
//将受体蛋白另存为
protein.pdb
//将配体分子进行合适的质子化后另存为
mol.mol2

2.为小分子生成top文件和gro文件
打开sobtop，把mol.mol2拖入sobtpo
//选择1，生成top文件
//选择3，尽可能使用GAFF力场
//选择0，进入下一步
//选择4，如果可能，预先构建成键参数，任意猜测缺少的选项
//回车，生成的top文件生成在sobtop软件根目录下
//回车，生成的itp位置限制文件在sobtop软件根目录下
//选择2，生成gro文件
//回车，生成的gro文件在sobtop软件根目录下
//回车，退出sobtop软件

//将sobtop软件中的
mol.gro
mol.itp
mol.top
//三个文件剪切到实验工作环境目录下

3.产生蛋白质连带一个离子的拓扑文件
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top
//选择AMBER99SB-ILDN力场和TIP3P水(注意，如果静电荷不为零: Total charge in system x.000 e 则需要添加抗衡离子)
//得到的topol——Protein_chain_        A.itp对应蛋白，topol_Ion_chain_A2.itp对应离子,posre开头的文件对应限制势itp

4.合并gro文件，另存为complex.gro文件
//将mol.gro加入到pdb2gmx产生的protein.gro的末尾，并将第二行的原子数改为（蛋白质原子数+小分子原子数）,建议作图确认结构合理性，注意文件中的空格和回车问题
//注意，protein.gro文件末尾的三个小数是晶格的坐标，不要删除或大幅修改

//蛋白质的限制势itp文件在pdb2gmx的时候已经产生，但小分子的没有，genrestr是对输入的结构产生坐标或距离限制势itp文件的工具，接下来运行命令，进行限制势的产生
gmx genrestr -f mol.gro -o posre_mol.itp
//选择组的时候选择0，system默认的位置限制势常数是1000kJ/mol/nm2，已经足够大
//将下列语句插入到mol.itp文件的末尾，注意，复制时连“#”井号一同复制，最好在末尾添加之前空一行，方便检查文件错误

#ifdef POSRES
#include "posre_mol.itp"
#endif

//这样当mdp中使用define = -DPOSRES的时候配体的位置也会被限制了

//把配体的itp文件引入整体的拓扑文件topol.top，把分子数也设置好
//注意，在引入的时候需要将小分子的mol.itp文件引入到蛋白质链之前，因为mol.itp最开头定义了[atomtypes]因此，这个itp要最优先被引入
//将下列语句插入到引入蛋白质的itp文件引入之前；
#include "mol.itp"
并在末尾的[molecules]中引入mol 1，将topol.top的格式与complex.gro中分子出现的顺序对应
//即topol.top中应该是类似这样的顺序：

; Include forcefield parameters
#include "amber99sb-ildn.ff/forcefield.itp"

; Include chain topologies
#include "mol.itp"
#include "topol_Protein_chain_A.itp"
#include "topol_Ion_chain_A2.itp"

; Include water topology
#include "amber99sb-ildn.ff/tip3p.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

; Include topology for ions
#include "amber99sb-ildn.ff/ions.itp"

[ system ]
; Name
CATIONIC TRYPSIN in water

[ molecules ]
; Compound        #mols
Protein_chain_A     1
Ion_chain_A2        1
mol 1

//注意，末尾的mol和1之间有一个空格，且下一步添加溶剂后有可能会缺失回车，需要纠正

5.设定盒子，加水，加离子，能量极小化
//设定盒子
gmx editconf -f complex.gro -o complex_box.gro -d 0.8 -bt cubic
//设定的盒子是立方盒子，可能会增加计算量，但是不容易产生边界相互作用的问题

//加水
gmx solvate -cp complex_box.gro -o complex_SOL.gro -p topol.top
//注意这一步加水后有可能topol.top文件最后一行的SOL可能会串行，需要手动添加回车，避免其与mol 1连在同一行，容易在后续处理中报错

//将mdp模板中的em.mdp,pr.mdp,md.mdp拷贝到当前目录
//产生临时tpr文件
gmx grompp -f em.mdp -c complex_SOL.gro -p topol.top -o em.tpr -maxwarn 1
//这里如果警告较多可以将maxwarn的数值改大一些

//加离子，使得整个体系变为电中性（与第三步中的电荷数相对应）
gmx genion -s em.tpr -p topol.top -o system.gro -neutral
//这里选择分组时选择SOL对应的分组
//产生的带有离子且电中性的体系为system.gro

//能量极小化
gmx grompp -f em.mdp -c system.gro -p topol.top -o em.tpr

gmx mdrun -v -deffnm em

6.限制性动力学100ps
gmx grompp -f pr.mdp -c em.gro -p topol.top -r em.gro -o pr.tpr

gmx mdrun -v -deffnm pr
//注意，对于复杂体系如果用常规2fs步长可能模拟刚开始就会崩溃，因此此处做限制性动力学期间，步长用较小的1fs以求稳妥。


7.产生索引文件
gmx make_ndx -f pr.gro
依次输入：
x1|x2|x3（定义“蛋白-离子-配体”组，新的组号接在原有的组号之前，组号为R）（lijilo注释：|为右大括号右侧的按键，用“shift”+“、”便可输入该符号）
!R(把其他部分定义为一个组，新的组号为R+1)
name R protein_lig(改名为蛋白_配体)
name R+1 envir（改名为环境）
q
//得到的index。ndx里的组名就和模板文件里面的md.mdp中的组对应了

8.常规动力学1ns
gmx grompp -f md.mdp -c pr.gro -p topol.top -o md.tpr -n index.ndx -maxwarn 10
//如果这里警告多，可以将maxwarn后的数值改大一些

gmx mdrun -v -deffnm md
//由于配体、离子与蛋白的相互作用较为明显，所以这里任务之中将他们一起作为一个控温组和同一个消除屏东转动的组
//如果需要跑较长时间的动力学，需要修改md.mdp文件


9.分析RMSD
gmx rms -f md.xtc -s md.tpr -o rmsd_protein.xvg -n index.ndx
//两次都选骨架部分（Protein）
//注意观测在模拟过程中骨架的RMSD波动是否较大,计算蛋白质的RMSD值

gmx rms -s md.tpr -f md.xtc -o rmsd_complex.xvg -n index.ndx
//第一次选择骨架部分（Protein）第二次选择配体（MOL）
//可以消除蛋白质整体的运动，观察小分子配体相对于蛋白质的运动

10.计算RMSF
gmx rmsf -s md.tpr -f md.xtc -res -n index.ndx -o rmsf_protein.xvg
//选择骨架部分（Protein）

11.计算蛋白的回旋半径（Radius of Gyration）
gmx gyrate -s md.tpr -f md.xtc -o gyrate.xvg
//选择蛋白质的主链（Backbone）

12.计算溶剂可及表面积（SASA）
gmx sasa -s md.tpr -f md.xtc -o sasa.xvg -tu ns
//选择骨架部分（Protein）

13.计算氢键数量
gmx hbond -s md.tpr -f md.xtc -n index.ndx -num hbond.xvg -tu ns
//依次选择骨架部分（Protein）和配体部分（MOL）

14.生成分子动力学模拟轨迹动画
gmx trjconv -f md.xtc -s md.tpr -o movie_noWater.pdb -pbc mol -center -ur compact -dt 100 -n index.ndx
//选择骨架部分（Protein）和非水分子（non-Water）

sobereva

这里漏了给配体分子设置原子电荷的步骤。如果给sobtop的mol2文件里没有自带原子电荷，必须让sobtop获得原子电荷，常规且推荐的做法是让sobtop载入Multiwfn计算的记录了RESP电荷的chg文件，http://sobereva.com/soft/Sobtop里有具体说明。

楼主似乎假定了当前研究的蛋白是带着离子的，这实际上属于少数情况。如果蛋白质本身不带离子，也就不会出现topol_Ion_chain_A2.itp，[ molecules ]里也不会有Ion_chain_A2。

另外提醒一下，如果蛋白质和配体结合方式未知，一般适合先用分子对接程序得到对接后的复合物结构，比起自己在可视化程序里手动摆靠谱得多。

还有个建议是用AMBER14SB代替GROMACS自带的更老的AMBER99SB-ILDN。

lijilo

sobereva 发表于 2024-5-21 15:41
这里漏了给配体分子设置原子电荷的步骤。如果给sobtop的mol2文件里没有自带原子电荷，必须让sobtop获得原子 ...

感谢老师指正，这边接下来就优化一下后续gmx的实验流程

zsu007

谢谢分享！

低调的板凳

lijilo 发表于 2024-5-22 09:32
感谢老师指正，这边接下来就优化一下后续gmx的实验流程

sobtop用chg文件导入的话千万要注意，在正常导入mol2文件之后就先把gro文件生成了，之后再导入multiwfn的chg文件，生成top和itp文件。不然脚本生成的chg文件中的原子坐标会覆盖之前mol2的原子坐标，导致生成的gro文件中原子坐标错误。

MD-Chen

谢谢分享。

xxx1025

请问为什么限制性动力学模拟过程中为什么不需要分组控温呢？而是正式模拟时候才进行了分组控温。两者都是用V-rescale热浴，不会担心限制性动力学模拟过程中出现“热溶剂，冷溶质”的情况吗？

haizei

我不需要做配体受体复合物，而是将寡糖分子加入到蛋白质盒子中，然后我发现蛋白质的三级结构全变成了无序卷曲，这是为什么，拓扑文件也正确修改了

9heihei6

请问一下，如果跟着楼主的步骤进行操作，但是区别是我是用ACPYPE Server处理的配体，不过也是得到了ITP/TOP/GRO三个文件，但是在操作过程中，em时却提示警告了，请问可能是什么原因呢？我如何修改呢？

9heihei6

9heihei6 发表于 2025-2-21 01:01
请问一下，如果跟着楼主的步骤进行操作，但是区别是我是用ACPYPE Server处理的配体，不过也是得到了ITP/TOP ...

并且我在这一步忽视掉这个警告，却在进行em时又出现了报错，请问是同一个问题么？如何解决呢？
WARNING: Listed nonbonded interaction between particles 3475 and 3482
at distance 2.611 which is larger than the table limit 2.050 nm.

This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.

IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.

lijilo

9heihei6 发表于 2025-2-21 01:08
并且我在这一步忽视掉这个警告，却在进行em时又出现了报错，请问是同一个问题么？如何解决呢？
WARNING: ...

看上去你这里有两个原子之间的距离太近了呢，你可以找一下这两个原子，手动删除一个？

薄乐

//将mdp模板中的em.mdp,pr.mdp,md.mdp拷贝到当前目录
这三个文件应该在哪里找到呢？

lijilo

薄乐 发表于 2025-6-13 15:09
//将mdp模板中的em.mdp,pr.mdp,md.mdp拷贝到当前目录
这三个文件应该在哪里找到呢？

请看我帖子最后的附件

wangxiaohai

我是蛋白和2个配体，构建的topol文件一直有这样的报错，该怎么修改啊，我直接用的sobtop生成的小分子拓扑文件，求求了

wangxiaohai

用作者的这个命令“3.产生蛋白质连带一个离子的拓扑文件
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top
//选择AMBER99SB-ILDN力场和TIP3P水(注意，如果静电荷不为零: Total charge in system x.000 e 则需要添加抗衡离子)
//得到的topol——Protein_chain_        A.itp对应蛋白，topol_Ion_chain_A2.itp对应离子,posre开头的文件对应限制势itp”也没有对应的Protein_chain_A.itp产生，被卡1天了都