原文：James (Wes) Barnett
翻译：谷歌/康文渊-湖南大学

在本教程中，我将向您展示如何在一个盒子的TIP4PEW水中创建一个包含一个OPLS力场甲烷的系统。

##设置
和以前一样，我们需要一个结构文件，一个拓扑文件和参数文件。我们将使用GROMACS工具gmx pdb2gmx从pdb文件生成拓扑。

用pdb2gmx设置残基

对于这个分子，我们将使用OPLS力场。 力场位于顶层力场目录（可能是/usr/share/gromacs/top或类似的东西）。

如果您不确定您的GROMACS安装位置在哪里，可以使用如下命令：

1
$ echo $ GMXPREFIX

如果您sourcing 了GROMACS配置文件，则会为您提供安装位置。在该目录中找到share/gromacs/top并进入它（例如，如果GMXPREFIX是/usr，则转到/usr/share/gromacs/top）。或者你可以简单地转到$GMXDATA/top。

我们来看看force field目录的内容：

1
2
$ cd oplsaa.ff
$ ls

你会看到几个文件，但我们现在只对其中的一些文件感兴趣。注意forcefield.itp。这是模拟中使用的主要文件。在里面你会看到一个[defaults]部分以及包含两个其他文件 - 一个用于键联系，另一个用于非键联系。我们也对atomtypes.atp感兴趣，它给出了难以理解的opls _ ####术语的描述以及aminoacids.rtp，它们给出了用于gmx pdb2gmx命令能够识别的氨基酸残基列表。

用你的文本编辑器打开atomtypes.atp,比如用vim打开它：

1
$ vim atomtypes.atp

转到opls_138的行。请注意，注释是否为alkane CH4。但是，请注意第二列中的mass - 这只是CH4组的碳，所以我们也需要氢。这是一个“全原子”模型 - 每个原子都被表示出来。相应的氢是opls_140。您可能还需要查看OPLS force field 论文的支持信息。论文中的参数应该与我们刚刚看到的参数相匹配。现在记下这两种原子类型并关闭文件。

让我们来看看这两种原子类型的ffnonbonded.itp：

1
2
$ grep opls_138 ffnonbonded.itp
$ grep opls_140 ffnonbonded.itp

在这里，我们看到原子类型的名称，键类型，质量，电荷，ptype，sigma和epsilon。记下每个类型的电荷 - 我们将需要它来构建我们新的氨基酸残基。作为一个方面说明，ffbonded.itp将使用键类型，键类型和二面角类型。

在继续之前，您可能希望将您的顶级力场目录复制到别的地方，例如您的主目录，因为我们将修改它并添加一些文件。将其复制到您的主目录中：

1
$ cp -r $GMXDATA/top $HOME/GMXLIB

你可能需要root权限执行它。现在将$ GMXLIB环境变量更改为：

1
$ export GMXLIB = $ HOME/GMXLIB

将上述内容添加到.bash_profile中以使其成为永久性的,然后进行执行：

1
$ cd $ GMXLIB

您现在处于刚刚拷贝的副本中，所有模拟将使用该目录而不是GROMACS默认目录中提供的目录。

现在进入oplsaa.ff并打开aminoacids.rtp。您会注意到文件中已有多个残基。我们将为我们的甲烷添加一个名为methane.rtp的新文件，其中包含我们称之为CH4的残基。关闭aminoacids.rtp。我们需要告诉gmx pdb2gmx我们的残差文件中原子的原子和键。我们也可以告诉它角度，但是我们会将它们排除在外，因为gmx pdb2gmx会为我们解决这个问题。现在需要在oplsaa.ff目录中创建以下内容并另存为methane.rtp：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
[ bondedtypes ]
; bonds  angles  dihedrals  impropers all_dihedrals nrexcl HH14 RemoveDih
  1      1       3          1         1             3      1    0
; Methane
[ CH4 ]
 [ atoms ]
   C      opls_138    -0.240     1
   H1     opls_140     0.060     1
   H2     opls_140     0.060     1
   H3     opls_140     0.060     1
   H4     opls_140     0.060     1
 [ bonds ]
   C      H1
   C      H2
   C      H3
   C      H4

关于上述文件的一些注意事项：[bondedtypes]来自aminoacids.rtp并且是必需的。 在[atoms]下的名字可以随便命名，但是需要匹配稍后构建pdb文件中的名字。 注意在第一列中我们给出了原子名称，然后我们给出了原子类型，电荷，然后是电荷组。 在[bonds]下，我们只是告诉它每个原子如何连接到其他原子。 在这种情况下，C与每个氢都有连接。我们可以选择添加[angles]，但是如前所述，GROMACS会为我们搞定。 现在关闭文件。 有关这方面的更多信息，请参阅第5.6节。

创建pdb文件并且运行gmx pdb2gmx

