分子模拟/表面活性剂Gromacs分子动力学模拟实践

来自KlniuWiki
跳转至: 导航搜索
分子模拟/表面活性剂Gromacs分子动力学模拟实践
学科类别 计算化学
来源属性 原创

1 前言

研究生三年,做了两年左右的分子动力学模拟,挺遗憾的是,虽然在我看来十分用心了,但结果却不是那么别人满意,毕业了,匆匆写了一篇论文,仅仅给自己一个交代,对老师三年的培养却无报答,虽然遗憾也很自责,却只能心里想想,现在确实没有心力和时间做这么大部头的东西了。

就两年的感觉来说,丁老师曾经评价我是一个完美主义者,不知是褒还是无奈的贬,工作了这么久,发现活在这个世界认真必定会输的,还好有一个信仰支撑着,会得到很多东西,但这些东西没有多少人会认可的,而且也无法与别人道,因此,只能埋在心里,自得其乐。

在毕业之前把分子模拟的一套东西想交给师妹,但结果也不甚理想,一方面是分子模拟的力场构建确实挺复杂;另外一方面是自己做的这一套东西,均是适合于自己的思维习惯,没有特别考虑别人使用时候的感受,因此,前几天师兄问下做分子模拟的经验,于是就想把这些经验赶紧整理出来,不然再过一段时间,真的就全部忘记了。公开此篇文章,一方面纪念这三年,另一方面希望能帮助一些人,我的这些经验很多都来源于论坛、邮件列表等等,实在是别人思想的加工,因此继续回报下去。

有一点需要重要说明的是:力场构建的基础是利用ATB生成的分子及力场合并起来的,ATB有一篇文章是说明其力场构建的过程的,但详细的参数如何选择、程序如何运行,我也不知道,当时我只取用了其优化分子结构的一部分经验,至于怎么去匹配现有的GROMOS53a6力场,却并不知晓,我只知道它的构建过程是经过论证的,因此就拿来使用了。在其FAQ中,有下面一段话:

The ATB only generates a complete topology for molecules containing <40 atoms. How can I generate a reliable topology for a larger molecule?

Although the ATB only generates a complete topology for molecules containing <40 atoms, it can be used to create an initial template for a larger ligand molecule. This template contains all atom-type information and bonded parameters. In some cases, initial estimates of the charges may also be given. Before use, these templates need to be examined and where necessary edited manually. In order to obtain more precise parameters, a larger molecule can be broken into smaller fragments and the parameters from these fragments can be used to complete the template of the larger molecule.

意思是说,对于大于40个原子的分子,可以拆开分子至小的片段,然后用小的分子的力场信息合并生成大的力场,原理即是如此,后面会详细说明。如果有人想研究完整的生成力场的方法,可以努力下去~~

说明:

  1. 这里的文章适用于表面活性剂分子,对于其他分子是否适用,主要是考察力场的适用性及模拟方法的合理性。
  2. 力场使用GROMOS53a6,分子力场的构建及模拟方法的选择与力场均有关系,请注意。
  3. 本文以表面活性剂分子4-十六烷基苯磺酸盐(4-C16)为例来说明整个过程。
  4. 有一些程序和脚本可以简化本文的一些操作,源代码和使用说明参见:分子模拟/Gromacs中53a6力场的构建

2 表面活性剂分子结构优化及力场构建

2.1 依赖的程序或网站

  • gromacs, 分子动力学模拟软件。
  • ATB,在线生成GROMOS53A6 拓扑文件。
  • GAMESS(US),量子力学软件,可以用来优化分子结构。安装方法参见:GAMESS(US),安装说明最后修订日期为2013年4月,于软件新版本不一定适用,仅参考。
  • Avogadro,可视化分子文件编辑软件,支持多种常见分子格式。
  • vmd, 可视化软件,支持pdb, gro, trj, xtc等文件
  • pymol, 与vmd类似,支持gro, pdb等
  • packmol, 分子堆砌软件,可以堆砌多种结构的文件。

上述软件,除了GAMESS(US)之外,都可以在archlinux官方或aur源中找到。

关于Gamess(US)的安装,可以参考如下安装说明:GAMESS(US)

2.2 原理

  1. 制造一个初始分子文件, pdb;
  2. 对此分子文件运用量子力学方法进行结构优化,获得一个优化过的分子文件(pdb)和对应各个原子的电荷;
  3. 将分子拆为多个片段的pdb;
  4. 将上述片段分别在ATB上生成拓扑文件;
  5. 利用我们自己优化过的pdb、电荷和ATB生成的片段拓扑文件,生成一个分子完整的拓扑文件。

2.2.1 制造初始pdb

