网络药理学：14、分子动力学模拟MD：水中的溶菌酶（蛋白前置处理、pdb2gmx得到蛋白拓扑文件、添加盒子、溶剂化、离子化、能量最小化、NVT温度平衡、NPT压强平衡、矫正轨迹）

原指引中文机翻见：https://blog.csdn.net/weixin_42486623/article/details/128756587
或者直接原指引网页自己机翻：http://www.mdtutorials.com/gmx/lysozyme/01_pdb2gmx.html

1.蛋白PDB前置处理

1.1去水

将下载的PDB格式文件进行去水操作：

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

也可以用pymol等软件去水。

1.2检查是否缺失氨基酸残基或原子

检查下载的PDB格式文件是否含有MISSING，无输出说明没有缺失：

cat 1aki.pdb | grep -A 10 "MISSING"  # 输出含有MISSING的行及其后面十行

也可以记事本打开后ctrl+F搜索。

以2wtt.pdb为例，查看缺失时是什么情况

其中MISSING RESIDUES意味着缺失氨基酸残基。
RES列是氨基酸残基名称，C列是所在链，SSSEQI是该氨基酸残基在PDB数据库的序列号。

MISSING ATTOM意味着缺失原子。

2.运行 pdb2gmx得到蛋白拓扑文件

2.1.pdb2gmx参数

在终端进入包含你的PDB文件的目录，然后使用以下命令运行 pdb2gmx：

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

参数说明：

• -f 【your_structure.pdb】：输入的PDB文件。
• -o 【processed.gro】：输出的拓扑文件，通常为 .gro 格式。
• -water 【spce】：指定使用的水模型，常用SPC/SPCE水模型。

其他常用参数如下：

• -p 【topol.top】：指定生成的拓扑文件名。
• -ignh：忽略 PDB 文件中的 H 原子，而是以选择的力场所期望的方式加上全氢（即所有氢原子）。请注意，该参数很有可能会忽略掉一些原蛋白中重要且特殊的氢原子，而加上普通的氢原子，导致结果有偏差。譬如导致蛋白配体结合时氢键有误。

2.2.选择力场

运行 pdb2gmx 后，系统会提示你选择一个力场。你会看到一个力场列表，选择适合你的模拟的力场（可以见上文的2.5）。
力场列表如下：

 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)

在这里我们根据教程选择15，输入15后回车。

2.3.gro、top、itp文件解读

根据结果显示，我们可以看到该蛋白质的净电荷为 +8e，且新生成文件：1AKI_processed.gro、topol.top、prosre.itp。

2.3.1.gro文件解读

1AKI_processed.gro：包含力场中定义的所有原子，用pymol打开和1AKI_clean.pdb是一模一样的，只不过 换成了GROMACS 格式的结构文件。

官方补充：许多用户认为 .gro 文件是必需的。这不是真的。GROMACS 可以处理许多不同的文件格式，其中 .gro 只是写入坐标文件的命令的默认值。这是一种非常紧凑的格式，但精度有限。例如，如果您更喜欢使用 PDB 格式，您需要做的就是指定一个带有 .pdb 扩展名的适当文件名作为输出。

2.3.2.itp文件解读

posre.itp：“position restraints”，即位置约束。在模拟开始时，使用位置约束可以帮助系统更快地达到平衡，尤其是在模拟开始时，防止结构发生不希望的大规模移动。

2.3.3.top文件解读

topol.top用Include关键字引用其他拓扑文件，我们输入以下命令进行查看

cat topol.top | grep -A 15 "Include"

输出如下：

可以看到一共有四个部分，分别是

• forcefield parameters：力场相关
• Posistion restraint file：位置约束相关
• water topology：水拓扑相关
• topology for ions：离子拓扑相关

力场相关：

