仓库源文，站点原文

layout: post title: GROMACS QM/MM教程4：使用连接原子 categories:

此链接提供了一个很好的QM/MM教程, 但它有点复杂, 而且缺少一些细节以致初学者很难重复整个教程. 因此, 我们利用一个简单的体系来帮助初学者一步一步地进行一个QM/MM模拟.

【陈国俊】 首先在教程前说一点闲话, 我从2016年寒假开始折腾GROMACS的QM/MM, 中间断断续续有小半年时间, 期间也看了用DFTB3进行QM/MM模拟的方法, 但是当时怎么编译安装都没能成功, 不是PLUMED出错, 就是GROMACS出错, 所以最后无奈放弃. 在其他量化软件里面, 除了GAUSSIAN我一个都不熟, 所以就开始折腾GAUSSIAN+GROMACS的QMMM方法. 一开始除了GROMACS官网里的那点简述和李老师及其他热心群友翻译的几篇教程, 可供参考的资料少之又少, 完全就是照葫芦画瓢, 犯了很多低级的错误, 比如新老版本参数不分等. 后来又开始仔细读手册, 在群里向群友们请教. 上次在北京参加李老师的培训班, 又和李老师讨论了这个事情, 最后终于调试成功, 所以整理了一下我的整个运行流程, 希望能给用得到的朋友带来一点帮助.

二点经验很重要:

  1. 手册不要离手, 出错了, 回去翻翻手册会找到很多新的启示;
  2. 要对拓扑文件有充分的了解, 这样才可以正确地编写相应的设置.

0. 准备工作

在开始此教程之前, 我们假定:

【陈国俊】 对于使用GAUSSIAN的情况, 本教程采用的是使用外挂脚本的方法, 所以不需要修改GAUSSIAN的源码. 我用的是李老师整理过的脚本. 使用这个脚本方便的地方在于, 你可以使用任意的QM方法, 而不仅仅限于GROMACS支持的那几种. 比如你可以使用半经验的PM6方法, 如果要考虑弱的色散相互作用, 你可以在自己准备的文件里添加关键词EmpiricalDispersion=GD3.

1. 构建模拟体系, 创建gro和top文件

随意构建一条小肽, 如ALA ALA ALA SER GLY ALA ALA ALA. 我们的目的是把SER和GLY作为QM区域进行QM/MM模拟. 小肽的结构文件命名为pep.pdb.