以4-十六烷基苯磺酸盐为例,其分子结构如下图,此处不考虑其携带的阳离子,离子是电离在水中的,而且在实际模拟过程中是通过genion添加的。

protected.png
4-C16分子式
2.2.1.1 绘制分子

PDB的制作可以使用分子软件来完成,Avogadro是一个比较厉害的分子编辑软件,3D显示。

绘制方法:

  • 选中工具栏上的“绘制工具(F8)”图标,然后左边会显示下一个原子的元素类型。
  • 在右边空白处拖动即可绘制。
  • 右边还有一些工具,比如选择等等,试试就知道是什么功能了。

另外一种方法是使用SMILES,可以写一个字符串来表示一个分子,然后导入分子编辑软件即可生成PDB,比如4-C16的SMILES:CCCCCCCCCCCCC(c1ccc(S(=O)(=O)O)cc1)CCC。

也有一些在线转换工具,可以将SMILES直接生成PDB,比如http://cactus.nci.nih.gov/translate/

需要注意的是:

  • 生成的分子默认是携带H原子的,可以选中要去除氢原子所连接的原子(可以多选),使用Avogadro菜单中构建->删除氢原子来删除,也可以单独选中氢原子删除。
  • Ctrl+A全选所有原子,再选择删除氢原子就可删除全部。
  • 对于一些氢原子,即使使用联合原子力场也是不能删除的,比如苯环上的氢,氢氧根键上的氢。
  • 当删除了不该删除的氢原子后,再通过构建->添加氢原子来添加。
  • 本文所使用的力场保留了所有的氢原子,但4-C16在生成分子时磺酸基上的氢原子会自动添加,需要删除。
2.2.1.2 优化结构

手动绘制的分子,其结构与分子真实结构相差很远,可以使用Avogadro来初步优化分子,这样可以极大缩短使用GAMESS(US)优化结构的时间。

方法:

  • Avogadro菜单栏:扩展->分子机理->建立力场,步骤数9999,其他默认即可,确定。
  • 扩展->优化几何结构,快捷键Ctrl+Alt+O。如果一次优化后的结构看起来仍不合理,可多次优化。
2.2.1.3 导出pdb

使用Avogadro另存,类型选择pdb。

2.3 分子结构优化与原子电荷计算

安装完成GAMESS(US)之后,GAMESS(US)需要配置文件inp即可完成分子结构的优化和电荷的计算,inp文件中包含了优化的参数及分子的信息。

2.3.1 GAMESS(US) inp文件的生成

可用Avogadro生成inp文件,用Avogadro打开pdb文件,打开“菜单栏->扩展->gamess->输入生成器",即可看到配置文件的内容,以4C-16为例:其内容如下:

!   File created by the GAMESS Input Deck Generator Plugin for Avogadro
 $BASIS GBASIS=STO NGAUSS=3 $END
 $CONTRL SCFTYP=ROHF RUNTYP=ENERGY MULT=2 $END

 $DATA 
