分子动力学仿真软件：GROMACS_（4）.构建初始拓扑模型

构建初始拓扑模型

在进行分子动力学模拟之前，构建一个准确的初始拓扑模型是至关重要的步骤。拓扑模型描述了分子系统的原子组成、化学键、角度、二面角等信息，以及这些信息如何影响分子的运动和相互作用。GROMACS 提供了多种工具和文件格式来帮助用户构建和验证这些拓扑模型。

1. 创建分子拓扑文件

分子拓扑文件在 GROMACS 中通常以 .top 文件的形式存在，包含了分子的所有必要信息。这些信息包括但不限于原子类型、化学键、角度、二面角、非键相互作用等。创建拓扑文件通常需要以下几个步骤：

1.1. 定义原子类型

原子类型是分子拓扑文件的基础，它描述了原子的属性，如质量、电荷、Lennard-Jones 参数等。这些参数通常在力场文件中定义，例如 gromos54a7.ff 力场中的 atomtypes.atp 文件。

例子：定义简单水分子的原子类型

假设我们使用 GROMOS54A7 力场来定义一个简单的水分子（H2O）。我们可以从力场文件中提取或定义如下原子类型：

; Simple water molecule atom types

[ atomtypes ]

; name    at.num  mass    charge  ptype   sigma      epsilon

OW        8      15.999  -0.83   A       0.316800   0.651900

HW        1      1.008   0.415   A       0.242920   0.000000

1.2. 定义分子

在拓扑文件中，我们需要定义分子的原子组成、化学键、角度、二面角等。这些信息通常在 [ molecules ] 部分中指定。

例子：定义简单水分子

; Water molecule topology

[ moleculetype ]

; name            nrexcl

water            3



[ atoms ]

; id  type  resnr  resid  atom  cgnr  charge

1    OW     1      SOL    O     1     -0.83

2    HW     1      SOL    H1    2     0.415

3    HW     1      SOL    H2    3     0.415



[ bonds ]

; id1  id2  funct  c0      c1

1     2    1      0.100   462750.0

1     3    1      0.100   462750.0



[ angles ]

; id1  id2  id3  funct  c0      c1

2     1    3    1      104.52  405.12

1.3. 定义力场参数

力场参数包括化学键、角度、二面角等的参数。这些参数通常在力场文件中定义，例如 gromos54a7.ff 力场中的 bondtypes.itp、angletypes.itp 等文件。

例子：定义简单的力场参数

; Force field parameters for simple water model

[ bondtypes ]

; i      j      func  b0      kb

OW      HW     1      0.100   462750.0



[ angletypes ]

; i      j      k      func  theta0   ktheta

HW      OW     HW     1      104.52   405.12

2. 使用工具生成拓扑文件

GROMACS 提供了一些工具来帮助用户生成和验证拓扑文件，例如 gmx pdb2gmx 和 gmx make_top。

2.1. 使用 gmx pdb2gmx 生成拓扑文件

gmx pdb2gmx 是一个强大的工具，可以从 PDB 文件生成 GROMACS 格式的拓扑文件。它会自动选择合适的力场，并生成相应的 .gro 和 .top 文件。

例子：从 PDB 文件生成水分子的拓扑文件

假设我们有一个包含水分子的 PDB 文件 water.pdb，可以使用以下命令生成拓扑文件：

gmx pdb2gmx -f water.pdb -o water.gro -p water.top -water tip3p

解释：

• -f water.pdb：输入的 PDB 文件。

• -o water.gro：输出的坐标文件。

• -p water.top：输出的拓扑文件。

• -water tip3p：选择水模型为 TIP3P。

2.2. 使用 gmx make_top 生成拓扑文件

gmx make_top 是一个更底层的工具，可以用于生成更复杂的拓扑文件。它需要用户手动指定分子的原子类型、化学键等信息。

例子：从手动输入生成水分子的拓扑文件

假设我们有一个包含水分子的坐标文件 water.crd，可以使用以下命令生成拓扑文件：

gmx make_top -f water.crd -o water.top -ff gromos54a7

解释：

• -f water.crd：输入的坐标文件。

• -o water.top：输出的拓扑文件。

• -ff gromos54a7：选择力场为 GROMOS54A7。

3. 验证拓扑文件

生成拓扑文件后，验证其正确性是非常重要的。GROMACS 提供了 gmx genconf 和 gmx editconf 等工具来帮助用户检查和修改拓扑文件。

3.1. 使用 gmx genconf 生成多个分子

gmx genconf 可以用于生成多个相同的分子，并检查它们的拓扑信息。

例子：生成 100 个水分子

假设我们已经生成了一个包含单个水分子的拓扑文件 water.top 和坐标文件 water.gro，可以使用以下命令生成 100 个水分子：

gmx genconf -f water.gro -nbox 10 10 10 -o 100waters.gro

