GROMACS教程：创建周期性体系的拓扑文件：以石墨烯为例

• 概述
• 1. 获取石墨烯的结构文件
• 2. 确定使用的力场
• 3. 准备n2t文件

• 2016-02-16 整理: 阮洋; 修订: 李继存
• 2017-01-24 补充: 陈建发

概述

有许多材料具有周期性结构, 下面是一些常见的例子.

还有其它一些类似的材料, 如SiO₂纳米管, BN纳米管, SiC纳米管以及石墨炔材料等, 都具有周期性结构. 在使用GROMACS对这些材料进行模拟时, 我们需要创建这些材料的拓扑文件, 也就是力场文件. 为此, GROMACS提供了一个根据坐标构建拓扑文件的工具g_x2top(对GROMACS 5.0以上版本, 使用gmx x2top), 非常适用于创建这些结构规律性很强的材料的拓扑文件. 下面我们以石墨烯为例来说明如何使用g_x2top来创建材料的拓扑文件并进行简单的模拟.

1. 获取石墨烯的结构文件

石墨烯是二维平面结构, 可认为其碳原子属于sp2杂化, 有关其结构的一些特点请参看石墨烯：建模, 几何性质及力场模拟.

创建石墨烯结构的方法有很多. 这里我们使用一个石墨烯在线创建工具来创建. 根据你的需要, 创建合适大小的石墨烯. 此工具给出的坐标是xyz格式的, 为了后面便于处理, 我们需要将其改写为GROMACS的gro格式. 在这里我不建议使用pdb格式, 因为pdb文件格式有不同版本, 而不同版本的g_x2top程序所支持的pdb文件格式稍有不同, 后面处理时可能出现比较奇怪的问题.

在保存为gro文件时, 需要非常注意的一点就是盒子的大小. 石墨烯是二维周期性的无限体系, 所以要考虑其周期性话, 盒子的长宽需要精确地等于实际的长宽, 否则在后面创建具有周期性的拓扑文件时就会出错. 如果使用pdb文件同样也需要注意这一点. 你可以使用VMD等软件对得到的构型进行周期性扩展, 检查盒子大小是否合适. 如果不考虑周期性, 只是将石墨烯视为孤立的分子, 那么盒子的大小只要不小于实际大小即可.

利用前面的在线工具, 我们创建一个10x5的石墨烯, 其长宽为24.248711x21.000000Å. 将得到的xyz坐标转成gro格式, 并将盒子的长宽设为前面的值(注意,xyz文件的单位为Å, 而gro文件的单位为nm), 对于盒子的高, 根据你要模拟的体系大小设置即可. 我们在这里将其设置为2 nm. 同时, 为了便于查看, 我们使石墨烯处于盒子中间. 你可以点击这里下载已经准备好的文件.

下面是放于盒子中的石墨烯

2. 确定使用的力场

在模拟时, 可以使用刚性的石墨烯模型, 也可以使用柔性的模型. 如果使用刚性的模型, 在模拟过程中石墨烯的结构保持不变, 所有不需要石墨烯自身的力场参数, 只需要石墨烯和体系中其他分子之间的相互作用参数. 如果使用柔性的模型, 在模拟过程中石墨烯的结构会发生变化, 这就需要能够描述石墨烯变形的力场参数. 在这种情况下, 模拟时还可以有两种处理方式. 一种是将石墨烯视为孤立的分子, 不考虑其周期性, 这样在模拟时边界碳原子由于缺少约束, 可能会卷曲. 如果发生了卷曲, 可以通过固定一些边界碳原子来解决. 另一种处理方式是将石墨烯视为具有周期性结构的无限体系, 这样就不存在边界碳原子, 也就不会发生石墨烯卷曲的问题. 请根据你自己的情况选择使用哪种模拟模式.

我们在这里将使用柔性的石墨烯模型, 而且考虑其周期性. 这样除石墨烯与其他分子的相互作用参数外, 我们还需要石墨烯自身的相互作用参数, 也就是能够描述石墨烯变形的力场.

根据自己的需要, 确定使用哪个石墨烯力场. 一般来说, 一个石墨烯的力场至少要考虑相邻两个碳原子之间的键, 相邻三个碳原子之间的键角, 相邻四个碳原子之间的二面角, 并且二面角可能有两种类型, 一类是顺式, 一类是反式. 更复杂的力场可能还要考虑1–3相邻碳原子之间的相互作用. 这种力场的一个例子可参见这篇论文. 当然, 文献上还有其他的一些石墨烯力场, 请根据自己的情况选用.