• forcefiled parameters：的下一行即用的力场文件，因为我们选择了OPLS-AA力场，所以这里我们可以看到#include "oplsaa.ff/forcefield.itp"
• moleculetype：两个参数，其中"Name" 定义了分子名称
• atoms：定义了蛋白质中的信息：
  • residue：氨基酸残基名称，如LYS
  • atom： 原子名称，如N
  • cgnr：电荷分组的编号，如这里前6行都是1，说明都是一组的。
  • charge：电荷。
  • mass：原子质量。

位置约束相关：
可以看到引用了我们pdb2gmx命令生成的另一个文件posre.itp

水拓扑相关：
因为我们在pdb2gmx命令中指定了水模型，所以可以看到spce.itp

离子拓扑相关：
引用了OPLS-AA力场中自带的ions.itp

---

同时top文件用方括号[]来定义不同的节，每个节包含了特定类型的信息。

• [ defaults ]：这个节定义了默认的力场参数，如原子类型之间的非键相互作用的默认参数。
• [ atomtypes ]：这个节列出了在模拟中使用的所有原子类型及其相关参数，如原子质量、电荷、Lennard-Jones参数等。
• [ moleculetype ]：这个节定义了一个分子类型，包括它的名字和它所包含的原子。每个分子类型下面通常会包含一个或多个 [ atoms ]、[ bonds ]、[ pairs ]、[ angles ]、[ dihedrals ] 等节，如下：
  • [ atoms ]：在这个节中，列出了分子中的所有原子，并为每个原子指定了其类型、电荷、原子质量以及它在分子中的位置。
  • [ bonds ]：这个节列出了分子中的所有键，包括形成键的原子索引和键的类型。
  • [ pairs ]：这个节用于指定需要特殊处理的非键相互作用，通常用于1-4相互作用。
  • [ angles ]：这个节列出了分子中的所有键角，包括形成键角的原子索引和键角的类型。
  • [ dihedrals ]：这个节列出了分子中的所有二面角，包括形成二面角的原子索引和二面角的类型。
• [ system ]：这个节定义了整个模拟系统的名称，并列出了组成系统的所有分子类型及其数量。
• [ molecules ]：这个节列出了系统中所有分子类型的数量，用于告诉GROMACS在模拟中应该包含多少个每种类型的分子。

---

最后我们输入命令查看整个系统的净电荷：

cat topol.top | grep "qtot" | tail -n 1

输出：

1960   opls_272    129    LEU     O2    682       -0.8    15.9994   ; qtot 8

这里qtot 8，代表着该蛋白质的净电荷为 +8e。

3.运行editconf添加盒子

运行命令：

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

参数解读：

• -f：输入文件名
• -o：输出文件名
• -c： 使盒子居中
• -d 【1.0】：使盒子和蛋白边缘距离1.0 nm
• -bt 【cubic】：定义一个立方体盒子。

运行完后看到gromacs reminds you + 一句名人名言，且新增一个文件IAKI_newbox.gro则是成功。该文件用pymol打开和IAKI_clean.pdb打开是一样的。

4.运行solvate实现溶剂化（即往盒子添加水分子）

4.1.运行solvate命令

运行命令：

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

参数解读：

• -o：输出文件名
• -cp 【1AKI_newbox.gro】：指定了包含溶质（通常是蛋白质或其它生物大分子）坐标的.gro文件。
• -cs 【spc216.gro】：指定了包含溶剂分子（通常是水）坐标的.gro文件。在这里的spc216.gro是指定的水模型自带的文件。
• -p 【topol.top】：指定生成的拓扑文件名。

4.2.结果pymol展示

运行完后看到gromacs reminds you + 一句名人名言，且新增一个文件IAKI_solv.gro则是成功。该文件用pymol打开如下：

很明显看到是个正方体盒子。

4.3.拓扑查看水分子

另外生成的'#topol.top.1#'是因为gromacs如果新生成的文件和旧文件重名，会将旧文件用#'包裹起来。

我们打开新的topol.top文件，看看里面新增了什么内容。

cat topol.top | grep -A 10 "molecules"

