眉山市网站建设_网站建设公司_营销型网站_seo优化

在蛋白质设计领域，我们不仅关心如何“创造”一个新的蛋白质结构或相互作用界面，更关心我们的设计在真实的、动态的、水环境下的表现如何。一个在静态结构上看似完美的设计，可能在生理环境中因动态柔性、溶剂效应或熵变等因素而失去预期功能。这时，分子动力学模拟便成为连接“静态设计蓝图”与“动态功能实景”不可或缺的桥梁。

然而，针对“小分子-蛋白质复合物”这一高频应用场景，网络上相关的模拟教程往往步骤割裂、版本混杂，让初学者，尤其是专注于蛋白质设计的研究者，在实践时屡屡受挫，难以将模拟真正转化为设计迭代的有效工具。在历经诸多“段错误”与未收敛的警告后，我也终于梳理出一条清晰、可复现的路径。

今天，我很荣幸能将这份融合了多篇专家文献精髓与亲身实战经验的完整操作指南分享给蛋白质设计社区的各位同仁。本指南旨在将分子动力学模拟流程标准化，帮助设计者系统性地评估：

1. 设计蛋白与配体结合的动态稳定性（RMSD, RMSF）。

2. 关键相互作用网络（如氢键）在模拟中的持续性。

3. 结合口袋在溶剂环境下的构象变化与柔性。希望它能为大家优化结合亲和力、验证突变效果、理解构效关系提供一套可靠的计算显微镜。

首要的感谢，必须献给sobereva老师（卢天）。他创建的Multiwfn、Sobtop等强大工具，以及二十年来在计算化学领域海量、无私的答疑解惑，为我们所有后来者照亮了前路。这份指南中关键的电荷计算与拓扑生成步骤，便深深依赖于他的工作。

同样诚挚的感谢，给予Computational Chemistry Commune论坛的lijilo作者，其详实的GROMACS操作帖是本文重要的学习蓝本。（http://bbs.keinsci.com/thread-45791-1-1.html）

我坚信，开源、共享与规范的传承是科学进步的基石。因此，在使用本指南涉及的软件与方法时，请务必遵守学术规范，正确引用相关作者的工作（文中已用红色高亮标出）。任何无视此提醒的行为，皆是对无数先行者心血的不尊重，也将损害整个社区的健康生态。

道阻且长，行则将至。愿这份指南能助各位在蛋白质设计的动态世界中，更从容地验证、优化与创造。

（由于个人编写，其中难免笔误，还请不吝赐教，还望海涵。）

需要提前安装的软件：

1.Gromacs（https://manual.gromacs.org/）

2.Multiwfn（http://sobereva.com/multiwfn/，最近好像崩溃了，可以过几天再下载）

3.sobtop（http://sobereva.com/soft/Sobtop/）

4.VMD（https://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD）

5.Pymol（https://pymol.org/）

6.Gaussian（正版需付费，自行寻找方法下载）

目录

1.分子对接（可选）

2 使用multiwfn和Gaussian为小分子设置原子电荷

3.为小分子生成top文件和gro文件

4.产生蛋白质的拓扑文件

5 合并gro文件，另存为complex.gro文件

6 生成小分子限制势

7 编辑top文件

8 设定盒子，加水，加离子，能量极小化

9 进行NVT平衡

10 进行NPT平衡

11 产生索引文件

12 进行分子动力学

13 计算分子动力学结果参数

14 生成分子动力学模拟pdb文件

1.分子对接（可选）

1.1 将蛋白与小分子放入cbdock2（http://183.56.231.194:8001/cb-dock2/index.php）中，进行分子对接。

1.2 下载结合位点合适的小分子-蛋白复合物 。

1.3 使用pymol打开protein_ligand_out_5.-8.8.complex.pdb，删除

小分子。

1.4 打开顶部下拉菜单选择“File > Export Molecule...”。

1.5 保存为protein.pdb。

1.6 再次使用pymol打开protein_ligand_out_5.-8.8.complex.pdb，在pymol中删除蛋白。

1.7 小分子加氢

1.8 另存小分子为mol.mol2。

1.9 创建一个工作文件夹（如work2），将protein.pdb，mol.mol2复制到文件夹。

注：如果使用cbdock2，请引用以下文章

Yang Liu, et al, CB-Dock2: improved protein-ligand blind docking by integrating cavity detection, docking and homologous template fitting. Nucleic Acids Research, 2022. Xiaocong Yang, et al, FitDock: protein-ligand docking by template fitting. Briefings In Bioinformatics, 2022.

不然后果自负

2 使用multiwfn和Gaussian为小分子设置原子电荷

2.1 将multiwfn的Multiwfn_3.8_dev_bin_Linux/examples/RESP/RESP2.sh文件复制到工作文件夹。

chmod _x ./RESP2.sh ./RESP2.sh mol.mol2

注1：RESP2中默认用的是Gaussian的g09版本，如果安装的是g16版本，请打开RESP2.sh文件将此处g09改为g16。

