Gromacs基于OPLS-AA力场的聚合物建模及模拟

本文内容包括：

1. OPLS-AA力场简介/高分子模拟常用力场；

2. 聚合物建模方法；

3. Gromacs添加非标准残基建模聚合物的方法；

4. 聚合物PDB文件生成；

5. 开始模拟。

第一步：

尽管Gromacs与Amber都包含适合蛋白质与核酸体系模拟的力场，但到目前为止罕有专为人工聚合物体系开发的力场。不同于蛋白质/核算具有有限种类的残基，聚合物原则上可以有无限种不同的单体，去为每一种单体都设定优化参数是不可能的。

因此在实际作业中，往往使用普通的全原子力场对高分子体系进行模拟，包括COMPASS（仅MS）、CVFF（仅MS），OPLS-AA, DREIDING，GAFF，MM3，AMBER等。其中OPLS-AA是作为高分子模拟中比较常用的一种，全称Optimized Potentials for Liquid Simulations，转为准确再现凝聚态物化性质开发。被Gromacs原生支持，并且也有成熟的拓扑文件生成工具（LigParGen，TPPmktop，GMXtop），使用起来较为便捷。

第二步：

高分子模拟的一大难点在于如何构建具有较高聚合度的分子以及拓扑文件。

对于聚合度低的情况（10-40个重复单元），最简单的方法是直接在GaussView中构建单体判断，重复拼接，保存为pdb格式然后将得到的分子直接丢到TPPmktop得到拓扑文件。

当聚合度较高时，以上方法就不再适用，一方面手动拼接上百个单元过于耗时，另一方面过大的大分子也不被TPPmktop等拓扑文件生成程序支持。此时推荐的方法是添加非标准残基来生成拓扑文件。残基的概念是来自于蛋白质，gromacs处理蛋白质建模的方法是通过pdb2gmx插件读取蛋白质pdb文件记载的残基名，匹配数据库对于蛋白质每个残基成键参数的定义来生成拓扑文件。对于人工高分子也可以用这种方法，通过自定义片段的信息（RTP文件），也就是添加非标准残基，来用pdb2gmx自动生成高分子的拓扑文件。

第三步：

本文所有的文本编辑（RTP，PDB文件等）都默认使用Notepad++，会使用到一些特殊功能，建议安装。

首先在GaussView中构建单体单元，本例中使用的是DMA单元（N,N-DiMethylAcrylamide）

全选原子，点击GaussView窗口菜单中的Custom Fragments（在DNA图标右边），选中一个Group右键点击New Fragment from active molecule，就得到了自定义的基团，然后构建一个三单元的聚合物用于获取头基、尾基及中间单元的RTP文件，保存为PDB格式。这里建议修改原子名称，将同元素原子用编号区别开，例如按照如下格式，记得要保持原有对齐，比原来多了几个字符就在其后删去多少空格（参考压缩文件1.pdb）：

C1 H11 H21 H31 C2 H21 C3 O N C4 H41 H42 H43 C5 H51 H52 H53……

（第二个重复单元也可以C1开始编号）

将所得pdb文件提交到TPPmktop ，点击Load PDB file for making all-atom OPLS topology右边的按钮提交文件，得到生成的itp与rtp文件。我们只需要rtp文件。

rtp文件包含了指认的原子类型、成键类型、improper二面角，其它未指认的均使用opls-aa的默认参数。

[ atoms ]从左到右包括了原子名称，指认的原子种类，电荷，所在碳原子编号

C1 opls_135 -0.180 1

[ bonds ]只包含了原子连接的信息，参数依照指认的原子类型读取自opls-aa参数库

C5 H43

[ impropers ]包含了涉及的原子与improper二面角类型，参数读取自opls-aa参数库

C2 O C3 H53 improper_O_C_X_Y

在这里我们需要依据这些数据构建自己的rtp文件，由于头基包含17个原子，从得到的rtp文件顺序复制前17个原子信息即可，注意核对电荷总和，对于中性的残基片段应为0；但bonds的顺序是乱的，这里我们就不使用得到的rtp文件，自己按照成键规则编辑，注意一项是，与下一个片段相连的键，需要在下个片段的原子前加上+，表示是下个片段的头原子，在下个片段时则在上个片段相连的原子前加上-号，表示与上个片段的末原子相连；improper项根据这个片段的原子从得到的rtp中复制，如果出现与下个原子相连的情况，也需要加上+号。