比起原来的拓扑文件，新增了最后一行SOL 10644，意味着增加了10644个水分子。

5.离子化

5.1.生成离子化所需的mdp文件（ions.mdp）

我们可以直接点击教程的“在此处”或“here”。

新打开页面的内容如下：

; 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

然后我们输入命令，新建文件并写入相关内容：

echo "省略内容" > ions.mdp

5.2.运行grompp添加离子

输入命令：

gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral

然后显示15个选项让我们选择，因为我们在3.4的步骤最后看到的是SOL 10644，那么在这里我们选择第13组SOL，输入13后回车，生成新文件1AKI_solv_ions.gro。
用pymol打开后如下，可以看到8个小球，实际上是8个离子：

我们再次查看拓扑文件：

cat topol.top | grep -A 10 "molecules"

可以看到是8个氯离子。

6.能量最小化

如果是用其他软件，譬如chem3D做能量最小化的话，可以看：https://www.bilibili.com/video/BV1om421G7nD?p=2

6.1.生成能量最小化所需的mdp文件（minim.mdp）

6.1.1.得到mdp文件内容

我们可以直接点击教程的“使用以下输入”或“here”。

新打开页面的内容如下：

; 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

然后我们输入命令，新建文件并写入相关内容：

echo "省略内容" > minim.mdp

6.2.2.mdp文件内容详解与修改

内容详解：

• integrator = steep: 这个参数指定使用最速下降法（steepest descent）进行能量最小化。
• emtol = 1000.0: 当系统中的最大力小于1000.0 kJ/mol/nm时，最小化过程将停止。
• emstep = 0.01: 这个参数设置了最小化过程中的步长。
• nsteps = 50000: 这是最小化过程中允许的最大步数。

其中emtol越大，nsteps越小的话，整个体系可能是越不稳定的。所以我们可以修改一下这两个参数实现更精准的模拟结果，将emtol改为900，将nsteps改为6000。

sed -i 's/emtol       = 1000.0/emtol       = 900.0/' minim.mdp
sed -i 's/nsteps      = 50000/nsteps      = 60000/' minim.mdp

然后输入命令：

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

6.2.能量最小化

gmx mdrun -v -deffnm em

参数详解：

• -v：显示进度，因为mdrun步骤需要消耗算力，所以过程可能较长，如果不加该参数，只会报错显示或成功时显示名人名言。
• -nb gpu：使用gpu加速

结果显示：

可以看到在第830步就达到最大力，其中整个体系的潜在能量Potential Energy = -5.8681488e+05，最大力Maximum force = 8.6097821e+02 on atom 736约为860，施加在第736个原子上。

同时我们得到文件如下：

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

6.3.得到能量最小化的xvg文件

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

我们看到选择列表如下：

输入10后回车，表示我们选择潜在能量。
再输入0，回车，表示终止输入。

打开qtgrace软件，拖拽新生成的potential.xvg如下：

注意：如果你和我一样，主机是win11，用wsl2生成的一个ubuntu实例，要将ubuntu实例中的文件拷贝到主机文件夹中，再拖拽到qtgrace打开，否则是打不开的，显示初始空白界面，而且不会报错，并且会影响到后续文件的打开！

横轴坐标就是步数，纵轴坐标为潜在能量。如果整体曲线趋势最后是趋于平缓的，那么证明能量最小化成功，否则失败，需要调整第7步的参数。

7.NVT温度平衡

7.1.生成NVT温度平衡所需的mdp文件（nvt.mdp）

7.1.1.得到mdp文件内容

我们可以直接点击教程的“在此处”或“here”。

新打开页面的内容如下：

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

然后我们输入命令，新建文件并写入相关内容：

echo "省略内容" > nvt.mdp

7.1.2.mdp文件内容详解

• nsteps= 50000：意味着模拟将运行50000步。由于每一步的时间步长（dt）是2飞秒，所以总的模拟时间是50000×0.002 fs=100 ps（皮秒）。