注2：默认运行较慢，可以在RESP2.sh的#Optimize geometry的$chk=gau.chk后加上使用内存和核心数，一般内存使用当前可用内存的2/3，线程则使用全部当前可用线程。

添加内容如下，其中%mem是内存最高运行用量，%nprocs是运行所用线程数，一般可以使用全部线程。

%mem=4GB %nprocs=12

ps：此步运行时间较长，作者运行了约30分钟（如果运行失败，可检查小分子是否加氢，见1.10步）

2.2 运行结束后会生成3个chg文件

注：如果使用multiwfn，请引用以下文章

Tian Lu*, Feiwu Chen, Multiwfn: A multifunctional wavefunction analyzer, J. Comput. Chem., 33, 580-592 (2012)

Tian Lu, A comprehensive electron wavefunction analysis toolbox for chemists, Multiwfn, J. Chem. Phys., 161, 082503 (2024)

不然后果自负

3.为小分子生成top文件和gro文件

3.1 进入事先下载，解压完毕的sobtop文件夹下

cd xxx/xxx/xxx/sobtop_1.0(dev5)

3.2 将文件夹下所有文件变为可执行文件

chmod +x *

3.3 运行sobtop

./sobtop

3.4 输入小分子文件

xxx/xxx/xxx/work2/mol.mol2

3.5 生成小分子top文件和gro文件

//选择7，分配原子电荷

//选择10，读取multiwfn文件输出文件

//输入mol.chg文件所在位置

//选择0，返回

//选择1，生成top文件

//选择3，尽可能使用GAFF力场

//选择0，进入下一步

//选择4，如果可能，预先构建成键参数，任意猜测缺少的选项

//回车，生成的top文件生成在sobtop软件根目录下

//回车，生成的itp位置限制文件在sobtop软件根目录下

//选择2，生成gro文件

//回车，生成的gro文件在sobtop软件根目录下

//回车，退出sobtop软件

3.6 复制sobtop文件夹下的mol.top文件，mol.tip文件，mol.fro文件到w

注：如果使用Sobtop，请引用以下文章

Tian Lu, Sobtop, Version [当前版本], http://sobereva.com/soft/Sobtop (accessed on 日 月 年)

Sobtop在进一步完善后会写详细的手册并发布正式版，然后专门写成介绍paper发在学术期刊上，届时请引用期刊上的文章。

不然后果自负

4.产生蛋白质的拓扑文件

注：回到工作文件夹（我的是work2）下

4.1 使用pdb2gmx产生蛋白质的top文件

gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top

4.2 输入AMBER99SB-ILDN protein对应的数字，回车。

4.3 输入TIP3P对应的数字，回车。

4.4 这一步最终生成三个文件（如果有多条链，会生成多个文件）

或

5 合并gro文件，另存为complex.gro文件

5.1 复制protein.gro为complex.gro

cp protein.gro complex.gro

5.2 将mol.gro的文件全部内容复制到complex.gro最后一行

5.3 删除掉多余内容

删除掉这些，不要有多余空行

mol 44

（这一步不知道是否正确，但不删掉有时会报错）

5.4 将9.50704 8.78086 7.86566移至最后

5.5 将complex.gro第二行数值改为原数值加上mol.gro中第二行数值。

12898+44=12942

6 生成小分子限制势

6.1 蛋白质的限制势itp文件在pdb2gmx的时候已经产生，但小分子的没有，genrestr是对输入的结构产生坐标或距离限制势itp文件的工具，接下来运行命令，进行限制势的产生

gmx genrestr -f mol.gro -o posre_mol.itp

6.2 选择组的时候选择0，system默认的位置限制势常数是1000kJ/mol/nm2，已经足够大，结果产生posre_mol.itp

6.3将下列语句插入到mol.itp文件的末尾，注意，复制时连“#”号一同复制，最好在末尾添加之前空一行，方便检查文件错误。

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

添加后：

6.4 保存，关闭

7 编辑top文件

7.1 打开topol.top文件

7.2 在文件首添加#include "mol.itp"

#include "mol.itp"

7.3在文件尾红框处添加mol 1

其实一个空格就可以，如果有强迫症可以把1调整到上面一样的位置

7.4 保存，关闭

8 设定盒子，加水，加离子，能量极小化

8.1 设定的盒子是立方盒子，可能会增加计算量，但是不容易产生边界相互作用的问题

gmx editconf -f complex.gro -o complex_box.gro -d 0.8 -bt cubic

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

gmx solvate -cp complex_box.gro -o complex_SOL.gro -p topol.top

8.3 将mdp模板中的em.mdp,pr.mdp,md.mdp拷贝到当前目录产生临时tpr文件

gmx grompp -f em.mdp -c complex_SOL.gro -p topol.top -o em.tpr -maxwarn 10

8.4 加离子，使得整个体系变为电中性这里选择分组时选择SOL对应的分组

