（四）gromacs教程 - 帮助中心

1.介绍

肿瘤坏死因子-α （TNF-α）是强效促炎细胞因子，在类风湿关节炎病变、银屑病关节炎（PsA）和银屑病中起作用。TNF-α发挥体内生物学活性的结构是对称的同源三聚体；本文通过gromacs 模拟TNF-α的蛋白质结构；

2.环境及使用软件

算家云

- Ubuntu 24.04
- gromacs-gpu-2025.1
- SPDBV

3.模拟

#### 1. 获取并处理pdb文件

从PDB 数据库下载肿瘤坏死因子-α同源三聚体结构的pdb 文件1TNF.pdb ，手动或者使用wget命令。

``wget https://files.rcsb.org/download/1TNF.pdb`


在使用gromacs前需要对pdb 文件进行处理
- 氢原子
- C端氧原子
- 结晶水
- 缺失原子/残基/侧链
- 二硫键
- 带电残基
使用SPDBV 对蛋白质进行处理，需要运行在windows上，用SPDBV 打开文件，它可以替pdb 文件替换缺失的侧链，注意SPDBV 可能会在添加的侧链前加上奇怪的控制字符, 而且这些符号只能在文本编辑器中手工去掉!)，此外，有一些晶体结果会有配体，需要根据实际的需要进行删除处理。

如，在我们下载的1TNF.pdb 文件添加氧原子类型 OXT。

- File-Open PDB file-1omb.pdb - Build-Add C-terminal Oxygen (OXT) - File-Save-Current Layer-1tnf_fix.pdb - 然后删除1tnf_fix.pdb文件中添加到文件末尾的以SPDBV开头的行.


注：一定要使用SPDBV软件打开蛋白结构，然后在保存，它会自动修复一些结构问题，不然后继会报错。
#### 2. pdb2gmx 获取拓扑文件
pdb2gmx 命令可利用pdb 文件创建GROMACS 的输入坐标和拓扑文件. 拓扑文件包含了所有力场参数(基于所选择的力场)

`gmx_mpi pdb2gmx -ignh -ff amber99sb-ildn -f 1tnf.pdb -o 1tnf.gro -p 1tnf.top -water tip3p`


运行命令后会获得：

说明：如果不指定 -ff 和 -water 会出现列表选择力场和水模型
#### 3. 创建模拟盒子
说明：-bt dodecahedron 创建一个菱形十二面体盒子，计算效率最高； -d 1.2 设定分子到盒子边缘的最小距离单位，注意-d 不能小于0.9

`gmx_mpi editconf -f 1tnf.gro -o 1tnf-box.gro -bt dodecahedron -d 1.2`


#### 4. 真空中能量最小化
若只需要进行真空中的模拟，完成此步骤后就可以直接到成品模拟了。
编辑em-vac.mdp 文件，里面会指定每种计算类型的参数。

`; 传递给预处理器的一些定义 define = -DFLEXIBLE ; 使用柔性水模型而非刚性模型, 这样最陡下降法可进一步最小化能量

; 模拟类型, 结束控制, 输出控制参数 integrator = steep ; 指定使用最陡下降法进行能量最小化. 若设为cg则使用共轭梯度法 emtol = 500.0 ; 若力的最大值小于此值则认为能量最小化收敛(单位kJ mol^-1^ nm^-1^) emstep = 0.01 ; 初始步长(nm) nsteps = 1000 ; 在能量最小化中, 指定最大迭代次数 nstenergy = 1 ; 能量写出频率 energygrps = System ; 要写出的能量组

; 近邻列表, 相互作用计算参数 nstlist = 1 ; 更新近邻列表的频率. 1表示每步都更新 ns_type = grid ; 近邻列表确定方法(simple或grid) coulombtype = PME ; 计算长程静电的方法. PME为粒子网格Ewald方法, 还可以使用cut-off rlist = 1.0 ; 短程力近邻列表的截断值 rcoulomb = 1.0 ; 长程库仑力的截断值 vdwtype = cut-off ; 计算范德华作用的方法 rvdw = 1.0 ; 范德华距离截断值 constraints = none ; 设置模型中使用的约束 pbc = xyz ; 3维周期性边界条件`

`gmx_mpi grompp -f em-vac.mdp -c 1tnf-box.gro -p 1tnf.top -o em-vac.tpr`


运行后得到：输入文件 em-vac.tpr 和 参数文件 em-vac.mdp
然后进行能量最小化

`gmx_mpi mdrun -v -deffnm em-vac`

