GROMACS ：Lysozyme in Water

1. 准备工作

首先，确保你已经安装了 GROMACS 并且能够正常运行它。本教程的目标是模拟溶于水的溶菌酶（Lysozyme）。我们将通过以下步骤完成从下载蛋白质结构到进行分子动力学模拟的全过程。

步骤 1: 获取蛋白质结构

你可以从 RCSB PDB 数据库 下载溶菌酶（Lysozyme）的三维结构文件，打开网站搜索1AKI编号，点开1AKI，可以选在下载的结构，选在pdb文件进行下载。假设你已下载了 lysozyme.pdb 文件。

步骤 2: 处理获取蛋白质结构

下载结构后，您可以使用 VMD、Chimera、PyMOL 等查看程序可视化结使用 vi、emacs （Linux/Mac） 或记事本 （Windows） 等纯文本编辑器处理删除。

grep是linux的处理命令，命令是将含有HOH的行删除并保存到1AKI_clean.pdb

grep -v HOH 1aki.pdb > 1AKI_clean.pdb

对于win中建议使用notepad++，直接删除处理

注意：这部分要进行处理不然后续会出现top和坐标文件不对应，报错。同时这部分在gromacs中的残疾结构不存在导出top文件会报错

步骤 3: 获取并选择力场

准备好输入到第一个 GROMACS 模块 pdb2gmx 中。pdb2gmx 的目的是生成三个文件：

1. 分子的拓扑结构（力场文件）。
2. 位置约束文件（主要是用来固定蛋白质的内部结构，防止变性，失去特性）。
3. 后处理的结构文件（gro和pdb坐标文件）。

:grey_exclamation:注：对于蛋白质这类物质可以使用这个命令生成，这是因为蛋白质是由固定的残基构成，残基种类是确定的，这个命令只是索引坐标中的残基，将相应的残基整理输入到一个top文件中。

在 GROMACS 中，为了生成模拟所需的拓扑文件，你需要选择适合蛋白质的力场。我们在这里使用 GROMOS 力场。

运行以下命令：

gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water spce

-f 1AKI_clean.pdb 输入那个文件
-o 1AKI_processed.gro 输出那个文件
-water spce 加入水分模型的选择，会同时生成水分子力场文件

输出如下，选择使用哪个力场，键入 15，然后按 Enter：

Select the Force Field:
From '/usr/local/gromacs/share/gromacs/top':
 1: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
 2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
 3: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
 4: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
 5: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
 6: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
 7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
 8: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
 9: GROMOS96 43a1 force field
10: GROMOS96 43a2 force field (improved alkane dihedrals)
11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
14: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9)
15: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)

该命令将生成蛋白质的处理后的结构文件 1AKI_processed.gro、topol.top 和 posre.itp。1AKI_processed.gro 是一个 GROMACS 格式的结构文件，其中包含力场中定义的所有原子。topol.top 文件是系统拓扑（稍后会详细介绍）。posre.itp 文件包含用于约束重原子位置的信息（稍后将详细介绍）。

2. 水和溶剂化

1. 使用 editconf 模块定义盒子尺寸。
2. 使用 solvate 模块（以前称为 genbox）将盒子装满水。

步骤 3: 创建模拟盒子

接下来，使用 GROMACS 创建一个容纳蛋白质的模拟盒子，并填充适当的水分子。使用 editconf 命令创建一个立方体的模拟盒子，并设置蛋白质与盒子边缘的距离为 1.0 nm：

推荐菱形十二面体，因为它的体积是相同周期距离的立方框的 ~71%，从而节省了需要添加以溶解蛋白质的水分子数量。

让我们用 editconf 定义这个盒子：

gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic

上述命令将蛋白质在框 （-c） 中居中，并将其放置在距框边缘至少 1.0 nm （-d 1.0） 的位置。框类型定义为多维数据集 （-bt cubic）。

注意：指定 1.0 nm 的溶质盒距离意味着蛋白质的任意两个周期性图像之间至少有 2.0 nm。这个距离对于模拟中常用的几乎任何截止方案都足够了，这个地方涉及截断距离，力随着距离增加而减小，当距离达到一定程度就可以忽略，节省计算资源。阶段距离要大于盒子边长的一半。这是因为盒子是周期性的，防止同一原子和自身相互作用

这将生成一个新的结构文件 1AKI_newbox.gro，其中包含了蛋白质和一个立方体盒子。

步骤 4: 添加水分子