• dt=0.002：单位是飞秒（fs），即每一步模拟的时间间隔是2飞秒。

• nstxout = 500: 每500步保存一次坐标数据，这相当于每1.0皮秒（ps）保存一次。

• nstvout = 500: 每500步保存一次速度数据，也相当于每1.0皮秒（ps）保存一次。

• nstenergy = 500: 每500步保存一次能量数据，同样相当于每1.0皮秒（ps）保存一次。

• nstlog = 500: 每500步更新一次日志文件，这意味着每1.0皮秒（ps）更新一次日志。

7.2.NVT温度平衡

输入命令：

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

这里生成的mdrun命令是可以中止的，如果你途中有事的话，可以ctrl c中止命令，结果会临时存储在一个nvt.cpt文件上，下次运行命令时添加参数-cpi nvt.cpt，就可以在上次运行的基础上续跑了。

同时我们得到文件如下：

• nvt.log：NVT温度平衡进程的 ASCII 文本日志文件
• nvt.edr：二进制能量文件
• nvt.trr：二进制全精度轨迹
• nvt.gro：NVT温度平衡结构

7.3.得到NVT温度平衡的xvg文件

gmx energy -f nvt.edr -o temperature.xvg

我们看到选择列表如下：

输入16后回车，表示我们选择温度Temperature。
再输入0，回车，表示终止输入。

打开qtgrace软件，拖拽新生成的temperature.xvg如下：

待更：这个图怪怪的，可能是因为我还没有处理。不过可以通过纵坐标看出来，它的波动幅度在296到306之间，实际上是很小的范围，所以也算是平衡了吧hhh。

8.NPT压强平衡

8.1.得到NPT压强平衡所需的mdp文件（npt.mdp）

我们可以直接点击教程的“在此处”或“here”。

新打开页面的内容如下：

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 

然后我们输入命令，新建文件并写入相关内容：

echo "省略内容" > npt.mdp

8.2.NPT压强平衡

输入命令：

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

8.3.得到NPT压强平衡的xvg文件

输入命令：

gmx energy -f npt.edr -o pressure.xvg

我们看到选择列表如下：

输入18后回车，表示我们选择压强Pressure。
再输入0，回车，表示终止输入。
会得到一个pressure.xvg，但是研究意义不大，所以不展开讨论了。

9.正式MD模拟

9.1.得到MD所需的mdp文件（md.mdp）

9.1.1.得到mdp文件内容

我们可以直接点击教程的“在此处”或“here”。

新打开页面的内容如下：

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 

然后我们输入命令，新建文件并写入相关内容：

echo "省略内容" > md.mdp

9.1.2.mdp文件内容详解

对于正式的文献研究，这里的nsteps需要到100000ps（皮秒），即100ns（纳秒）。

9.2.正式MD模拟

输入命令：

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 -v # cpu计算，大概需要4-5分钟。
gmx mdrun -deffnm md_0_1 -nb gpu # gpu计算

10.矫正轨迹

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

参数详解：

• -s 【md_0_1.tpr】: 指定拓扑文件（.tpr），这是包含了模拟系统的所有参数的文件，如原子坐标、力场参数、盒子大小等。于第九步正式MD模拟中生成。
• -f 【md_0_1.xtc】: 指定输入轨迹文件（.xtc），这是包含了模拟过程中所有时间点的原子坐标的压缩轨迹文件。于第九步正式MD模拟中生成。
• -o 【md_0_1_noPBC.xtc】: 指定输出轨迹文件的名称（.xtc），这个文件将包含经过处理后的轨迹。
• -pbc 【mol】: mol 表示将以分子为单位进行周期性边界条件的处理，即确保分子在输出轨迹中是完整的，而不是被盒子边界切断。
• -center: 这个参数将系统中心置于模拟盒子的中心。

我们看到选择列表如下：

选择 1 （“Protein”）后回车，表示我们选择蛋白质作为要居中的组
再选择 0 （“System”）后回车，表示我们将整个体系作为输出。