<figure><script>var Mol1=new ChemDoodle.TransformCanvas3D('Mol-1',642,396);Mol1.specs.shapes_color='#fff';Mol1.specs.backgroundColor='black';Mol1.specs.set3DRepresentation('Ball and Stick');Mol1.specs.projectionPerspective_3D=false;Mol1.specs.compass_display=true; /*//Mol1.specs.atoms_resolution_3D=15; //Mol1.specs.bonds_resolution_3D=15; //Mol1.specs.crystals_unitCellLineWidth=1.5;*/ Mol1.nextFrame=function(delta){var matrix=[];ChemDoodle.lib.mat4.identity(matrix);var change=delta*Math.PI/15000;ChemDoodle.lib.mat4.rotate(matrix,change,[1,0,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,1,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,0,1]);ChemDoodle.lib.mat4.multiply(this.rotationMatrix, matrix)}; Mol1.startAnimation=ChemDoodle._AnimatorCanvas.prototype.startAnimation;Mol1.stopAnimation=ChemDoodle._AnimatorCanvas.prototype.stopAnimation;Mol1.isRunning=ChemDoodle._AnimatorCanvas.prototype.isRunning;Mol1.dblclick=ChemDoodle.RotatorCanvas.prototype.dblclick;Mol1.timeout=5;Mol1.handle=null; var Fmol='78\npep\nN 3.326 1.548 -0.000\nH 3.909 0.724 -0.000\nC 3.970 2.846 -0.000\nH 3.672 3.400 -0.890\nC 3.577 3.654 1.232\nH 3.877 3.116 2.131\nH 4.075 4.623 1.206\nH 2.497 3.801 1.241\nC 5.486 2.705 -0.000\nO 6.009 1.593 -0.000\nN 6.365 3.947 -0.000\nH 5.944 4.447 0.770\nC 7.803 3.771 -0.000\nH 8.039 2.707 -0.000\nC 8.435 4.410 -1.232\nH 8.215 5.478 -1.241\nH 9.514 4.262 -1.206\nH 8.027 3.948 -2.131\nC 8.434 4.409 1.230\nO 7.735 4.980 2.064\nN 9.940 4.342 1.434\nH 10.244 4.675 0.530\nC 10.501 4.962 2.617\nH 9.701 5.405 3.210\nC 11.236 3.938 3.477\nH 12.046 3.493 2.900\nH 11.646 4.431 4.358\nH 10.540 3.159 3.788\nC 11.494 6.055 2.246\nO 11.727 6.310 1.067\nN 12.205 6.840 3.338\nH 12.547 6.071 3.897\nC 13.142 7.872 2.943\nH 13.190 7.923 1.855\nC 12.715 9.235 3.477\nH 11.731 9.488 3.082\nH 12.672 9.202 4.565\nO 13.662 10.225 3.067\nH 13.393 11.084 3.403\nC 14.537 7.579 3.478\nO 14.749 6.573 4.150\nN 15.682 8.537 3.181\nH 15.248 9.405 3.463\nC 16.995 8.223 3.708\nH 16.947 7.290 4.269\nH 17.701 8.116 2.885\nC 17.475 9.333 4.632\nO 16.773 10.320 4.839\nN 18.834 9.225 5.308\nH 19.403 9.013 4.501\nC 19.255 10.303 6.181\nH 18.479 11.068 6.215\nC 19.503 9.797 7.598\nH 20.284 9.037 7.582\nH 19.817 10.627 8.231\nH 18.584 9.365 7.996\nC 20.544 10.940 5.682\nO 21.089 10.528 4.661\nN 21.168 12.097 6.450\nH 21.175 11.694 7.377\nC 22.395 12.676 5.943\nH 22.694 12.156 5.033\nC 22.212 14.155 5.621\nH 21.922 14.691 6.525\nH 23.149 14.564 5.242\nH 21.435 14.270 4.866\nC 23.519 12.557 6.962\nO 23.320 12.023 8.051\nN 24.902 13.102 6.636\nH 24.646 14.031 6.332\nC 25.945 12.970 7.633\nH 25.545 12.468 8.514\nC 27.117 12.153 7.097\nH 27.534 12.648 6.220\nH 27.884 12.070 7.867\nH 26.770 11.157 6.822\nC 26.481 14.332 8.051\nO 26.027 15.360 7.553\n'; Mol1.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol1.startAnimation();Mol1.stopAnimation();function setProj1(yesPers){Mol1.specs.projectionPerspective_3D=yesPers;Mol1.setupScene();Mol1.repaint()}function setModel1(model){Mol1.specs.set3DRepresentation(model);Mol1.setupScene();Mol1.repaint()}function setSpeed1(){Mol1.timeout=500-document.getElementById('spd1').value;Mol1.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol1.startAnimation()}</script><br><span class='meta'>视图: <input type='radio' name='group2' onclick='setProj1(true)'>投影 <input type='radio' name='group2' onclick='setProj1(false)' checked=''>正交    速度: <input type='range' id='spd1' min='1' max='500' onchange='setSpeed1()'/><br>模型: <input type='radio' name='model' onclick='setModel1('Ball and Stick')' checked=''>球棍 <input type='radio' name='model' onclick='setModel1('van der Waals Spheres')'>范德华球 <input type='radio' name='model' onclick='setModel1('Stick')'>棍状 <input type='radio' name='model' onclick='setModel1('Wireframe')'>线框 <input type='radio' name='model' onclick='setModel1('Line')'>线型   <input type='checkbox' onclick='Mol1.specs.atoms_displayLabels_3D=this.checked;Mol1.repaint()'>名称<br>左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动</span><br><figurecaption>Fig.1</figurecaption></figure>

我们将使用AMBER99SB-ildn力场, 为获得拓扑文件, 运行下面的命令并选择AMBER99SB-ildn力场以及TIP3P水模型:

<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">gmx</span> pdb2gmx <span style="color:#666">-f</span> pep.pdb <span style="color:#666">-o</span> pep.gro <span style="color:#666">-p</span> topol.top <span style="color:#666">-ignh</span> </pre></div>

然后添加盒子:

<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">gmx</span> editconf <span style="color:#666">-f</span> pep.gro <span style="color:#666">-o</span> conf.gro <span style="color:#666">-d</span> 1 </pre></div>

这样我们获得了三个文件: 结构文件conf.gro, 拓扑文件topol.top和位置限制文件porse.itp.

2. 添加连接原子(LA, Link Atom)信息, 修改gro和top文件

这里我们不解释为什么要在QM/MM模拟中引入连接原子, 你自己很容易找到答案.

下面的图片展示了我们如何把整个体系分为QM部分和MM部分.