使用 solvate 命令将水分子填充到模拟盒子中。我们使用 spc216.gro 水模型，但你也可以选择其他水模型（如 TIP3P 或 SPC）：

gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top

（-cp） 输入上步骤生成的带有盒子的文件，溶剂的配置 （-cs） 是标准 GROMACS 安装的一部分。我们使用的是 spc216.gro，这是一个通用的平衡 3 点溶剂模型。还有一些其他模型 spc216.gro 作为 SPC、SPC/E 或 TIP3P 水的溶剂配置。-o 1AKI_solv.gro -p topol.top，这两个都是输出输出 1AKI_solv.gro和topol.top。

该命令会在你的模拟盒子中添加水分子，并更新拓扑文件 topol.top。

3. 添加离子

pdb2gmx 的输出告诉我们，该蛋白质的净电荷为 +8e（基于其氨基酸组成，Total charge 8.000 e）。为了中和系统的电荷并使其接近生理条件，我们需要向系统中添加适量的离子（如钠离子和氯离子）。

步骤 5: 中和系统电荷

首先，使用 gmx grompp 命令生成预处理文件(ions.md问价如下)：

; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching 
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = cutoff    ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr

mdp文件是用来设置规则条件的，比如什么压力、温度、步长等等
-f ions.mdp -c 1AKI_solv.gro -p topol.top这三个都是输入命令
-o 输出的时二进制文件，这个文件提交给电脑运算进行操作，跑模拟的时候使用

然后，使用 genion 命令添加离子并中和系统的电荷。此命令会自动选择添加钠离子（Na+）和氯离子（Cl-）来确保系统电中性：

gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral

屏幕上会输出选择页面，选择第 13 组 “SOL” 进行嵌入离子，意思是添加离子替换SOL分子，SOL在这里代表水分子。
请注意，-pname NA 和 -nname CL 参数表示添加的离子类型，阳离子和阴离子，-neutral 参数确保最终系统的电荷为零，最终会添加8g个氯离子来中和体系。

4. 能量最小化

pdb2gmx 生成的那个 posre.itp 文件吗？我们现在就使用它。posre.itp 的目的是对蛋白质的重原子（任何不是氢的东西）施加位置应变力。允许移动，但前提是克服了大量的能量损失。位置约束的实用性在于，它们允许我们平衡蛋白质周围的溶剂，而不会增加蛋白质结构变化的变量，蛋白质最重要的就是分子的结构如果变形就会变性，没有特有性质。这个文件被包含在topol.top中。

在进行分子动力学模拟之前，必须对系统进行能量最小化，以消除不合理的原子间相互作用。这一步会降低最终md的不合理性，降低计算时间。甚至是报错

步骤 6: 能量最小化

首先，准备一个能量最小化的参数文件（例如 minim.mdp）。然后，运行以下命令进行能量最小化：

文件内容如下：

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

gmx grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

确保坐标文件和top文件之间相互匹配
-v 将运行过程打印出来：它使 mdrun 变得冗长，以至于它在每一步都将其进度打印到屏幕上。
-deffnm 标志将定义输入em.tpr和输出的文件名em.*。

• em.log：EM 进程的 ASCII 文本日志文件
• em.edr：二进制能量文件
• em.trr：二进制全精度轨迹
• em.gro：能量最小化结构

运行完成后，将生成能量最小化后的结构文件em.gro。

em.edr 文件包含 GROMACS 在 EM 期间收集的所有能量项。您可以使用 GROMACS 能量模块分析任何 .edr 文件，为了查看模拟过程参数是否合适：

gmx energy -f em.edr -o potential.xvg

键入 “10 0” 以选择 Potential （10）;零 （0） 终止输入。

5. 生成动力学模拟

步骤 7: 温度和压力耦合

在开始分子动力学模拟之前，通常需要先进行温度和压力的平衡。我们将首先进行 NVT（恒容恒温）平衡，然后进行 NPT（恒压恒温）平衡。这两部也是预平衡体系。

1) NVT 平衡（恒容恒温）

准备 NVT 平衡的参数文件（如 nvt.mdp），然后运行以下命令：
nvt.mdp文件内容如下：

title                   = OPLS Lysozyme NVT equilibration 
define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 50000     ; 2 * 50000 = 100 ps
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Nonbonded settings 
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
tau_t                   = 0.1     0.1           ; time constant, in ps
ref_t                   = 300     300           ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl                  = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 300       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr

gmx mdrun -deffnm nvt

2) NPT 平衡（恒压恒温）

接下来，进行 NPT 平衡，使用适当的参数文件（如 npt.mdp）：

npt.mdp文件如下：

title                   = OPLS Lysozyme NPT equilibration 
define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 50000     ; 2 * 50000 = 100 ps
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = yes       ; Restarting after NVT 
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Nonbonded settings 
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
tau_t                   = 0.1     0.1           ; time constant, in ps
ref_t                   = 300     300           ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
pcoupltype              = isotropic             ; uniform scaling of box vectors
tau_p                   = 2.0                   ; time constant, in ps
ref_p                   = 1.0                   ; reference pressure, in bar
compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
refcoord_scaling        = com
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = no        ; Velocity generation is off 

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr

gmx mdrun -deffnm npt

6. 生产模拟

步骤 8: 进行生产模拟

一旦平衡完成，就可以进行生产阶段的分子动力学模拟。准备生产模拟的参数文件（例如 md.mdp），然后运行：

md.mdp文件如下：

title                   = OPLS Lysozyme NPT equilibration 
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 500000    ; 测试建议调小2 * 500000 = 1000 ps (1 ns)
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 0         ; suppress bulky .trr file by specifying 
nstvout                 = 0         ; 0 for output frequency of nstxout,
nstfout                 = 0         ; nstvout, and nstfout
nstenergy               = 5000      ; save energies every 10.0 ps
nstlog                  = 5000      ; update log file every 10.0 ps
nstxout-compressed      = 5000      ; save compressed coordinates every 10.0 ps
compressed-x-grps       = System    ; save the whole system
; Bond parameters
continuation            = yes       ; Restarting after NPT 
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighborsearching
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = Protein Non-Protein   ; two coupling groups - more accurate
tau_t                   = 0.1     0.1           ; time constant, in ps
ref_t                   = 300     300           ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPT
pcoupltype              = isotropic             ; uniform scaling of box vectors
tau_p                   = 2.0                   ; time constant, in ps
ref_p                   = 1.0                   ; reference pressure, in bar
compressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Dispersion correction
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel                 = no        ; Velocity generation is off 

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
gmx mdrun -deffnm md_0_1
# 使用gpu加速
gmx mdrun -deffnm md_0_1 -nb gpu

生产模拟将会输出原子轨迹文件（.xtc）和其他模拟结果。

7. 数据分析

1 轨迹处理

第一个是 trjconv，它用作后处理工具，用于去除坐标、校正周期性或手动更改轨迹（时间单位、帧频率等）。在本练习中，蛋白质可能看起来“破碎”或可能“跳”到盒子的另一侧，因为盒子是周期性的，就是好多盒子堆积起来，虽然只显示一个盒子。要考虑此类行为，请发出以下内容：

gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center

选择 1 （“Protein”） 作为要居中的组，选择 0 （“System”） 作为输出。我们将根据这个 “修正” 的轨迹进行所有分析。让我们先看看结构稳定性。

2 rms分析

gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns

选择 4（“Backbone”）作为最小二乘拟合和 RMSD 计算的组。
-s md_0_1.tpr：指定了输入的 TPR 文件（模拟的启动文件）。这里的文件名是 md_0_1.tpr，它包含了模拟的拓扑和其他参数。
-f md_0_1_noPBC.xtc：指定了输入的轨迹文件（这里是 md_0_1_noPBC.xtc）。该文件包含了模拟过程中分子的位置坐标（每个时间点的结构）。
-o rmsd.xvg：输出文件的名称。rmsd.xvg 是将会包含计算出的 RMSD 值的文件。
-tu ns：设置输出的时间单位为纳秒（ns），即使输入的轨迹是以皮秒（ps）为单位进行存储的。-tu 标志将以 ns 为单位输出结果，即使轨迹以 ps 为单位编写。这样做是为了使输出清晰（特别是如果你有一个长时间的模拟 - 1e+05 ps 看起来不如 100 ns）

蛋白质的回转半径是其紧凑性的量度。如果蛋白质稳定折叠，它可能会保持相对稳定的 Rg 值。如果蛋白质去折叠，它的 Rg 会随着时间的推移而变化。让我们在模拟中分析溶菌酶的回转半径：

蛋白质的回转半径

gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg

xvg文件查看可以使用qtgrace查看，上述运行文件在这个链接，文件较大删除轨迹文件