分析梳理--分子动力学模拟的常规步骤一（Gromacs）

作者，Evil Genius

现在大家羡慕的工作还是公务员、事业编，以前不理解，现在明白了。

群里说做科研不能有后顾之忧容易养闲人，实际上不可能，现行制度是低薪高奖励模式，例如博士进科研，月薪6000，这是保底的，有科研成果，比如发了一篇cancer cell，奖金20万，或者直接升级，没办法真的当闲人。

今天我们需要继续分子动力学，现在看文章基本上都有做一些药物预测，找一找靶点，形成一个闭环。

那么分子动力学，也是一个很大的工程，我们只能慢慢来了。

1.评估体系:明确三个问题:做什么、为什么做、怎么做

2.选择工具:

软件的选择:GROMACS

力场的选择:力场用来描述体系中最小单元间的相互作用，是对实验性质或量子化学计算结果拟合后生成的经验式。常见的有三类力场:全原子力场、联合力场、粗粒化力场。

3.初始结构:通过实验数据或者某些工具得到体系内的每一个分子的初始结构坐标文件之后，把这些分子按照一定的规则或是随机的放在一起，从而得到整个体系的初始结构，这也是模拟的输入文件。

4.输入参数:得到了结构输入文件，还需要力场参数输入文件，也就是针对体系的力场文件。这通常由所选用的力场决定，包括电荷，键合参数和非键参数等势能函数的输入参数

5.确定盒子:体系的大小通常由选用的盒子大小决定。必须对可行性与合理性做出评估，从而确定体系的大小，这依赖于具体的体系。

6.能量最小化:由于初始构象可能会存在两个原子靠得太近的情况(称之为badcontact),所以需要在正式模拟开始的第一步对体系进行能量最小化。

比较常用的能量最小化方法有两种:

最速下降法:快速移除体系内应力的好方法，但是接近能量极小点时收敛比较慢;

共轭梯度法:在能量极小点附近收敛效率高一些.

一般做能量最小化时都是先利用最速下降法进行优化，完成之后再对得到的构象利用共轭梯度法优化一次，这样做能有效地保证后续模拟的进行。

7.平衡模拟:需要设置适当的模拟参数，并且保证这些参数的设置与力场的构造过程相一致。

举个简单的例子，GROMOS力场是用范德华势双截断来定义范德华参数的，如果你用GROMOS力场的话也应该用双截断来处理范德华相互作用。

常见的模拟思路是:先在NVT下限制住溶质(剂)做限制性模拟，这是一个升温的过程，当温度达到你设定的温度后，接着做NPT模拟，此过程将调整体系的压强进而使体系密度收敛。

判断体系是否达到平衡的几种方式:

看能量(势能，动能和总能)是否收敛

看体系的压强，密度等等是否收敛

看体系的RMSD是否达到你能接受的范围

其他经验

8.成品模拟:经过一段时间的平衡模拟，在确定体系已经完全弛豫之后，就可以开始采集数据了。运行足够长时间的模拟以确定感兴趣的现象或是性质能够被观测到，并且务必确保此现象出现的可重复性。

9.数据分析:数据拿到手后，很容易通过一些可视化软件得到轨迹动画，但这并不能拿来发文章。真正的工作才刚刚开始--分析数据。你所感兴趣的现象或性质只是表面，隐含在它们之中的机理才是文章的主题。

流程图

GROMACS文件类型(功能分类)

参数文件

mdp运行参数文件,用作gmx grompp与gmx convert-tpr 的输入文件

mp2用作gmx xpm2ps 的输入文件

结构文件

gro GROMACS格式

g69 GROMOS-96格式

pdb 蛋白质数据库格式

通用结构格式:gro，g96，pdb，tpr，tpb或tpa

结构+质量(db):tpr，tpb，tpa，gro，g96或pdb ，用作分析工具的输入文件，当使用gro或pdb 时，从质量数据库中读取近似质量拓扑文件

top 体系拓扑(文本)

itp 包含拓扑(文本)

rtp 残基拓扑(文本)

ndx 素引文件

拓扑文件

拓扑文件包含体系中原子类型、电荷、成键和非键项及参数、各种分子数目等信息。

.top和.itp都叫拓扑文件(topology)

.top:作为体系的主拓扑文件，是grompp的输入文件之一。

.itp:included topology，被.top 通过#include 方式引用,储存力场参数和定义分子信息。itp文件不是必须的，也可以直接把itp文件里的内容直接写到top的合适位置里。拓扑文件里可以开头用分号进行注释，每个字段用[]标识。

top文件

itp文件

首先我们来看第一步：生成蛋白拓扑结构

看一看代码gmx pdb2gmx

主要是理解参数即可

Options to specify input files:

 -f      [<.gro/.g96/...>]  (protein.pdb)
           Structure file: gro g96 pdb brk ent esp tpr