确定了QM区域以后, 我们需要将连接原子添加到结构文件conf.gro, 并在拓扑文件topol.top中定义连接原子的原子类型及其参数, 同时修改与其相关的成键作用参数.

修改结构文件conf.gro

如上图所示, 由于我们划分QM区域与MM区域时破坏了原先的成键, 所以我们要在QM和MM区域的边界处中添加两个连接原子(假定其名称为LA), 因此体系的总原子数应比原来增加2. 打开结构文件conf.gro, 首先将第二行中体系的总原子数由81改为83. 然后我们还要添加这两个LA原子的坐标.

值得注意的是, 如果你的体系中有多个分子, 那么添加LA时就要注意添加的位置. 原则上要添加到LA所属分子的最后一个原子的后面. 由于我们的体系中只有一个分子, 所以直接将LA添加到最后一个原子的后面即可. 但要注意格式与前面的一致.

<table class="highlighttable"><th colspan="2" style="text-align:left">conf.gro</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span>GROningen Mixture of Alchemy <span style="color: #AA22FF; font-weight: bold">and</span> Childrens' Stories <span style="color: #666666">83</span> <span style="color: #666666">1</span>ALA N <span style="color: #666666">1</span> <span style="color: #666666">1.122</span> <span style="color: #666666">1.070</span> <span style="color: #666666">1.170</span> <span style="color: #666666">1</span>ALA H1 <span style="color: #666666">2</span> <span style="color: #666666">1.149</span> <span style="color: #666666">1.019</span> <span style="color: #666666">1.088</span> <span style="color: #666666">1</span>ALA H2 <span style="color: #666666">3</span> <span style="color: #666666">1.149</span> <span style="color: #666666">1.019</span> <span style="color: #666666">1.252</span> <span style="color: #666666">1</span>ALA H3 <span style="color: #666666">4</span> <span style="color: #666666">1.022</span> <span style="color: #666666">1.082</span> <span style="color: #666666">1.170</span> <span style="color: #666666">1</span>ALA CA <span style="color: #666666">5</span> <span style="color: #666666">1.186</span> <span style="color: #666666">1.200</span> <span style="color: #666666">1.170</span> <span style="color: #666666">1</span>ALA HA <span style="color: #666666">6</span> <span style="color: #666666">1.153</span> <span style="color: #666666">1.243</span> <span style="color: #666666">1.086</span> <span style="color: #666666">1</span>ALA CB <span style="color: #666666">7</span> <span style="color: #666666">1.147</span> <span style="color: #666666">1.280</span> <span style="color: #666666">1.293</span> <span style="color: #666666">1</span>ALA HB1 <span style="color: #666666">8</span> <span style="color: #666666">1.192</span> <span style="color: #666666">1.369</span> <span style="color: #666666">1.291</span> <span style="color: #666666">1</span>ALA HB2 <span style="color: #666666">9</span> <span style="color: #666666">1.048</span> <span style="color: #666666">1.294</span> <span style="color: #666666">1.294</span> <span style="color: #666666">1</span>ALA HB3 <span style="color: #666666">10</span> <span style="color: #666666">1.174</span> <span style="color: #666666">1.231</span> <span style="color: #666666">1.376</span> <span style="color: #666666">1</span>ALA C <span style="color: #666666">11</span> <span style="color: #666666">1.338</span> <span style="color: #666666">1.186</span> <span style="color: #666666">1.170</span> <span style="color: #666666">1</span>ALA O <span style="color: #666666">12</span> <span style="color: #666666">1.390</span> <span style="color: #666666">1.074</span> <span style="color: #666666">1.170</span> <span style="color: #666666">....</span>(<span style="#FF0000">略</span>) <span style="color: #666666">8</span>ALA N <span style="color: #666666">71</span> <span style="color: #666666">3.279</span> <span style="color: #666666">2.225</span> <span style="color: #666666">1.834</span> <span style="color: #666666">8</span>ALA H <span style="color: #666666">72</span> <span style="color: #666666">3.298</span> <span style="color: #666666">2.269</span> <span style="color: #666666">1.746</span> <span style="color: #666666">8</span>ALA CA <span style="color: #666666">73</span> <span style="color: #666666">3.384</span> <span style="color: #666666">2.212</span> <span style="color: #666666">1.933</span> <span style="color: #666666">8</span>ALA HA <span style="color: #666666">74</span> <span style="color: #666666">3.338</span> <span style="color: #666666">2.166</span> <span style="color: #666666">2.009</span> <span style="color: #666666">8</span>ALA CB <span style="color: #666666">75</span> <span style="color: #666666">3.501</span> <span style="color: #666666">2.130</span> <span style="color: #666666">1.880</span> <span style="color: #666666">8</span>ALA HB1 <span style="color: #666666">76</span> <span style="color: #666666">3.571</span> <span style="color: #666666">2.123</span> <span style="color: #666666">1.950</span> <span style="color: #666666">8</span>ALA HB2 <span style="color: #666666">77</span> <span style="color: #666666">3.469</span> <span style="color: #666666">2.039</span> <span style="color: #666666">1.854</span> <span style="color: #666666">8</span>ALA HB3 <span style="color: #666666">78</span> <span style="color: #666666">3.539</span> <span style="color: #666666">2.176</span> <span style="color: #666666">1.799</span> <span style="color: #666666">8</span>ALA C <span style="color: #666666">79</span> <span style="color: #666666">3.437</span> <span style="color: #666666">2.348</span> <span style="color: #666666">1.975</span> <span style="color: #666666">8</span>ALA OC1 <span style="color: #666666">80</span> <span style="color: #666666">3.510</span> <span style="color: #666666">2.354</span> <span style="color: #666666">2.044</span> <span style="color: #666666">8</span>ALA OC2 <span style="color: #666666">81</span> <span style="color: #666666">3.392</span> <span style="color: #666666">2.451</span> <span style="color: #666666">1.925</span> <span style="color: #666666">9</span>QMM LA1 <span style="color: #666666">82</span> <span style="color: #666666">0.000</span> <span style="color: #666666">0.000</span> <span style="color: #666666">0.000</span> <span style="color: #666666">9</span>QMM LA2 <span style="color: #666666">83</span> <span style="color: #666666">0.000</span> <span style="color: #666666">0.000</span> <span style="color: #666666">0.000</span> <span style="color: #666666">4.54900</span> <span style="color: #666666">3.43200</span> <span style="color: #666666">3.08000</span> </pre></div> </td></tr></table>