gmx genion -s em.tpr -p topol.top -o system.gro -neutral -maxwarn 10

选择熔剂，即SOL对应的group数字，输入后回车。

8.5 能量极小化

gmx grompp -f em.mdp -c system.gro -p topol.top -o em.tpr -maxwarn 10 gmx mdrun -v -deffnm em

（有些时候会出现不平衡的警告，可以尝试使用em_conservative.mdp,该mpd文件步长更短，耗时更长，但更不容易崩溃）

9 进行NVT平衡（带位置约束）

创建NVT平衡输入文件：

# 生成tpr文件 gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr # 运行NVT平衡 gmx mdrun -v -deffnm nvt -nt 20

笔者的电脑运行大约需要15分钟

10 进行NPT平衡（带位置约束）

# 生成tpr文件 gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -p topol.top -o npt.tpr # 运行NPT平衡 gmx mdrun -v -deffnm npt -nt 20

笔者的电脑运行大约需要15分钟

11 产生索引文件

11.1产生索引文件

gmx make_ndx -f npt.gro

11.2 将蛋白、配体、离子定义为新组

输入蛋白、配体、离子前的数字，用“|”间隔开，即1|13|18回车。（由于配体、离子与蛋白的相互作用较为明显，所以这里任务之中将他们一起作为一个控温组和同一个消除屏东转动的组，）

11.3 命名为protein_lig

输入name 20 protein_lig，回车，（注意，这里的20是指19+1,请根据自己的情况填写）

11.4 将剩余的部分定义为新组

!20

11.5 命名为envir

输入name 21 envir，回车，（注意，这里的21是指19+1+1,请根据自己的情况填写）。

11.6 退出后记

输入q，回车。

12 进行分子动力学

12.1 生成md.tpr文件。（如果这里警告多，可以将maxwarn后的数值改大一些）

注：如果需要更改跑动力学的时间，需要修改md.mdp文件，默认为1ns，如果要跑100ns，在nsteps 后加两个0即可。

gmx grompp -f md.mdp -c pr.gro -p topol.top -o md.tpr -n index.ndx -maxwarn 10

12.2 运行

gmx mdrun -v -deffnm md

13 计算分子动力学结果参数

13.1分析RMSD两次都选蛋白部分（Protein），注意观测在模拟过程中骨架的RMSD波动是否较大,计算蛋白质的RMSD值。

gmx rms -f md.xtc -s md.tpr -o rmsd_protein.xvg -n index.ndx

13.2 分析RMSD

第一次选择蛋白部分（Protein）第二次选择配体（MOL），可以消除蛋白质整体的运动，观察小分子配体相对于蛋白质的运动。

gmx rms -s md.tpr -f md.xtc -o rmsd_complex.xvg -n index.ndx

13.3 计算RMSF选择蛋白部分（Protein）。

gmx rmsf -s md.tpr -f md.xtc -res -n index.ndx -o rmsf_protein.xvg

13.4 计算蛋白的回旋半径（Radius of Gyration）选择蛋白质的骨架（Backbone）。

gmx gyrate -s md.tpr -f md.xtc -o gyrate.xvg

13.5 计算溶剂可及表面积（SASA）

选择蛋白部分（Protein）。

gmx sasa -s md.tpr -f md.xtc -o sasa.xvg -tu ns

13.6 计算氢键数量依次选择蛋白部分（Protein）和配体部分（MOL）。

gmx hbond -s md.tpr -f md.xtc -n index.ndx -num hbond.xvg -tu ns

14 生成分子动力学模拟pdb文件

选择蛋白（protein）和非水分子（non-water）

pymol movie_noWater.pdb

后记

蛋白不同，小分子不同，会遇到很多bug，如果想完全学会用Gromacs进行分子动力学模拟，仍需翻阅论坛及各种学习网站，灵活处理警告与报错，学习之路，道阻且长啊！

参数可视化还有很多内容，最近事务繁忙，后面有机会在更新吧！

Reference

分子动力学模拟过程参考

（pr.mdp md.mdp em.mdp文件参考）

http://bbs.keinsci.com/thread-45791-1-1.html

nvt.mdp mpt.mdp文件参考

https://www.cnblogs.com/w-guangyu/p/7788651.html

https://www.cnblogs.com/w-guangyu/p/7788246.html

本文所用mdp文件通过百度网盘分享：

链接: https://pan.baidu.com/s/1fdl73CwmkbxNWqrOfJhYzw?pwd=yej4 提取码: yej4

配置参考

操作系统：Ubuntu 24.04.3 LTS

硬件信息：i5-14600kf, rtx5060ti 16G

蛋白设计交流QQ群：438723520（QQ群今后仅分享文献资料，模型文件请通过百度网盘获取）

作者：bioinforiver

2025年12月21日

本内容仅供学习交流，未经许可不许转载。

（ 本人系原作者，原文见：mp.weixin.qq.com/s/woJjm-sn4SitO302-TdwRQ）