解释：

• -f water.gro：输入的坐标文件。

• -nbox 10 10 10：在 x、y、z 方向上各生成 10 个分子。

• -o 100waters.gro：输出的坐标文件。

3.2. 使用 gmx editconf 检查和修改拓扑文件

gmx editconf 可以用于检查和修改分子的坐标和拓扑信息。例如，可以用来确认生成的多个分子是否正确排列。

例子：检查生成的 100 个水分子

gmx editconf -f 100waters.gro -o 100waters_box.gro -box 10 10 10

解释：

• -f 100waters.gro：输入的坐标文件。

• -o 100waters_box.gro：输出的坐标文件。

• -box 10 10 10：设置模拟盒子的大小为 10 nm × 10 nm × 10 nm。

4. 添加溶剂和离子

在构建复杂的分子系统时，添加溶剂和离子是常见的步骤。GROMACS 提供了 gmx solvate 和 gmx grompp 等工具来帮助用户完成这一任务。

4.1. 使用 gmx solvate 添加溶剂

gmx solvate 可以将溶剂分子添加到现有的分子系统中，确保系统充分溶剂化。

例子：将水分子添加到蛋白质系统中

假设我们有一个蛋白质系统的坐标文件 protein.gro 和拓扑文件 protein.top，可以使用以下命令添加水分子：

gmx solvate -cp protein.gro -cs spc216.gro -o protein_solvated.gro -p protein.top

解释：

• -cp protein.gro：输入的蛋白质坐标文件。

• -cs spc216.gro：溶剂分子的坐标文件。

• -o protein_solvated.gro：输出的溶剂化后的坐标文件。

• -p protein.top：输入的拓扑文件。

4.2. 使用 gmx grompp 添加离子

gmx grompp 是 GROMACS 中用于预处理模拟输入文件的工具。它可以用来添加离子，以中和系统的电荷或调整离子浓度。

例子：添加钠离子和氯离子

假设我们已经生成了一个溶剂化的系统文件 protein_solvated.gro 和拓扑文件 protein_solvated.top，可以使用以下命令添加钠离子和氯离子：

gmx grompp -f ions.mdp -c protein_solvated.gro -p protein_solvated.top -o ions.tpr

gmx genion -s ions.tpr -o protein_solvated_ions.gro -p protein_solvated.top -pname NA -nname CL -neutral -conc 0.15

解释：

• gmx grompp -f ions.mdp -c protein_solvated.gro -p protein_solvated.top -o ions.tpr：生成包含离子的 TPR 文件。

• gmx genion -s ions.tpr -o protein_solvated_ions.gro -p protein_solvated.top -pname NA -nname CL -neutral -conc 0.15：从 TPR 文件中生成含离子的坐标文件。

  • -s ions.tpr：输入的 TPR 文件。

  • -o protein_solvated_ions.gro：输出的含离子的坐标文件。

  • -p protein_solvated.top：输入的拓扑文件。

  • -pname NA：正离子的名称。

  • -nname CL：负离子的名称。

  • -neutral：确保系统电中性。

  • -conc 0.15：设置离子浓度为 0.15 M。

5. 优化初始结构

在进行分子动力学模拟之前，优化初始结构是非常重要的步骤。这可以通过能量最小化来实现，确保系统没有不合理的高能构象。

5.1. 能量最小化

能量最小化可以通过 gmx grompp 和 gmx mdrun 来实现。通常使用 em.mdp 文件来指定能量最小化的参数。

例子：进行能量最小化

假设我们已经生成了一个含离子的系统文件 protein_solvated_ions.gro 和拓扑文件 protein_solvated_ions.top，可以使用以下命令进行能量最小化：

1. 创建能量最小化输入文件 em.mdp：

; Energy minimization parameters

integrator  = steep

emtol       = 1000.0

emstep      = 0.01

nsteps      = 50000

2. 生成 TPR 文件：

gmx grompp -f em.mdp -c protein_solvated_ions.gro -p protein_solvated_ions.top -o em.tpr

3. 运行能量最小化：

gmx mdrun -v -deffnm em

解释：

• gmx grompp -f em.mdp -c protein_solvated_ions.gro -p protein_solvated_ions.top -o em.tpr：生成能量最小化的 TPR 文件。

• gmx mdrun -v -deffnm em：运行能量最小化。

6. 分子动力学模拟的初始准备

在进行分子动力学模拟之前，还需要进行一些初始准备，例如生成合适的初始速度、设置温度和压力等。

6.1. 生成初始速度

初始速度可以通过 gmx gen velocities 来生成。这一步通常在系统能量最小化后进行，确保系统从合理的初始状态开始模拟。

例子：生成初始速度

假设我们已经生成了一个能量最小化后的系统文件 em.gro 和拓扑文件 em.top，可以使用以下命令生成初始速度：