获得：日志文件 em-vac.log, 全精度轨迹文件 em-vac.trr, 能量文件 em-vac.edr, 能量最小化后的结构文件 em-vac.gro


!
#### 5. 向盒子中填充溶剂及离子，并进行能量最小化

`gmx_mpi solvate -cp em-vac.gro -cs spc216.gro -p 1tnf.top -o 1tnf-sol.gro`

-cp指定需要填充水分子的体系, 带模拟蛋白盒子, -cs指定使用SPC水模型进行填充, spc216是GROMACS统一的三位点水分子结构, -p修改体系的拓扑文件, 加入相应水分子的物理参数。


溶液体系的能量最小化文件em-sol.mdp

`define = -DFLEXIBLE


integrator      = steep
emtol           = 250.0
nsteps          = 5000
nstenergy       = 1
energygrps      = System

nstlist = 1 ns_type = grid coulombtype = PME rlist = 1.0 rcoulomb = 1.0 rvdw = 1.0 constraints = none pbc = xyz`


生成运行文件：

`gmx_mpi grompp -f em-sol.mdp -c 1tnf-sol.gro -p 1tnf.top -o em-sol.tpr`


会获得模拟用的输入文件em-sol.tpr ，向其中添加离子，来中和体系中的净电荷

`gmx_mpi genion -s em-sol.tpr -o 1tnf-em.gro -neutral -conc 0.15 -p 1tnf.top`

-neutral选项保证体系总的净电荷为零, 体系呈电中性, -conc选项设定需要的离子浓度(这里为0.15 M).genion默认使用的盐为NaCl. 如果你需要使用不同的阳离子和阴离子, 可使用 -pname(阳离子)和 -nname(阴离子)分别指定阴阳离子的名称(根据相应力场的ions.itp文件中离子的设定), 还可以使用 -pn和 -nn分别指定添加的阴离子数目. 运这个命令时, 会提示选择一个连续的溶剂分子组, 选择13 (SOL), 回车,


重新生成tpr 文件并执行能量最小化

`gmx_mpi grompp -f em-sol.mdp -c 1tnf-em.gro -p 1tnf.top -o 1tnf-em.tpr gmx_mpi mdrun -v -deffnm 1omb-em-sol`


!
#### 6. 位置限制预平衡模拟
对整个体系进行位置限制性模拟, 也就是对溶剂和离子进行弛豫同时保持蛋白质原子的位置不变. 在位置限制性模拟中会限制(或部分冻结)大分子中的原子位置, 而允许溶剂分子运动, 这样做像是将大分子浸入水中, 可以使体系进一步平衡.
水分子的弛豫时间约为10 ps, 大体系(大的蛋白质或脂)可能需要更长的平衡时间, 因此我们要进行超过10ps的位置限制性模拟(至少要长一个数量级). 我们将进行两个阶段的位置限制性模拟: 100 ps的NVT系综平衡和100 ps的NPT系综平衡. 模拟时我们使用的温度为300 K, 它接近于大多数实验条件的室温. 有些人会在310 K下进行模拟, 因为这个温度更接近于体温或生理温度.

NVT位置限制性模拟的参数文件 nvt-pr-md.mdp

`; 预处理选项 define = -DPOSRES ; 告诉GROMACS运行位置限制性模拟


; 运行控制参数
integrator               = md
dt                       = 0.002 ; 时间步长(单位为ps, 我们使用了2 fs). 只用于动力学积分器(如md), 能量最小化时不需要
nsteps                   = 50000 ; 模拟步数(总模拟时间为nsteps*dt)
; 输出控制参数
nstxout                  = 500   ; 输出模拟坐标的频率(nstxout=500且dt=0.002, 所以每1 ps输出一次)
nstvout                  = 500   ; 速度保存频率
nstenergy                = 500   ; 能量保存频率
nstlog                   = 500   ; log文件输出频率
energygrps               = system
; 近邻列表参数
nstlist                  = 5
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.0
; 静电和VDW参数
coulombtype              = PME  ; 长程静电相互作用的计算方法
pme_order                = 4    ; 三次插值
fourierspacing           = 0.16 ; FFT间隔
rcoulomb                 = 1.0  ; 计算静电作用的截断值(单位nm)
vdw-type                 = Cut-off
rvdw                     = 1.0
; 温度耦合部分非常重要, 必须正确填写.
tcoupl                   = v-rescale            ; 随机重新调整速度
tc-grps                  = Protein  Non-Protein ; 与控温器耦合的组(模型中的每个原子或残基都用一定的索引组表示), 对蛋白和非蛋白使用不同的组分开控制
tau_t                    = 0.1      0.1         ; 温度耦合的时间常数(单位ps). 必须每个tc_grps指定一个, 且顺序对应
ref_t                    = 300      300         ; 代表耦合的参考温度(即动力学模拟的温度, 单位K). 每个tc_grp对应一个ref_t
; 色散校正
DispCorr                 = EnerPres ; 校正VDW截断
; 不使用压力耦合
pcoupl                   = no                    ; NVT中不能使用压力耦合
; 初始速度选项
gen_vel                  = yes    ; 根据Maxwell分布随机产生速度
gen_temp                 = 300    ; 当你改变温度时, 别忘了改变gen_temp变量以生成速度
gen_seed                 = -1     ; 随机数生成器的种子

; 键约束选项 constraints = all-bonds ; 使用LINCS算法约束所有键 continuation = no ; 第一次运行 constraint_algorithm = lincs ; 约束算法 lincs_iter = 1 ; LINCS精度 lincs_order = 4 ; LINCS阶数, 与精度有关`


