分子动力学模拟三聚氰胺溶液性质

简单介绍了如何利用GROMACS+acpype模拟三聚氰胺水溶液

分子结构准备

对于有机小分子,直接在ChemSpider就好了,下载3D结构的mol格式文件。然后用任意分子可视化软件(比如avogadro)打开分子结构文件,检查结构是否完整和正确,并另存为mol2文件。

acpype使用

目前有两种方式,一种是在本机上利用acpype脚本(当然需要安装ambertools工具)直接产生适用于GROMACS的OPLSAA力场的拓扑文件(实验性质的);另一种方式是在acpype在线网站上提交。两者目的相同。我这里直接在线提交mol2结构文件。如下图:

在这里我选择的是直接利用antechamber产生的bcc电荷方式,整个分子净电荷为0,原子类型使用GAFF力场对应的原子。点击提交,对于这种小分子很快会得到结果。下载文件压缩包,将里面的如下几个文件拉出来就可以了:

GROMACS模拟

溶液体系构建

首先得构建一个含有三聚氰胺分子的溶液体系,这里使用了5个三聚氰胺分子,填充盒子大小为5.5 x 5.5 x 5.5 $nm^3$:

  • 第一步,将5个目标分子置于盒子中:

    1
    gmx insert-molecules -ci melamine_GMX.gro -nmol 5 -box 5.5 5.5 5.5 -o box.gro
  • 第二步,填满水分子。(本例采用SPC/E水模型,因此这里用spc216.gro

    1
    gmx solvate -cp box.gro -cs spc216.gro -o sol.gro

共添加水分子5294个,体系总原子15957。
体系如下图所示:

  • 第三步,处理一下体系拓扑。
    因为acpype处理以后会得到itptop文件,其中itp中包含的是三聚氰胺分子原子信息、电荷等拓扑文件,而top文件是整个体系的拓扑文件,因此我们应该在top中添加上水分子的相关部分,因为我整个体系想用OPLS-AA力场模拟,因此不得不注释掉它自己参数的那个[ defaults ]部分,因为这个头部是用于Amber力场的。整个改动如下图:

模拟开始

包括能量极小化,预平衡和成品模拟。在这里我提供了一次性完成整个流程的bash脚本,点击下载

分析步骤

  • 比如分析三聚氰胺氮原子周围水分子的氢原子径向分布情况,可以通过gmx rdf命令去完成。下面给出了本例中计算的rdf与文献中的结果比对情况图: 由于本例子模拟时间很短,因此并不与论文完全相符,但基本趋势相同。
    详略。可以看原文献 A molecular dynamics study on the melamine aqueous solution

注意事项

  • 这里的例子只是起说明作用,特别是通过acpype产生的参数用于OPLS-AA力场模拟并不一定合理,尽管参数相似。对于OPLS-AA力场的小分子参数生成工具有很多,也包括文献依据。
  • 可能有许多人在使用acpype产生了拓扑文件以后,当gmx grompp时总是报错:[ atomtypes ] Incalid order for directive atomtypes
     在本例中,我没有更改itp中的[ atomtypes ],但也没有报错,这是因为我是直接将三聚氰胺分子放在最前面,而后添加的水,因此顺序是正确的。
     如果你是模拟蛋白-配体之类的体系,通过acpype产生配体的拓扑文件,想直接把itp引用到蛋白top中,那么你就必须把产生的itp文件中的[ atomtypes ]部分移动到蛋白top文件最前面,仅次于#include "XXX/forcefield.itp",即放到这一行的后面即可。