在这里为简化起见, 我们使用OPLS-AA力场中的opls_145原子类型(苯的C原子)作为石墨烯中C原子的原型, 且只考虑碳碳键, 键角和二面角相互作用. 对一些小分子, 当没有现成的力场可用时, 选用已有力场中的类似原子类型进行构建是一个可行的方法. 当然, 这需要你去查看已有的力场文件, 并清楚里面已有的各种原子类型及其相应的参数.

3. 准备n2t文件

创建拓扑文件时需要判断原子之间的成键关系. g_x2top判断成键关系时使用了最简单的方法, 距离准则. 如果两个原子之间的距离小于或等于设定的值, 就认为这两个原子之间存在键合相互作用. 根据一个原子周围连接的键数与每个键所连接的原子类型就可以判定这个原子的类型. 判断成键关系的键长和原子周围的键数, 是由n2t文件读入的(n2t代表name to type), 所以在使用g_x2top工具前我们需要准备一个n2t文件. 为此, 我们需要知道此文件的格式.

我们先来看看OPLS-AA力场自带的n2t文件, 它位于GROMACS主目录/share/top/oplsaa.ff/atomname2type.n2t, 文件很小, 全部内容如下:

C    opls_157    -0.18   12.011   4    H 0.108   H 0.108   H 0.108   C 0.150
C    opls_158    -0.12   12.011   4    H 0.108   H 0.108   C 0.150   C 0.150
C    opls_157     0.145  12.011   4    H 0.108   O 0.141   H 0.108   C 0.150
C    opls_157     0.145  12.011   4    H 0.108   H 0.108   O 0.108   C 0.150
C    opls_158     0.205  12.011   4    H 0.108   C 0.150   O 0.140   C 0.150
C    opls_158    -0.06   12.011   4    H 0.108   C 0.150   C 0.150   C 0.150
C    opls_145    -0.12   12.011   3    C 0.150   C 0.150   H 0.108
C    opls_145    -0.12   12.011   3    C 0.133   C 0.150   O 0.140
C    opls_235     0.5    12.011   3    C 0.133   N 0.140   O 0.140
C    opls_235     0.5    12.011   3    C 0.133   N 0.132   C 0.150
C    opls_235     0.5    12.011   3    C 0.133   N 0.132   H 0.108
O    opls_236    -0.5    15.9994  1    C 0.123
N    opls_238    -0.5    14.0067  3    H 0.108   C 0.140   C 0.150
N    opls_238    -0.5    14.0067  3    H 0.108   C 0.132   N 0.123
N    opls_238     0      14.0067  2    C 0.140   N 0.123
H    opls_140     0.06    1.008   1    C 0.108
H    opls_155     0.418   1.008   1    O 0.095
H    opls_241     0.30    1.008   1    N 0.095
O    opls_154    -0.683  15.9994  2    C 0.140   H 0.095
P    opls_393     1      30.97376 4    O 0.180   O 0.180   O 0.180   O 0.180
O    opls_395    -1      15.9994  2    P 0.180   H 0.095
O    opls_394    -1      15.9994  1    P 0.180
N    opls_237     0      14.0067  4    C 0.140   C 0.140   C 0.140   C 0.140

可以看到, 每个原子一行, 第一列为构型中的原子名称, 第二列为力场中原子类型, 第三列为原子的电荷, 第四列为原子的摩尔质量, 第五列为原子周围连接的其他原子数, 后面的列则分别给出所连接的每个原子的名称以及键长. 关于此文件格式的其他说明请参考GROMACS手册.

OPLS-AA自带的n2t文件内容很少, 不能满足我们的需要, 为此, 我们需要对其进行修改, 增加石墨烯碳原子的设置. 对未修饰的石墨烯, 不考虑周期性时, 可能的碳原子类型有三种: 分别对应于与周围一个, 两个和三个原子相连接的碳原子. 前两种是边界碳原子, 最后一种是正常的碳原子. 当考虑周期性时, 就需要考虑一种碳原子类型了. 此外, 我们假定石墨烯中碳原子的电荷为0, 当不考虑边界碳原子时, 这是个很好的近似. 当然, 文献中有些模拟的石墨烯表面是带电的. 如果你也需要这样做, 可以修改相应的电荷列.

在我们的石墨烯结构中, 碳碳键长为0.142 nm. 我们希望被修改的参数，同时尽量不影响的其他模拟。比较安全的做法是将GROMACS主目录/share/top/oplsaa.ff/目录拷贝到工作目录，所有的修改都在拷贝的atomname2type.n2t文件中进行。运行GROMACS时，程序会自动加载当前目录下的force field文件. 还有一种选择，就是在GROMACS主目录/share/top/oplsaa.ff/目录中拷贝一个atomname2type.n2t文件副本，将其命名为graphene.n2t，所有的修改都在graphene.n2t中进行.