Title
C1
C     6.0   -12.77100    -6.17700     0.04300
C     6.0   -11.32700    -6.00700    -0.00900
C     6.0   -10.57700    -4.86300    -0.36900
C     6.0    -9.22400    -4.58400    -0.45700
C     6.0    -8.54000    -3.44300    -0.81500
C     6.0    -7.18600    -3.21800    -0.88600
C     6.0    -6.51300    -2.07600    -1.24400
C     6.0    -5.16000    -1.85000    -1.31500
C     6.0    -4.49300    -0.70700    -1.67200
C     6.0    -3.13100    -0.48200    -1.74300
C     6.0    -2.41600     0.61300    -2.13800
C     6.0    -1.04500     0.88700    -2.27400
C     6.0    -0.44700     2.09100    -2.36200
C     6.0    -1.02800     3.15400    -1.60600
C     6.0    -1.86900     2.92200    -0.50600
C     6.0    -2.36800     3.93100     0.28000
C     6.0    -1.86100     5.23900     0.11000
S    16.0    -2.49700     6.65300     1.02200
O     8.0    -3.97300     6.67800     0.61000
O     8.0    -2.26100     6.24000     2.48500
O     8.0    -1.77500     7.79700     0.61400
C     6.0    -0.86700     5.45600    -0.84200
C     6.0    -0.54000     4.46500    -1.74500
C     6.0     0.86900     2.22900    -2.67200
C     6.0     1.70400     3.35100    -3.00800
C     6.0     3.14300     3.60400    -3.03000
H     1.0   -13.00200    -7.17400     0.35700
H     1.0   -13.18600    -6.00400    -0.92800
H     1.0   -13.18800    -5.47900     0.73800
H     1.0   -10.99400    -6.76600    -0.68600
H     1.0   -11.14000    -5.98400     1.04400
H     1.0   -10.93200    -4.11900     0.31300
H     1.0   -10.72400    -4.97600    -1.42300
H     1.0    -8.85900    -5.32100    -1.14100
H     1.0    -9.04800    -4.50500     0.59600
H     1.0    -8.89400    -2.69800    -0.13300
H     1.0    -8.73500    -3.48100    -1.86600
H     1.0    -6.83000    -3.96100    -1.56800
H     1.0    -6.98600    -3.18500     0.16500
H     1.0    -6.86900    -1.33300    -0.56200
H     1.0    -6.71300    -2.10900    -2.29500
H     1.0    -4.80200    -2.59200    -1.99800
H     1.0    -4.95600    -1.82200    -0.26500
H     1.0    -4.85000     0.03400    -0.98900
H     1.0    -4.69800    -0.73400    -2.72200
H     1.0    -2.77800    -1.26200    -2.38500
H     1.0    -2.96600    -0.36900    -0.69200
H     1.0    -2.74900     1.37400    -1.46400
H     1.0    -2.59600     0.47900    -3.18500
H     1.0    -0.75200     0.37300    -3.16500
H     1.0    -0.72800     0.62000    -1.28800
H     1.0    -0.77100     2.22200    -3.37300
H     1.0    -2.12400     1.94800    -0.27900
H     1.0    -3.09900     3.73700     0.98200
H     1.0    -0.37300     6.36100    -0.87200
H     1.0     0.07600     4.69200    -2.54100
H     1.0     0.98800     1.57900    -3.51400
H     1.0     1.22100     2.13800    -1.66600
H     1.0     1.34000     4.11800    -2.35700
H     1.0     1.62100     3.21100    -4.06600
H     1.0     3.48900     3.78400    -2.03300
H     1.0     3.34500     4.46100    -3.63700
H     1.0     3.64900     2.75300    -3.43500
 $END

$DATA以上的为优化参数,以下的为原子信息。

你也可以使用Avogadro上面的选项选择来生成想要的参数,完成后导出inp文件,为方便说明,此处使用文件名4c-16.inp。

或者使用以下参数来替换inp文件中$DATA以上的内容:

 $CONTRL   RUNTYP=OPTIMIZE SCFTYP=RHF DFTTYP=B3LYP
            COORD=CART ICHARG=-1                  $END
 $SCF       DIRSCF=.T. CONV=1.0E-08                   $END
 $SYSTEM   TIMLIM=50000 MWORDS=128               $END
 $BASIS     GBASIS=N31 NGAUSS=6 NDFUNC=1          $END
 $STATPT    NSTEP=200  OPTTOL=1.0E-06               $END
 ! Calculate electric potential
 $ELPOT     IEPOT=1 WHERE=PDC                      $END
 $PDC       PTSEL=CONNOLLY CONSTR=CHARGE       $END

 $DATA

上述参数表示:优化时使用B3LYP/6-31+g(d)基组,并运用Kollman-Singh方法计算电荷;定义电荷为分子所带电荷数,此处为-1;使用笛卡尔坐标类型;SCF的收敛标准设定至10-8,几何优化精度至10-6 Hartree/Bohr。

Warning: 记得根据你的分子类型修改ICHARGE电荷数量。

参数的选择参考了ATB所使用的方法,参见: MALDE A K, ZUO L, BREEZE M, et al. An Automated force field Topology Builder (ATB) and repository: version 1.0[J]. Journal of Chemical Theory and Computation, American Chemical Society, 2011, 7(12): 4026–4037.

2.3.2 结构优化及电荷计算

完成后运行rungms即可:

rungms 4c-16.inp 01 6 >& 4-c16-gam.log

此处,4c-16为输入文件,01为gamess的版本号,6为线程数,最后把结果导出到4-c16-gam.log

完成后使用Avogadro打开log文件,导出为pdb文件,如果要使用联合原子力场,需要删除其中不必要的氢原子。

为简化上述过程,有一个程序可以直接通过pdb文件来完成优化,其接受一个pdb文件和一些选项,生成gamess的结果log文件,和优化后的pdb文件。详细使用方法参见:分子模拟/Gromacs中53a6力场的构建#分子优化

2.3.3 提取电荷

2.4 电荷优化

2.5 完整分子拓扑文件模板生成

2.6 分子分割

2.7 分子片段拓扑文件获取

2.8 拓扑文件合并

2.8.1 电荷

2.8.2 pairs

2.8.3 angles

2.8.4 dihedrals

3 分子模拟流程

3.1 模型构建

3.1.1 Gromacs方法

3.1.2 Packmol

3.2 模拟命令

4 结果分析

protected.png
4C-16分子结构