Options to specify output files:

 -o      [<.gro/.g96/...>]  (conf.gro)
           Structure file: gro g96 pdb brk ent esp
 -p      [<.top>]           (topol.top)
           Topology file
 -i      [<.itp>]           (posre.itp)
           Include file for topology
 -n      [<.ndx>]           (index.ndx)      (Opt.)
           Index file
 -q      [<.gro/.g96/...>]  (clean.pdb)      (Opt.)
           Structure file: gro g96 pdb brk ent esp

Other options:

 -chainsep <enum>           (id_or_ter)
           Condition in PDB files when a new chain should be started (adding
           termini): id_or_ter, id_and_ter, ter, id, interactive
 -merge  <enum>             (no)
           Merge multiple chains into a single [moleculetype]: no, all,
           interactive
 -ff     <string>           (select)
           Force field, interactive by default. Use -h for information.
 -water  <enum>             (select)
           Water model to use: select, none, spc, spce, tip3p, tip4p, tip5p,
           tips3p
 -[no]inter                 (no)
           Set the next 8 options to interactive
 -[no]ss                    (no)
           Interactive SS bridge selection
 -[no]ter                   (no)
           Interactive termini selection, instead of charged (default)
 -[no]lys                   (no)
           Interactive lysine selection, instead of charged
 -[no]arg                   (no)
           Interactive arginine selection, instead of charged
 -[no]asp                   (no)
           Interactive aspartic acid selection, instead of charged
 -[no]glu                   (no)
           Interactive glutamic acid selection, instead of charged
 -[no]gln                   (no)
           Interactive glutamine selection, instead of charged
 -[no]his                   (no)
           Interactive histidine selection, instead of checking H-bonds
 -angle  <real>             (135)
           Minimum hydrogen-donor-acceptor angle for a H-bond (degrees)
 -dist   <real>             (0.3)
           Maximum donor-acceptor distance for a H-bond (nm)
 -[no]una                   (no)
           Select aromatic rings with united CH atoms on phenylalanine,
           tryptophane and tyrosine
 -[no]ignh                  (no)
           Ignore hydrogen atoms that are in the coordinate file
 -[no]missing               (no)
           Continue when atoms are missing and bonds cannot be made, dangerous
 -[no]v                     (no)
           Be slightly more verbose in messages
 -posrefc <real>            (1000)
           Force constant for position restraints
 -vsite  <enum>             (none)
           Convert atoms to virtual sites: none, hydrogens, aromatics
 -[no]heavyh                (no)
           Make hydrogen atoms heavy
 -[no]deuterate             (no)
           Change the mass of hydrogens to 2 amu
 -[no]chargegrp             (yes)
           Use charge groups in the .rtp file
 -[no]cmap                  (yes)
           Use cmap torsions (if enabled in the .rtp file)
 -[no]renum                 (no)
           Renumber the residues consecutively in the output
 -[no]rtpres                (no)
           Use .rtp entry names as residue names

我们来看看参数

-f 蛋白质的文件，这个一般我们采用的都是PDB文件，不过gro g96 pdb brk ent esp tpr这些格式的文件大家也要了解。

-o输出文件，一般我们输出gro格式

原子信息行的格式

其它参数

立场和溶剂选项

在分子对接和分子动力学模拟中，力场（Force Field）和水模型（Water Model）的选择需要根据具体的研究目标和计算资源来决定。不存在“绝对最好”的组合，只有“最合适”的选择。

目前，学术界的黄金标准是结合主流通用力场与三点水模型（如TIP3P），这一组合在计算精度和效率之间取得了最佳平衡。

主流力场推荐

CHARMM 和 AMBER 是目前应用最广、验证最充分的经典力场。此外，GAFF 和 OPLS 在药物设计中也有出色表现。

• CHARMM (Chemistry at HARvard Macromolecular Mechanics)
  • 特点与适用场景：在蛋白质和小分子的处理上都非常平衡，对配体结合亲和力的预测准确度高。其水模型（如TIPS3P）在蛋白质-配体体系中的相互作用能计算表现良好。
  • 常用版本：CHARMM36m、CGenFF (针对小分子)
• AMBER (Assisted Model Building with Energy Refinement)
  • 特点与适用场景：在蛋白质和核酸模拟领域有广泛的应用基础，对蛋白质结构的描述非常稳定。常与GAFF（通用琥珀力场）配合使用来处理药物小分子。
  • 常用版本：ff14SB (蛋白质)、ff99SB*ILDN、GAFF (小分子)
• 其他优秀选项
  • OPLS (Optimized Potentials for Liquid Simulations)：在工业界药物研发中非常流行，OPLS3e 版本在亲和力预测方面表现出色，甚至被一些研究认为比部分公开力场更精准。
  • OpenFF (Open Force Field Initiative)：新一代的、开源的、自动优化的力场。其 OpenFF Sage 版本性能已经可以与GAFF和CGenFF媲美，是一个很有潜力的新兴选择。