C    opls_145   0.00      12.011  1    C 0.142
C    opls_145   0.00      12.011  2    C 0.142   C 0.142
C    opls_145   0.00      12.011  3    C 0.142   C 0.142   C 0.142

这三行分别定义了三种不同的碳原子类型, 并且指定碳原子之间的键长. 其中第三行说明如果体系中的任何一个碳原子, 与周围3个碳原子的距离都在0.142 nm左右, 就把它当作opls_145原子类型, 其电荷为0, 摩尔质量为12.011. 第一行和第二行的意义类似, 即说明与一个碳原子和与两个碳原子成键的碳原子被当成什么类型, 电荷多少, 摩尔质量多少.

需要说明的是, 上面的原子类型名称也可以使用自己定义的名称, 而不是使用力场中已有的原子类型, 这种情况下, g_x2top仍然能够生成拓扑文件, 只是要对拓扑文件进行更多的修改. 理论上只要明白拓扑文件的含义, 再借助g_x2top就可以创建需要的拓扑文件了.

运行下面的命令

g_x2top -f gra.gro -o gra.top -ff oplsaa -pbc -name graphene -kb 400000 -kt 600 -kp 150

上面命令中-ff选项指明了使用的力场. -pbc选项指明计算成键时考虑周期性边界条件. 如果你将石墨烯视为孤立分子, 请使用-nopbc选项. -name选项指定了拓扑文件中分子的名称. 默认情况下, g_x2top会自动为键合相互作用加上力场参数, 但这种自动加入的力场参数并不是根据力场中相应原子类型间的键合参数加入的. 所以, 我们使用了选项-kb, -kp, -kt, 它们分别指定键, 键角和二面角的力常数, 它们的默认值为400000, 400, 5. 如果不想让程序自动加入力常数, 可以使用-noparam选项. 关于g_x2top的选项, 请参考其文档.(【陈建发 注】GROMACS 5.0.X版本存在bug, 运行会导致segmentation fault, 其他版本运行正常. 建议不要使用5.0系列版本.)

运行命令后, 屏幕上会输出类似下面的内容:

Opening force field file C:/Users/Jicun/E_Prf/Gromacs-4.6/share/top/oplsaa.ff\atomname2type.n2t
Opening force field file C:/Users/Jicun/E_Prf/Gromacs-4.6/share/top/oplsaa.ff\graphene.n2t
There are 26 name to type translations in file C:/Users/Jicun/E_Prf/Gromacs-4.6/share/top/oplsaa.ff
Generating bonds from distances...
atom 200
There are 1 different atom types in your sample
Generating angles and dihedrals from bonds...
Before cleaning: 1200 pairs
Before cleaning: 1200 dihedrals
There are  300 Ryckaert-Bellemans dihedrals,    0 impropers,  600 angles
           900 pairs,      300 bonds and   200 atoms
Total charge is 0, total mass is 2402.2

WARNING: topologies generated by g_x2top can not be trusted at face value.
         Please verify atomtypes and charges by comparison to other
         topologies.

可以看到, g_x2top在确定原子类型时, 使用了OPLS-AA自带的atomname2type.n2t文件, 也使用了我们创建的graphene.n2t文件. 对含200个碳原子的石墨烯, 程序计算得到键数300, 键角数600, 二面角数1200. 这些值与理论值相等(当碳原子是为N时, 若考虑周期性, 则键数为1.5N, 键角数3N, 二面角数6N), 这说明, g_x2top正确地确定了原子的类型和成键情况.

到这一步, 我们已经得到了一个比较原始的石墨烯的拓扑文件gra.top(下载), 里面列出了所有的键合相互作用.

在使用之前, 我们先检查一下前面得到的拓扑文件, 看看其中的信息是否完整, 需不需要修改.

打开得到的拓扑文件, 如果你熟悉GROMACS拓扑文件的格式, 就会发现, 里面列出了所有的键, 键角, 和二面角, 并根据我们的设定给出了势函数的参数. 你可以在这个拓扑文件的基础之上进行修改, 以满足自己的需要. 例如, 我们使用了最简单的简谐势来描述键, 键角和二面角, 因此只需要给出力常数即可. 如果你要使用更复杂的势能函数, 只有力常数是不够的, 这就需要你对相应的相互作用项进行修改.

在文件中, [ bonds ]部分内容类似下面这样:

1     2     1 1.400000e-001 4.000000e+005 1.400000e-001 4.000000e+005
1    20     1 1.400000e-001 4.000000e+005 1.400000e-001 4.000000e+005
1   182     1 1.410000e-001 4.000000e+005 1.410000e-001 4.000000e+005
2     3     1 1.400000e-001 4.000000e+005 1.400000e-001 4.000000e+005
2    21     1 1.400000e-001 4.000000e+005 1.400000e-001 4.000000e+005
3     4     1 1.400000e-001 4.000000e+005 1.400000e-001 4.000000e+005
... ...

前两列为原子编号, 第三列为势函数类型, 第四列为碳碳的平衡键长, 第五列为力常数, 与我们设定的值一致, 最后两列对我们使用的势函数是多余的, 可以忽略或直接删除. [ angles ]和[ dihedrals ]部分的内容与此类似.

上面使用的力常数在某些情况下并不重要. 如果你只是用它们来维持石墨烯的近似刚性结构, 只要将它们设置为较大的值就可以了, 但值太大时可能会导致振动频率过高, 致使模拟不稳定, 或需要使用很小的时间步长. 如果你想使石墨烯更刚性(柔性), 只要增大(减小)力常数即可. 值得注意的是, 如果对你的模拟而言, 石墨烯的力学性质非常重要, 那你最好不要随意设置这些力常数, 而是需要使用能给出好的力学性质的石墨烯力场. 如果没有合适的, 可以试试前面提到个力场, 它能给出与实验值比较符合的力学性质和导热系数.

如果你不需要考虑石墨烯中碳原子之间的范德华相互作用, 可以将[ moleculetype ]部分的nrexcl改为1. 默认使用的3会考虑相邻3条键以上的两个碳原子之间的范德华相互作用.

如果你想要将碳原子设置为力场中没有的原子类型, 并需要加入特殊的力场参数, 则需要定义新的原子类型以及相互作用参数. 首先在[ atomtypes ]部分定义新的原子类型, 然后在 [ nonbond_params ]里面添加范德华相互作用参数. 如果前面运行g_x2top时使用了-noparam, 还需要设定键合参数. 可使用下面几种方法:

1. 直接在拓扑文件中每一个键合项的后面加入相应力场参数, 类似于g_x2top自动在添加力场参数的方式. 使用支持列编辑的文本编辑工具, 比较容易完成.
2. 不在拓扑文件中添加, 而是在ffbonded.itp文件中的[ bondtypes ], [ angletypes ]和[ dihedraltypes ]中加入相应的项. 这种方法会对力场文件进行修改, 不建议.
3. 将上面的内容单独写进一个.itp文件, 之后#include进拓扑文件中. 这种方法比较省事也易于修改.
4. g_x2top可以生成一个.rtp文件, 如果将.rtp文件加入到力场中的.rtp数据库中, 就可以利用pdb2gmx命令生成拓扑文件. 这需要原子类型存在于.atp文件中, 否则使用pdb2gmx命令会出错. 这种方法需要熟悉GROMACS的力场处理方式, 只适合有经验的用户.

得到了拓扑文件后, 我们再创建一个vac.mdp文件(下载), 并搭配前面的gro文件就可以进行真空中的模拟了. 运行下面的命令进行模拟:

grompp -f vac.mdp -c gra.gro -p gra.top -o vac.tpr
mdrun -v -deffnm vac

如果你要填充溶剂分子或其他分子, 请参考其他的教程. 下面是一个添加水分子进行水溶液中模拟的简单示例.

• 先填充水分子: genbox -cp gra.gro -cs spc216.gro -o gra_wat.gro, (下载gra_wat.gro)

• 再将拓扑文件gra.top另存为gra_wat.top并进行修改: 在文件的开始处包含使用的水分子模型#include "oplsaa.ff/spce.itp", 在文件的结束处添加填充水分子的名称和数目SOL 252, (下载gra_wat.gro)

• 然后创建文件nvt.mdp (下载)

• 最后执行下面的命令:

  grompp -f nvt.mdp -c gra_wat.gro -p gra.top -o nvt.tpr
  mdrun -v -deffnm nvt
  

本文只对创建石墨烯的拓扑文件进行了介绍, 按照这个思路, 只要有了确定的力场参数以及电荷信息, 同样可以创建BN, SiC, CNT, SiN, PNT等材料的拓扑文件. 对于有机小分子有时也可以使用这种方法.

1. 使用GROMACS模拟碳纳米管的一个简单例子
2. Amber与Gromacs读入碳纳米管的方法
3. GROMACS Carbon Nanotube
4. Modeling Carbon Nanotubes with GROMACS
5. help on Graphene Nano Sheets
6. Single Wall Carbon Nanotubes in 4.0.3 GROMACS