最终效果（这里只列举了其中一个残基，完整见压缩文件DMA.rtp）：

[ BEG ] ；残基名称，自定义

[ atoms ]

C1 opls_135 -0.180 1

H11 opls_140 0.060 1

H12 opls_140 0.060 1

H13 opls_140 0.060 1

C2 opls_137 -0.060 2

H21 opls_140 0.060 2

C3 opls_235 0.500 3

N opls_239 -0.140 4

O opls_236 -0.500 5

C4 opls_243 -0.110 6

C5 opls_243 -0.110 7

H41 opls_140 0.060 6

H42 opls_140 0.060 6

H43 opls_140 0.060 6

H51 opls_140 0.060 7

H52 opls_140 0.060 7

H53 opls_140 0.060 7

[ bonds ]

C1H11

C1H12

C1H13

C1C2

C2H21

C2+C1 ；与下个片段的头原子相连

C2C3

C3O

C3N

NC4

NC5

C4H41

C4H42

C4H43

C5H51

C5H52

C5H53

[ impropers ]

C2 N C3 O improper_O_C_X_Y

C3 C4 N C5 improper_Z_N_X_Y

将编辑好的rtp文件放到 你的gromcas文件夹/share/gromacs/top/oplsaa.ff内即可使用

第四步：

对于构建较高分子量的聚合物，通过脚本复制重复单元来生成高分子链是可行的方法。编辑之前得到的1.pdb，删去头基尾基除了碳的所有原子（这个在notepad中删除更快点），然后修改头尾碳为氯原子（用于标识），并用GaussView修改C-Cl键长为原来的一半，然后另存为gjf文件（如果gjf文件带有残基名，那么先存为mol2再存为gjf除去残基名），用notepad打开gjf文件，只复制原子坐标信息，粘贴到新文件中，在首行插入片段的原子数，再留一空行，保存为1.xyz，然后使用以上链接的脚本（或者压缩包的polymerize.sh）生成聚合物xyz文件。

直接生成的聚合物存在几何构型不合理的问题，往往会被VMD判定存在不合理的成键，但在能量最小化中可以解决。

要使用pdb2gmx生成拓扑文件，我们需要得到带有残基号的PDB文件。所以下一步是用得到的xyz文件得到带有残基名的pdb文件。

pdb文件有一点麻烦的是文件里面的每个记录都有着严格的格式。每个记录中的字段, 如标识, 原子名称, 原子序号, 残基名称, 残基序号等, 不仅要按照严格的顺序书写, 而且每个字段所占的字符串长度, 及其所处的位置都是严格规定好的。在这里我使用自己编辑的excel文档xyz2pdb.xlsx与notepad来转换xyz格式并确保得到的pdb文件格式合法

如图所示，在P列左端的为工作区，用于记载原始数据，黑粗闸线框住的是输出区

具体的使用方法如下：

首先另开表格导入高分子的xyz文件（包含原子名称，坐标），将原子名称一列复制到O列，坐标复制到KLM三列。C列是原子序号。E列的原子名参考RTP的命名方式。G列用于修改残基名，自己定义。I列是判断残基序号的，你需要根据自己体系的残基原子数来修改I列，这里我在MID残基处使用=IF(I18<10," “,(IF(I18<100,” “,” ")))函数来判断残基序号。

其它空列（BDFHJ）并非无用，包含了为了对齐pdb格式所用的空格及函数，根据自己需要修改（一般不需要动）。

修改完后，直接复制整个黑粗闸线区的内容到notepad上并转换tab为空格（在编辑>空白字符操作一栏，这步很重要）即可得到标准pdb格式的文件。

第五步：

做好以上工作后，就可以使用pdb2gmx指令得到top拓扑文件与gro坐标文件；不过由于指认二面角的函数类型有误，所得到的拓扑文件还需要经过修改。

打开拓扑文件，修改[ dihedrals ]项的函数类型为“3”（第五个个数字），improper[ dihedrals ]项的函数类型为“4”，这个操作在notepad上批量修改较为方便（例如替换“1 i”为“4 i”，含空格）。

gro文件可能也需要修改，当聚合度很高时，某个坐标数值也会变得很大，自动创建的模拟盒子大小也会不足够，这时修改gro文件最后一行为一个合适的数值（三个数字均为最长值的两倍），进行能量最小化与真空NVT模拟后，笔直的高分子就会蜷缩为一团，这时候再调整模拟盒子大小为你期望的值，然后用solvatebox 指令添加溶剂。