💧 水模型选择指南

水模型的选择对模拟结果的准确性至关重要。下面是几种主流模型的对比：

最广泛的水模型还是3点水模型(TIP3P和SPC)来源于上世纪80年代，四点水模型TIP4P是OPLS力场御用的(用TIP3P，SPC也完全可以)，OPLS2.0后御用SPC,然后还有五点水，六点水，七点水模型，比较小众，比如七点水唯一的亮点是可以在较大温度范围内较准确地重现密度和介电常数。

💡 如何根据你的目标做选择？

1. 我该选择哪个力场？

• 如果你刚开始接触分子模拟，或进行常规的蛋白质-配体对接：从 CHARMM36m 或 AMBER ff14SB 开始是最稳妥的选择。它们是社区内最主流、文档最丰富的力场，出现问题也更容易找到解决方案。
• 如果你的研究重点是药物筛选，需要精确计算多个类似物分子的结合自由能：可以关注 OPLS3e 或 OpenFF Sage，它们在药物分子的亲和力排序上表现优异。

1. 我该选择哪个水模型？

• 除非有特殊理由，首选 TIP3P 模型。它结合了不错的速度和足够的精度，被绝大多数主流力场（如CHARMM和AMBER）作为默认选项，并且有大量的研究数据可以参考。
• 当你的研究课题深入到分析结合口袋内水分子的能量贡献，并且计算资源充足时，可以尝试 TIP4P-Ew 以获得更准确的水结构信息

-ignh:忽略PDB文件中的氢原子，对NMR结构非常有用.否则，如果存在氢原子，它们的名称和顺序必须与GROMACS力场中的完全一样.对氢原子存在各种不同的约定，所以处理起来有时让人很头疼!如果你需要保留初始氢原子的坐标但需要重命名

-ter:交互式地为N末端和C末端分配电荷态

-inter:交互式地为Glu(谷氨酸),Asp(天冬氨酸),Lys(赖氮酸),Arg(精氨酸)和His指定电荷态，选择涉及二硫键的Cys(胱氨酸)

我们常用的命令就是

# 基本用法（交互式选择力场和水模型）
gmx pdb2gmx -f protein.pdb -o conf.gro -p topol.top

# 指定力场和水模型，忽略PDB中的氢
gmx pdb2gmx -f protein.pdb -o conf.gro -p topol.top -ff amber99sb -water tip3p -ignh

会输入三个文件

1.分子拓扑文件topol.top

2.位置限制文件posre.itp

3.后处理结构文件 pro_processed.gro

位置限制文件是什么？

位置限制文件（通常名为 posre.itp）是分子动力学模拟中一个非常重要的辅助文件。它的核心作用是在模拟过程中，将体系中的某些原子（通常是蛋白质主链的重原子）“栓系”在它们的初始位置附近，允许小幅度振动，但限制其发生大范围的移动。

🎯 核心作用：为什么要加“位置限制”？

在分子动力学模拟的标准流程中，位置限制文件主要在两个阶段发挥关键作用：

体系平衡阶段（最常用）：

问题：当你将一个蛋白质放入充满溶剂（水）的盒子中，并添加离子后，溶剂分子和离子在初始状态下可能处于能量上不合理的、非常拥挤的位置。

后果：如果放开所有原子进行模拟，溶剂分子会以极大的力量“撞击”蛋白质，可能导致蛋白质骨架扭曲、侧链被撞飞，甚至整个结构崩溃。

解决方案：在平衡模拟（NVT和NPT） 开始时，对蛋白质的重原子施加位置限制。这样，溶剂分子和离子可以自由驰骋、找到自己的平衡位置，而蛋白质则能稳定地保持在原始构象附近。待溶剂环境稳定后，再在后续的生产模拟（Production Run） 中放开限制。

部分结构固定模拟：

你可能只想研究某个活性位点环区的柔性，而将蛋白质的其他部分视为刚性背景。这时，可以只为需要固定的部分生成位置限制文件，其他部分则自由运动。

⚙️ 工作原理：它是如何起作用的？

位置限制通过向体系的势能方程中添加一个额外的谐振势能项（Harmonic Restraint） 来实现。

数学表达：V_restraint = 1/2 * k * (r - r_ref)^2

V_restraint：施加的限制势能。

k：力常数（Force Constant），其大小决定了限制的强度。pdb2gmx 默认生成的值是 1000 kJ/(mol·nm²)。

r：原子在模拟过程中的实时位置。

r_ref：原子的参考位置，通常就是该原子在输入结构文件（如 conf.gro）中的初始坐标。

产生的效果：原子偏离其初始位置越远，势能就越高，受到的回拉力就越大。这就好比用一根根橡皮筋将每个原子固定在了它们的起点。

生活很好，有你更好。