gmx gen velocities -f em.gro -o em_vel.gro -temp 300

解释：

• -f em.gro：输入的能量最小化后的坐标文件。

• -o em_vel.gro：输出的包含初始速度的坐标文件。

• -temp 300：设置初始温度为 300 K。

6.2. 设置温度和压力

设置系统的温度和压力可以通过 NVT 和 NPT 模拟来实现。NVT 模拟保持系统的体积和温度恒定，而 NPT 模拟保持系统的压力和温度恒定。

例子：进行 NVT 模拟

1. 创建 NVT 模拟输入文件 nvt.mdp：

; NVT simulation parameters

integrator  = md

dt          = 0.002

nsteps      = 50000

tinit       = 0

init_t      = 0

tcoupl      = v-rescale

tc_grps     = system

tau_t       = 1.0

ref_t       = 300

2. 生成 TPR 文件：

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

3. 运行 NVT 模拟：

gmx mdrun -v -deffnm nvt

解释：

• gmx grompp -f nvt.mdp -c em_vel.gro -p em.top -o nvt.tpr：生成 NVT 模拟的 TPR 文件。

• gmx mdrun -v -deffnm nvt：运行 NVT 模拟。

例子：进行 NPT 模拟

1. 创建 NPT 模拟输入文件 npt.mdp：

; NPT simulation parameters

integrator  = md

dt          = 0.002

nsteps      = 50000

tinit       = 0

init_t      = 0

tcoupl      = v-rescale

tc_grps     = system

tau_t       = 1.0

ref_t       = 300

pcoupl      = Parrinello-Rahman

pcoupltype  = isotropic

tau_p       = 2.0

ref_p       = 1.0

compressibility = 4.5e-5

2. 生成 TPR 文件：

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

3. 运行 NPT 模拟：

gmx mdrun -v -deffnm npt

解释：

• gmx grompp -f npt.mdp -c nvt.gro -p em.top -o npt.tpr：生成 NPT 模拟的 TPR 文件。

• gmx mdrun -v -deffnm npt：运行 NPT 模拟。

7. 调整和优化拓扑文件

有时需要对生成的拓扑文件进行调整和优化，以确保系统的稳定性和准确性。GROMACS 提供了多种工具来帮助用户完成这一任务。

7.1. 使用 gmx check 检查拓扑文件

gmx check 可以用于检查拓扑文件的完整性，确保所有必需的信息都已正确指定。

例子：检查拓扑文件

gmx check -f em.tpr

解释：

• -f em.tpr：输入的 TPR 文件。

7.2. 使用 gmx analyse 分析系统属性

gmx analyse 可以用于分析系统的各种属性，例如能量、温度、压力等，确保系统在模拟过程中表现正常。

例子：分析能量最小化结果

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

解释：

• -f em.edr：输入的能量文件。

• -o em_energy.xvg：输出的能量分析文件。

7.3. 使用 gmx trjconv 转换轨迹文件

gmx trjconv 可以用于转换轨迹文件，例如将 PDB 文件转换为 GRO 文件，或将 GRO 文件转换为其他格式。

例子：转换轨迹文件

假设我们有一个轨迹文件 md.xtc 和对应的坐标文件 md.gro，可以使用以下命令将其转换为 PDB 文件：

gmx trjconv -f md.xtc -s md.gro -o md.pdb -pbc mol

解释：

• -f md.xtc：输入的轨迹文件。

• -s md.gro：输入的坐标文件。

• -o md.pdb：输出的 PDB 文件。

• -pbc mol：设置周期性边界条件为分子。

8. 高级拓扑文件调整

对于更复杂的分子系统，可能需要进行更高级的拓扑文件调整，例如添加自定义的力场参数、定义特殊的相互作用等。这些调整有助于确保系统的稳定性和准确性。

8.1. 添加自定义的力场参数

有时力场文件中没有包含某些特定分子的参数，用户需要手动添加这些参数。这些参数可以包括原子类型、化学键、角度、二面角等的定义。

例子：添加自定义的力场参数

假设我们需要添加一个自定义的原子类型 CX，以及相应的化学键、角度和二面角参数。可以按照以下步骤进行：

1. 编辑力场文件 ffnonbonded.itp：

; Custom atom type

[ atomtypes ]

; name    at.num  mass    charge  ptype   sigma      epsilon

CX        6      12.011  0.00    A       0.350      0.100

2. 编辑力场文件 bondtypes.itp：

; Custom bond parameters

[ bondtypes ]

; i      j      func  b0      kb

CX      CX     1      0.150   500000.0

3. 编辑力场文件 angletypes.itp：

; Custom angle parameters

[ angletypes ]

; i      j      k      func  theta0   ktheta

CX      CX     CX     1      110.00   400.0

4. 编辑力场文件 dihedraltypes.itp：