然后生成模拟输入文件，在运行模拟

`gmx_mpi grompp -f nvt-pr-md.mdp -c 1tnf-em.gro -r 1tnf-em.gro -p 1tnf.top -o 1tnf-nvt-pr.tpr gmx_mpi mdrun -deffnm 1tnf-nvt-pr`

NPT位置限制性模拟的参数文件 npt-pr-md.mdp

`define = -DPOSRES


integrator               = md
dt                       = 0.002
nsteps                   = 50000
nstxout                  = 500
nstvout                  = 500
nstfout                  = 500
nstenergy                = 500
nstlog                   = 500
energygrps               = System
nstlist                  = 5
ns-type                  = Grid
pbc                      = xyz
rlist                    = 1.0
coulombtype              = PME
pme_order                = 4
fourierspacing           = 0.16
rcoulomb                 = 1.0
vdw-type                 = Cut-off
rvdw                     = 1.0
Tcoupl                   = v-rescale
tc-grps                  = Protein  Non-Protein
tau_t                    = 0.1      0.1
ref_t                    = 300      300
DispCorr                 = EnerPres
; 压力耦合
Pcoupl                   = Parrinello-Rahman ; Parrinello-Rahman控压器.
Pcoupltype               = Isotropic         ; isotropic 指盒子可以平均地向各个方向(x, y,z)膨胀或压缩以维持一定的压力. 进行膜模拟时需要用semiisotropic.
tau_p                    = 2.0               ; 压力耦合的时间常数(单位ps).
compressibility          = 4.5e-5            ; 溶剂的压缩系数(4.5e-5为水在300 K和标准大气压下的压缩系数).
ref_p                    = 1.0               ; 压力耦合的参考压力(单位bar, 1大气压约为0.983 bar).
refcoord_scaling         = com
gen_vel                  = no                ; 不产生速度

constraints = all-bonds continuation = yes constraint_algorithm = lincs lincs_iter = 1 lincs_order = 4`


同样的操作

`gmx_mpi grompp -f npt-pr-md.mdp -c 1tnf-nvt-pr.gro -r 1tnf-nvt-pr.gro -p 1tnf.top -o 1tnf-npt-pr.tpr gmx_mpi mdrun -deffnm npt-nvt-pr-1omd`


!
说明：在mdrun，需要GPU加速的，可通过 -nb gpu 指定，默认是自动

- 1：用显卡做NB计算（-nb gpu -pme cpu -bonded cpu ） - 2：用显卡做NB+PME计算（-nb gpu -pcme gpu -bonded cpu ）多rank并行还需加上 -npme 1 - 3：用显卡做NB+BF计算（-nb gpu -pme cpu -bonded gpu ） - 4：用显卡做NB+PME+BF计算（-nb gpu -pme gpu -bonded gpu ）多rank并行还需加上 -npme 1 - 使用GPU计算能量分组（energygrps），在mdp 文件中，如果energygrps 分组了，gpu运算是不支持的

`bash energygrps = Protein Non-Protein # 计算蛋白质和非蛋白质的能量`


要取消能量分组，修改为：

`bash strong`

或直接删除 energygrps 行（GROMACS 会默认使用 System）。


#### 7. 成品模拟

NPT成品模拟使用的参数文件 npt-md-10ns.mdp