现在我们准备创建pdb文件。 有几个程序可以创建分子结构文件。例如 Avogadro。 另一种方法是使用AmberTools软件包中的“xleap”。 将这个文件保存为methane.pdb。 你的文件应该看起来像这样。 将其保存在您的主目录中的某个位置，但不在$ GMXLIB中的任何位置。

在methane.pdb中将LIG更改为CH4。还要将第一个H更改为H1，将第二个更改为H2等等。PDB文件是固定格式，因此请将每列的开头保留在相同的位置。CONNECT和MASTER记录也不需要，因此可以删除它们。同时继续并将UNNAMED更改为METHANE。你修改过的文件应该是这样的：

1
2
3
4
5
6
7
8
COMPND    METHANE
AUTHOR    GENERATED BY OPEN BABEL 2.3.2
HETATM    1  C   CH4     1      -0.370   0.900   0.000  1.00  0.00           C  
HETATM    2  H1  CH4     1       0.700   0.900   0.000  1.00  0.00           H  
HETATM    3  H2  CH4     1      -0.727   0.122   0.643  1.00  0.00           H  
HETATM    4  H3  CH4     1      -0.727   0.731  -0.995  1.00  0.00           H  
HETATM    5  H4  CH4     1      -0.727   1.845   0.351  1.00  0.00           H  
END

保存文件为methane.pdb.

现在我们可以使用gmx pdb2gmx创建GROMACS .conf和.top文件：

1
$ gmx pdb2gmx -f methane.pdb

系统会提示您选择力场。选择OPLS。如果您在两个不同的力场目录之间有一个选项，请选择您所复制目录中的OPLS。对于水模型选择TIP4PEW。如果您发现GROMACS找不到残基CH4的错误，您可能会使用错误的力场。

将创建三个文件：conf.gro，posre.itp和topol.top。conf.gro是我们的文件，只包含一个甲烷，topol.top是系统的拓扑文件，posre.itp是我们溶质（甲烷）的可选位置限制文件。我们不会使用那个。在topol.top文件中注意到有一个[angles]部分。您还需要在topol.top中重命名化合物。看看并探索每个文件。 GROMACS手册的第5章将帮助您更多地了解拓扑文件。

注意：topol.top和methane.pdb将在其他教程中再次使用。

对于那些使用gmx pdb2gmx生成大型蛋白质拓扑的人来说，事情会变得更加复杂。这仅仅是一个简单的例子，我们可能真的可以在其他地方找到这种拓扑。

溶剂体系

我们的结构文件和拓扑文件到目前为止只有我们的甲烷。我们需要通过使用gmx solvate来添加水：

1
$ gmx solvate -cp conf.gro -o conf.gro -cs tip4p -p topol.top -box 2.3 2.3 2.3

参数文件

我们将使用前一教程中的相同文件。

模拟

我们将使用与上次相同的流程。这假定您的mdp文件位于名为mdp的目录中：

1
2
3
4
5
6
7
8
9
10
$ gmx grompp -f mdp/min.mdp -o min -pp min -po min
$ gmx mdrun -deffnm min
$ gmx grompp -f mdp/min2.mdp -o min2 -pp min2 -po min2 -c min -t min
$ gmx mdrun -deffnm min2
$ gmx grompp -f mdp/eql.mdp -o eql -pp eql -po eql -c min2 -t min2
$ gmx mdrun -deffnm eql
$ gmx grompp -f mdp/eql2.mdp -o eql2 -pp eql2 -po eql2 -c eql -t eql
$ gmx mdrun -deffnm eql2
$ gmx grompp -f mdp/prd.mdp -o prd -pp prd -po prd -c eql2 -t eql2
$ gmx mdrun -deffnm prd

您可以把它当作脚本来使用，在头部添加如下内容:

1
2
3
#!/bin/bash
set -e

然后执行

1
$ chmod +x run

运行脚本:

1
$ ./run

其中set -e表示如果有错误立马停止

分析

我们来计算一下称为径向分布函数的东西。首先，我们需要创建一个索引文件：

1
$ gmx make_ndx -f conf.gro
1
2
3
> a C
> a OW
> q

然后运行 gmx rdf

1
$ gmx rdf -f prd.xtc -n index.ndx

在提示时选择C作为参考组。 然后选择OW。 然后输入CTRL-D结束。 数据的第1列和第3列中的结果图将通过在gnuplot中执行以下操作来绘制：

1
> plot 'rdf.xvg' u 1:3 w l

它应该看起来像这样：

总结

在本教程中，我们学习了如何创建一个用于gmx pdb2gmx的残基模板文件（.rtp）。 我们为OPLS甲烷创建了一个结构，并为其生成了一个拓扑。我们用gmx solvated对它周围进行了加水。在此之后，我们进行了仿真，就像上次一样。最后，我们使用gmx rdf发现了C-OW径向分布函数。

问题

如果教程中存在错误或错误，或者如果有什么不清楚的地方，请打开一个问题，我很乐意解决它。 如果您有关于GROMACS的一般问题，请参阅GROMACS文档。 您还可以将关于GROMACS的一般问题通过查阅邮件列表来解决。