上面添加到两行内容中, 第一列是残基名称, 并不重要, 可以随便用其他名字来替代9QMM; 第二列是原子名称, 我们使用前面指定的名称LA. 这个名称也可以随意写, 只要拓扑文件中指向的原子类型与定义的LA原子类型相同即可. 但与拓扑文件中不一致时grompp会出现警告; 第三列是原子序号, 根据前面的编号继续即可(实际上也可以随意写, 只要总原子数目对即可); 最后三列是原子的坐标, 我们可以全部指定为零, 因为LA属于虚拟位点, 对于虚拟位点的坐标, GROMACS在运行时不会使用结构文件中的值, 而是根据拓扑文件中的信息和结构文件中定义虚拟位点的原子的坐标自动计算. 更具体的说明见后面修改拓扑文件部分.

如果你不喜欢使用零作为LA的坐标, 那么可参考下面的方法来确定LA的坐标(也可参考此链接的adding link atoms部分.

残基3ALA的C31原子和残基4SER的N33原子之间的连接原子应位于这两个原子之间, 它与N33原子之间的距离应正比于C31和N33原子之间的键长. 由于进行QM/MM计算时我们将连接原子视为氢原子, 它距所连氮原子的的距离并不是固定值, 而是根据原结构中的C-N距离来确定. 通常情况下C(=O)-N键长为0.152 nm, N-H键长为0.101 nm, 其比例 $k=0.101/0.152=0.664$. 这样将C31和N33之间的实际键长乘以0.664就可以得到N33和LA1之间的距离. 再根据C31和N33原子的坐标, 就可以计算出LA1的坐标了. 用数学公式表示的话, 可写为 $\vec R\text{LA1}=\vec R\text{N33}+k\vec R_\text{N33-C31}$.

同理, C49和N51之间的LA2原子的坐标也可以计算出来, 只不过LA2是连接在碳原子上的. 通常情况下C(=O)-N键长为0.152 nm, C-H键长为0.108 nm, 其比例为0.108/0.153=0.706.

修改拓扑文件topol.top

首先需要向力场中添加LA的原子类型. 我们既可以将LA类型直接添加到GROMACS自带的力场文件中(或复制一套放于工作目录下, 更稳妥), 也可以将其添加到我们前面得到的topol.top中. 一般建议使用后者, 因为使用前者容易污染GROMACS自带的力场文件. 为此, 我们需要修改前面得到的拓扑文件topol.top. 首先添加LA的原子类型, 将下面的内容添加到#include(力场头文件)语句的下面. 注意, 添加到位置非常重要, 一定要保证[ atomtypes ]位于展开内容之后[ defaults ]的后面.

<table class="highlighttable"><th colspan="2" style="text-align:left">topol.top</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span><span style="color: #AA22FF; font-weight: bold">[ atomtypes ]</span> <span style="#FF0000">LA</span> <span style="#FF0000">1</span> <span style="#FF0000">0</span> <span style="#FF0000">0</span> <span style="#FF0000">A</span> <span style="#FF0000">0</span> <span style="#FF0000">0</span> <span style="color: #008800; font-style: italic">; LA for QM/MM</span> </pre></div> </td></tr></table>

添加了LA的原子类型后, 还应在拓扑文件中分子类型部分添加LA以保持拓扑与结构文件一致. 首先, 我们在拓扑文件[ atoms ]部分的最后增加下列内容:

<table class="highlighttable"><th colspan="2" style="text-align:left">topol.top</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"> <span></span><span style="#FF0000">82</span> <span style="#FF0000">LA</span> <span style="#FF0000">9</span> <span style="#FF0000">QMM</span> <span style="#FF0000">LA1</span> <span style="#FF0000">82</span> <span style="#FF0000">0.000</span> <span style="#FF0000">0</span> <span style="#FF0000">83</span> <span style="#FF0000">LA</span> <span style="#FF0000">10</span> <span style="#FF0000">QMM</span> <span style="#FF0000">LA2</span> <span style="#FF0000">83</span> <span style="#FF0000">0.000</span> <span style="#FF0000">0</span> </pre></div> </td></tr></table>

同样, 第一列是LA原子在所属分子类型中的编号, 第二列对应的力场中的原子类型(注意与原子名称的不同), 第三列是残基编号, 第四列是残基名称(可随意, 不影响结果), 第五列是原子名称(可随意, 但最好保持与结构文件中的一致, 否则grompp会警告. 注意与原子类型的意义不同), 第六列是原子序号, 然后是原子电荷(LA原子不带电, 所以为零), 最后一列是原子的质量(虚拟位点质量为零).

接下来我们需要修改分子的成键相互作用部分, 将QM和MM原子之间的成键类型改为5, 表示两个原子之间是无相互作用的连接. 具体来说, 就是[ bonds ]部分中, 改为31 33 5和49 51 5.

在一些教程中, 还会建议下面的修改, 但根据我的测试, 下面的修改对结果并无影响, 因此可以省略.

最后, 我们需要在拓扑文件中定义LA原子所在虚拟位点. 所谓定义虚拟位点, 就是指定虚拟位点的约束条件. 在之后的模拟过程中, GROMACS会按照定义好的约束条件来更新虚拟位点的坐标. 我们之前说过, 在结构文件中添加LA原子时可以不用给出具体坐标, 简单指定为零即可, 就是因为GROMACS会自动计算.

在[ dihedrals ]后面添加[ virtual_site2 ](对应GROMACS老版本中的[ dummies2 ]. 注意添加的位置):

<table class="highlighttable"><th colspan="2" style="text-align:left">topol.top</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2 3 4</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span><span style="color: #AA22FF; font-weight: bold">[ virtual_sites2 ]</span> <span style="color: #008800; font-style: italic">; Site from funct a</span> <span style="#FF0000">82</span> <span style="#FF0000">33</span> <span style="#FF0000">31</span> <span style="#FF0000">1</span> <span style="#FF0000">0.664</span> <span style="#FF0000">83</span> <span style="#FF0000">49</span> <span style="#FF0000">51</span> <span style="#FF0000">1</span> <span style="#FF0000">0.706</span> </pre></div> </td></tr></table>

其中, site是虚拟位点在分子结构中的编号, from表示根据哪些原子计算虚拟位点. 在我们的示例中是根据31和33来计算LA1, 根据49和51来计算LA2. funct为计算虚拟位点的函数类型, 这里是1, a为约束参数, 对应于我们前面的计算值. 手册 4.7 虚拟相互作用位点对虚拟位点的设置有详细的介绍, 可以参考.

再添加[ constraints ]部分:

<table class="highlighttable"><th colspan="2" style="text-align:left">topol.top</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2 3</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span><span style="color: #AA22FF; font-weight: bold">[ constraints ]</span> <span style="#FF0000">33</span> <span style="#FF0000">31</span> <span style="#FF0000">2</span> <span style="#FF0000">0.152</span> <span style="#FF0000">49</span> <span style="#FF0000">51</span> <span style="#FF0000">2</span> <span style="#FF0000">0.152</span> </pre></div> </td></tr></table>

这里的0.152是C, N间的键长,

到这里GROMACS文件全部修改完成了.

3. 能量最小化

在进行QM/MM模拟之前, 可以先进行一下简单的能量最小化, 用以检查前面的输入文件是否正确, 同时也可以得到虚拟位点的坐标.

<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">gmx</span> grompp <span style="color:#666">-f</span> grompp_em.mdp <span style="color:#666">-c</span> conf.gro <span style="color:#666">-p</span> topol.top <span style="color:#666">-o</span> em.tpr <span style="color:#A2F">gmx</span> mdrun <span style="color:#666">-v</span> <span style="color:#666">-deffnm</span> em </pre></div>

完成后会得到能量最小化之后的结构em.gro. 用可视化软件查看下, LA是否处于正确的位置(如果发现显示的成键连接很乱, 不用管它, 只是显示的问题).

<figure><script>var Mol2=new ChemDoodle.TransformCanvas3D('Mol-2',642,396);Mol2.specs.shapes_color='#fff';Mol2.specs.backgroundColor='black';Mol2.specs.set3DRepresentation('Ball and Stick');Mol2.specs.projectionPerspective_3D=false;Mol2.specs.compass_display=true; /*//Mol2.specs.atoms_resolution_3D=15; //Mol2.specs.bonds_resolution_3D=15; //Mol2.specs.crystals_unitCellLineWidth=1.5;*/ Mol2.nextFrame=function(delta){var matrix=[];ChemDoodle.lib.mat4.identity(matrix);var change=delta*Math.PI/15000;ChemDoodle.lib.mat4.rotate(matrix,change,[1,0,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,1,0]);ChemDoodle.lib.mat4.rotate(matrix,change,[0,0,1]);ChemDoodle.lib.mat4.multiply(this.rotationMatrix, matrix)}; Mol2.startAnimation=ChemDoodle._AnimatorCanvas.prototype.startAnimation;Mol2.stopAnimation=ChemDoodle._AnimatorCanvas.prototype.stopAnimation;Mol2.isRunning=ChemDoodle._AnimatorCanvas.prototype.isRunning;Mol2.dblclick=ChemDoodle.RotatorCanvas.prototype.dblclick;Mol2.timeout=5;Mol2.handle=null; var Fmol='83\nProtein\nN 11.24 10.70 11.70\nH 11.52 10.17 10.89\nH 11.52 10.17 12.52\nH 10.24 10.81 11.71\nC 11.90 12.03 11.71\nH 11.58 12.55 10.81\nC 11.46 12.82 12.94\nH 11.94 13.80 12.93\nH 10.38 12.95 12.93\nH 11.76 12.28 13.84\nC 13.44 11.90 11.68\nO 13.88 10.76 11.68\nN 14.23 12.99 11.67\nH 13.84 13.92 11.72\nC 15.72 12.91 11.72\nH 16.00 11.86 11.71\nC 16.30 13.55 10.45\nH 17.39 13.41 10.46\nH 15.89 13.09 9.56\nH 16.10 14.63 10.45\nC 16.36 13.54 12.98\nO 15.62 14.12 13.77\nN 17.70 13.47 13.17\nH 18.29 13.06 12.47\nC 18.38 14.12 14.31\nH 17.63 14.58 14.94\nC 19.11 13.07 15.16\nH 19.52 13.55 16.05\nH 18.41 12.29 15.48\nH 19.93 12.63 14.60\nC 19.37 15.23 13.90\nO 19.59 15.45 12.71\nN 20.07 16.03 14.99\nH 19.94 15.78 15.96\nC 21.08 17.06 14.66\nH 21.13 17.11 13.57\nC 20.61 18.43 15.17\nH 20.58 18.41 16.26\nH 19.61 18.63 14.78\nO 21.48 19.44 14.71\nH 21.27 20.27 15.15\nC 22.49 16.72 15.18\nO 22.65 15.72 15.87\nN 23.52 17.53 14.87\nH 23.32 18.41 14.40\nC 24.88 17.37 15.41\nH 25.58 17.28 14.58\nH 24.92 16.44 15.98\nC 25.33 18.50 16.34\nO 24.61 19.48 16.52\nN 26.69 18.40 17.03\nH 27.30 17.61 16.84\nC 27.19 19.48 17.89\nH 26.42 20.24 17.93\nC 27.39 18.93 19.31\nH 27.69 19.75 19.96\nH 26.46 18.51 19.68\nH 28.17 18.18 19.30\nC 28.49 20.15 17.37\nO 29.00 19.70 16.35\nN 29.03 21.19 18.04\nH 28.62 21.50 18.91\nC 30.31 21.83 17.64\nH 30.65 21.34 16.73\nC 30.08 23.31 17.31\nH 31.01 23.73 16.93\nH 29.30 23.41 16.56\nH 29.80 23.85 18.22\nC 31.46 21.68 18.67\nO 31.20 21.15 19.75\nN 32.69 22.14 18.36\nH 32.85 22.66 17.51\nC 33.82 22.13 19.31\nH 33.47 21.67 20.23\nC 34.98 21.29 18.77\nH 35.74 21.21 19.54\nH 34.64 20.29 18.51\nH 35.42 21.79 17.90\nC 34.32 23.53 19.73\nO 35.22 23.56 20.60\nO 33.83 24.53 19.15\nCl 19.60 15.50 14.27\nCl 26.29 18.43 16.82\n'; Mol2.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol2.startAnimation();Mol2.stopAnimation();function setProj2(yesPers){Mol2.specs.projectionPerspective_3D=yesPers;Mol2.setupScene();Mol2.repaint()}function setModel2(model){Mol2.specs.set3DRepresentation(model);Mol2.setupScene();Mol2.repaint()}function setSpeed2(){Mol2.timeout=500-document.getElementById('spd2').value;Mol2.loadMolecule(ChemDoodle.readXYZ(Fmol));Mol2.startAnimation()}</script><br><span class='meta'>视图: <input type='radio' name='group2' onclick='setProj2(true)'>投影 <input type='radio' name='group2' onclick='setProj2(false)' checked=''>正交    速度: <input type='range' id='spd2' min='1' max='500' onchange='setSpeed2()'/><br>模型: <input type='radio' name='model' onclick='setModel2('Ball and Stick')' checked=''>球棍 <input type='radio' name='model' onclick='setModel2('van der Waals Spheres')'>范德华球 <input type='radio' name='model' onclick='setModel2('Stick')'>棍状 <input type='radio' name='model' onclick='setModel2('Wireframe')'>线框 <input type='radio' name='model' onclick='setModel2('Line')'>线型   <input type='checkbox' onclick='Mol2.specs.atoms_displayLabels_3D=this.checked;Mol2.repaint()'>名称<br>左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动</span><br><figurecaption>Fig.2</figurecaption></figure>

4. 创建索引文件

在运行QM/MM模拟时, GROMACS需要知道哪些原子应使用QM进行处理. 所以就需要我们指定处于QM区域内的原子, 并将其单独列为一组, 保存在索引文件中. 对于小体系, 可以直接手写索引文件. 对复杂点的体系使用gmx make_ndx可以轻松完成. 要注意的是, 在生成QM原子组时, 要把虚拟原子也包括进来, 否则的话, GAUSSIAN运行时会缺少连接原子, 构型不对, 得到错误的结果.

要创建QM原子组, 对我们的体系, 可使用下列命令:

<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">gmx</span> make_ndx <span style="color:#666">-f</span> em.gro <span style="color:#666">-o</span> index.ndx <span style="color:#A2F">></span> r 4 5 | r QMM <span style="color:#A2F">></span> name 14 QMatoms <span style="color:#A2F">></span> q </pre></div>

这样我们就得到了一个索引文件index.ndx. 其中包含了QM原子组QMatoms:

<table class="highlighttable"><th colspan="2" style="text-align:left">index.ndx</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2 3</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span><span style="color: #AA22FF; font-weight: bold">[ QMatoms ]</span> <span style="#FF0000">33</span> <span style="#FF0000">34</span> <span style="#FF0000">35</span> <span style="#FF0000">36</span> <span style="#FF0000">37</span> <span style="#FF0000">38</span> <span style="#FF0000">39</span> <span style="#FF0000">40</span> <span style="#FF0000">41</span> <span style="#FF0000">42</span> <span style="#FF0000">43</span> <span style="#FF0000">44</span> <span style="#FF0000">45</span> <span style="#FF0000">46</span> <span style="#FF0000">47</span> <span style="#FF0000">48</span> <span style="#FF0000">49</span> <span style="#FF0000">50</span> <span style="#FF0000">82</span> <span style="#FF0000">83</span> </pre></div> </td></tr></table>

5. 创建运行参数文件

由于示例体系很简单, 也不考虑溶剂, 所以我们省略预平衡阶段, 直接进行MD模拟. 进行MD前, 需要首先设置好GAUSSIAN运行所需的环境变量. 然后在grompp.mdp文件中指定QM/MM选项, 选项中的参数取决于你的体系. grompp.mdp文件中的其他设置与通常模拟没有区别.

我们指定的QM/MM选项如下:

<table class="highlighttable"><th colspan="2" style="text-align:left">grompp.mdp</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2 3 4 5 6 7 8</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span>QMMM <span style="color: #666666">=</span> yes<span style="color: #008800; font-style: italic"> ; QM/MM模拟开关打开</span> QMMM<span style="color: #666666">-</span>grps <span style="color: #666666">=</span> QMatoms<span style="color: #008800; font-style: italic"> ; QM原子所属组</span> QMmethod <span style="color: #666666">=</span> B3LYP<span style="color: #008800; font-style: italic"> ; QM的理论方法</span> QMbasis <span style="color: #666666">=</span> <span style="color: #666666">6-31</span>G<span style="color: #666666">*</span><span style="color: #008800; font-style: italic"> ; QM使用的基组</span> QMMMscheme <span style="color: #666666">=</span> normal<span style="color: #008800; font-style: italic"> ; QM/MM的类型(只有一个QM组用normal即可)</span> QMcharge <span style="color: #666666">=</span> <span style="color: #666666">0</span><span style="color: #008800; font-style: italic"> ; QM原子的净电荷</span> QMmult <span style="color: #666666">=</span> <span style="color: #666666">1</span><span style="color: #008800; font-style: italic"> ; QM原子的自旋多重度</span> SAsteps <span style="color: #666666">=</span> <span style="color: #666666">0</span><span style="color: #008800; font-style: italic"> ; GMX2016.4需要此选项</span> </pre></div> </td></tr></table>

这里需要额外说明的是, GROMACS支持的QMmethod只有几种, 如果换用其他基组grompp时还会报错, 所以我们随意指定GROMACS支持的选项即可. 我们在这里定义的QMmethod和QMbasis实际上并不会用到, 实际运行时使用的设置来自我们提前准备好的GAUSSIAN模板文件_Gkwd.gjf. 这样我们才可以使用GAUSSIAN支持的任意理论方法和基组. 根据GROMACS QM/MM教程2：编译设置及简单运行的说明准备_Gkwd.gjf. 可以参考如下示例:

<table class="highlighttable"><th colspan="2" style="text-align:left">_Gkwd.gjf</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2 3 4</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span>%chk<span style="color: #666666">=</span>input <span style="color: #008800; font-style: italic">#T HF NoSymm</span> Force <span style="color: #B8860B">Punch</span><span style="color: #666666">=</span>Derivatives Charge <span style="color: #B8860B">Prop</span><span style="color: #666666">=(</span>Field, Read<span style="color: #666666">)</span> </pre></div> </td></tr></table>

6. 运行模拟

到目前为止, 我们已经准备好了运行模拟所需的所有文件. 输入下面的命令进行模拟:

<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">gmx</span> grompp <span style="color:#666">-f</span> grompp.mdp <span style="color:#666">-c</span> em.gro <span style="color:#666">-p</span> topol.top <span style="color:#666">-n</span> index.ndx <span style="color:#666">-o</span> nvt.tpr <span style="color:#666">-maxwarn</span> 3 <span style="color:#A2F">gmx</span> mdrun <span style="color:#666">-deffnm</span> nvt </pre></div>

如果不出意外, 你应该看到类似下面的输出:

并且轨迹文件和日志文件也在更新. 这说明模拟已经成功地运行了. 祝贺你. 不过也别高兴的太早, 因为QM/MM运行过程中经常会出现构型问题, 以及GAUSSIAN Link 502自洽迭代不收敛的错误, 你可能还需要修改mdp选项, 并增加一些GAUSSIAN控制SCF收敛的选项才能顺利地模拟下去. 由于本教程仅为了示例如何运行一个QM/MM模拟, 所以我们就不再继续讨论这些细节了. 此外, 要进行一次有意义的模拟, 你需要学习更多分子模拟的知识, 这也超出了此教程的范围.

7. 结语

【陈国俊】 在结束部分, 还有点闲话想说一下. 根据我的测试, 这种QM/MM方法计算速度极慢, 只适用于小体系的短时间模拟. 即便对于本教程使用的这样小的示例体系, 在我的电脑上(i7-6700k), 使用PM6D3这种半经验方法模拟100 ps都需要两三个小时, 而使用B3LYP/6-31G*这种并非很高精度的方法时测试的速度是44 day/ns. 我估计运行如此之慢的原因一是因为GROMACS不能并行, 二是GAUSSIAN本身运算的速度就很慢, 三是这两者中间还需要一个脚本进行数据交换, 其速度可想而知. 如果想做大体系长时间的模拟, 还是趁早换其他的门路, 不要在这折腾浪费时间了.