`integrator = md dt = 0.002 nsteps = 5000000 ; 10 ns


nstxout                 = 50000
nstvout                 = 50000
nstfout                 = 50000
nstenergy               = 50000
nstlog                  = 50000
energygrps              = System
nstlist                  = 5
ns-type                  = Grid
pbc                      = xyz
rlist                    = 1.0
coulombtype              = PME
pme_order                = 4
fourierspacing           = 0.16
rcoulomb                 = 1.0
vdw-type                 = Cut-off
rvdw                     = 1.0
Tcoupl                   = v-rescale
tc-grps                  = Protein  Non-Protein
tau_t                    = 0.1      0.1
ref_t                    = 300      300
DispCorr                 = EnerPres
Pcoupl                   = Parrinello-Rahman
Pcoupltype               = Isotropic
tau_p                    = 2.0
compressibility          = 4.5e-5
ref_p                    = 1.0
gen_vel                  = no

constraints = all-bonds continuation = yes constraint_algorithm = lincs lincs_iter = 1 lincs_order = 4`

`gmx_mpi grompp -f npt-md-10ns.mdp -c 1tnf-npt-pr.gro -r 1tnf-npt-pr.gro -p 1tnf.top -o npt-md-10ns.tpr gmx_mpi mdrun -deffnm npt-md-10ns


后台运行

nohup gmx_mpi mdrun -deffnm npt-md-10ns &
tail -f nohup.out


![[Pasted image 20250417090338.png]]
模拟完成后会获得npt-md-10ns的各种结果文件，以trr结尾的是模拟的轨迹文件，用于分析模拟结果。可以使用trjconv 命令压缩轨迹文件，使用-pbc nojump 让所有原子处于盒子中。

`gmx_mpi trjconv -f npt-md-10ns.trr -s npt-md-10ns.tpr -o npt-md-10ns.xtc -pbc nojump -ur compact -center``

会将轨迹文件压缩成xtc 文件；在提示中两次都选择：0 （system）

#### 8. 参数逐行解释

| 参数 | 值 | 说明 |
| ------------------------------ | ------------------- | --------------------------------- |
| integrator | md | 使用分子动力学算法 |
| dt | 0.002 ps | 积分步长（2 fs） |
| nsteps | 5,000,000 | 总步数（10 ns） |
| nstxout | 500 | 每500步（1 ps）输出一次坐标 |
| nstvout | 500 | 每500步输出一次速度（通常可禁用） |
| nstfout | 500 | 每500步输出一次受力（通常可禁用） |
| nstenergy | 500 | 每500步输出能量数据 |
| nstlog | 500 | 每500步更新日志文件 |
| energygrps | System | 仅计算体系总能量（GPU兼容） |
| nstlist | 5 | 每5步更新邻居列表（过低，建议20） |
| ns-type | Grid | 使用网格搜索法生成邻居列表 |
| pbc | xyz | 三维周期性边界条件 |
| rlist | 1.0 nm | 短程相互作用截断半径（建议≥1.2） |
| coulombtype | PME | 粒子网格Ewald法处理长程静电 |
| pme_order | 4 | PME插值阶数（默认值） |
| fourierspacing | 0.16 nm | PME网格间距（建议0.12） |
| rcoulomb | 1.0 nm | 静电截断半径（应与rlist一致） |
| vdw-type | Cut-off | 范德华力截断处理 |
| rvdw | 1.0 nm | 范德华力截断半径（建议≥1.2） |
| Tcoupl | v-rescale | 速度重缩放温度耦合 |
| tc-grps | Protein Non-Protein | 温度耦合分组（需确保组名存在） |
| tau_t | 0.1 0.1 ps | 温度弛豫时间 |
| ref_t | 300 300 K | 参考温度 |
| DispCorr | EnerPres | 对能量和压力进行长程校正 |
| Pcoupl | Parrinello-Rahman | 压力耦合算法 |
| Pcoupltype | Isotropic | 各向同性压力控制 |
| tau_p | 2.0 ps | 压力弛豫时间 |
| compressibility | 4.5e-5 bar⁻¹ | 体系压缩率（水溶液默认值） |
| ref_p | 1.0 bar | 参考压力 |
| gen_vel | no | 不生成初始速度（续跑时使用） |
| constraints | all-bonds | 约束所有键（建议改为h-bonds） |
| constraint_algorithm | lincs | 使用LINCS算法约束键长 |
| lincs_iter | 1 | LINCS迭代次数（建议2-4） |
| lincs_order | 4 | LINCS插值阶数 |