; Custom dihedral parameters

[ dihedraltypes ]

; i      j      k      l      func  phi0   kphi

CX      CX     CX     CX     1      0.00   100.0

5. 在拓扑文件中引用自定义的力场参数：

在拓扑文件 topol.top 中，需要包含自定义的力场参数文件：

; Include force field parameters

#include "ffnonbonded.itp"

#include "bondtypes.itp"

#include "angletypes.itp"

#include "dihedraltypes.itp"



; Molecule type

[ moleculetype ]

; name            nrexcl

custom_molecule  3



[ atoms ]

; id  type  resnr  resid  atom  cgnr  charge

1    CX     1      RES    C1    1     0.0

2    CX     1      RES    C2    2     0.0

3    CX     1      RES    C3    3     0.0



[ bonds ]

; id1  id2  funct  c0      c1

1     2    1      0.150   500000.0

2     3    1      0.150   500000.0



[ angles ]

; id1  id2  id3  funct  c0      c1

1     2     3    1      110.00  400.0



[ dihedrals ]

; id1  id2  id3  id4  funct  c0      c1

1     2     3     4    1      0.00   100.0

8.2. 定义特殊的相互作用

在某些情况下，可能需要定义特殊的相互作用，例如特定原子之间的排斥力或吸引力。这些特殊的相互作用可以通过在拓扑文件中添加 dihedrals、pairs 或 settles 等节来实现。

例子：定义特殊的排斥力

假设我们需要在两个特定的原子之间定义一个特殊的排斥力，可以使用 pairs 节来实现：

1. 在拓扑文件中添加 pairs 节：

; Molecule type

[ moleculetype ]

; name            nrexcl

custom_molecule  3



[ atoms ]

; id  type  resnr  resid  atom  cgnr  charge

1    CX     1      RES    C1    1     0.0

2    CX     1      RES    C2    2     0.0

3    CX     1      RES    C3    3     0.0



[ bonds ]

; id1  id2  funct  c0      c1

1     2    1      0.150   500000.0

2     3    1      0.150   500000.0



[ angles ]

; id1  id2  id3  funct  c0      c1

1     2     3    1      110.00  400.0



[ dihedrals ]

; id1  id2  id3  id4  funct  c0      c1

1     2     3     4    1      0.00   100.0



[ pairs ]

; id1  id2  funct  c0      c1

1     3    1      0.350   100000.0

2. 解释 pairs 节：

   • 1 3 1 0.350 100000.0：定义了原子 1 和原子 3 之间的特殊相互作用，使用 Lennard-Jones 势能函数（funct 1），并设置了相应的参数 c0 和 c1。

8.3. 使用 gmx make_ndx 定义索引组

索引组可以帮助用户在模拟过程中更方便地选择和分析特定的分子或原子。gmx make_ndx 是一个交互式工具，可以用来定义这些索引组。

例子：定义索引组

假设我们有一个包含蛋白质和水分子的系统文件 protein_solvated.gro，可以使用以下命令定义索引组：

gmx make_ndx -f protein_solvated.gro

在交互模式中，可以输入以下命令来定义索引组：

13 & 14  ; 选择蛋白质和水分子的索引

name 13 Protein  ; 将索引 13 命名为 Protein

name 14 Water     ; 将索引 14 命名为 Water

q             ; 退出交互模式

解释：

• 13 & 14：选择索引 13 和 14 的分子。

• name 13 Protein：将索引 13 命名为 Protein。

• name 14 Water：将索引 14 命名为 Water。

• q：退出交互模式。

8.4. 调整周期性边界条件

周期性边界条件（PBC）对于确保分子系统在模拟过程中保持正确的物理特性非常重要。GROMACS 提供了多种 PBC 设置方法，用户可以根据需要进行调整。

例子：调整周期性边界条件

假设我们有一个系统文件 protein_solvated.gro，可以使用以下命令调整周期性边界条件：

gmx editconf -f protein_solvated.gro -o protein_solvated_pbc.gro -bt dodecahedron -box 10 10 10

解释：

• -f protein_solvated.gro：输入的坐标文件。

• -o protein_solvated_pbc.gro：输出的调整后的坐标文件。

• -bt dodecahedron：设置周期性边界条件为十二面体。

• -box 10 10 10：设置模拟盒子的大小为 10 nm × 10 nm × 10 nm。

9. 总结

构建一个准确的初始拓扑模型是分子动力学模拟成功的关键。通过定义原子类型、化学键、角度、二面角等信息，以及使用 GROMACS 提供的工具生成和验证拓扑文件，可以确保模拟系统的稳定性和准确性。此外，对于更复杂的系统，进行高级的拓扑文件调整和优化也是非常必要的。希望本文档能帮助用户更好地理解和使用 GROMACS 来构建和优化初始拓扑模型。