layout: post title: 使用GROMACS计算分子间相互作用 categories:
GMX将分子间相互作用分为范德华相互作用和库伦相互作用. 范德华相互作用是短程相互作用, 所以大部分能量都以 EVdw-SR 体现, 而长程部分以色散校正的形式体现, 即 EDisper.-corr.. 库伦相互作用即静电相互作用, 是长程相互作用. GMX处理静电时可以采用两种不同的方式, 截断和PME. 前者用于孤立的团簇体系, 后者用于周期性体系. 如果使用截断方式, 库伦相互作用就是简单的对势加和, 分解也很简单. 如果采用PME方式, 截断距离内的短程库伦 ECoulomb-SR 和倒易空间中的长程库伦 ECoul.-recip. 加起来才是整个体系的库伦相互作用. 这种方法也是能够写成对势累加形式的, 但过于麻烦, 所以一般程序不会支持, 因此也就导致没法直接将整个库伦作用分解到分子之间的贡献. 所以, 使用GMX计算两种或者两种以上分子的分子间相互作用时, 利用能量组功能, 可以直接得到分子间的范德华相互作用, 但对于分子间的库伦作用, 如果你使用了PME的话, 就只能通过相互作用的定义进行计算.
理论是这样, 但实际大多数时候, 我们需要计算的, 或者说我们感兴趣的都是一些孤立分子之间的库伦相互作用, 而不是这些分子处于周期性盒子中, 具有一定浓度情况下的库伦相互作用. 所以你在要进行PME的库伦相互作用能分解之前, 先想一想, 这是不是你需要的. 如果不是, 那你使用简单截断方法, 将截断值设为无穷大, 就可以得到结果了.
如果你铁了心要进行PME库伦相互作用的分解, 那就继续向下看.
体系中有A, B两种分子, 离子和水, 要计算A, B两种分子之间的相互作用
gmx trjconv将A, B的轨迹提取出来rerun 方法, 重跑一遍轨迹, 计算AB体系中的相互作用.首先需要AB体系的 AB.tpr文件, 所以先抽取一帧只含AB分子的构型
编辑拓扑文件, 只留下A和B, 保存为 AB.top, 然后生成tpr文件, 并重跑
此外, 要获得只含A和B的tpr文件, 使用gmx convert-tpr -n可能更方便.
gmx energy 分析计算相互作用的能量具体计算方法如下表所示：
<table id='tab-0'><caption> <input type='button' id='tab-0_tog' value='折叠表格' onclick="togtab('tab-0', this.value)"></caption><tr> <th rowspan="1" colspan="1" style="text-align:center;">能量项</th> <th rowspan="1" colspan="1" style="text-align:center;">AB</th> <th rowspan="1" colspan="1" style="text-align:center;">A</th> <th rowspan="1" colspan="1" style="text-align:center;">B</th> <th rowspan="1" colspan="1" style="text-align:center;">A-B</th> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">EVdw-SR</td> <td rowspan="1" colspan="1" style="text-align:center;">V1</td> <td rowspan="1" colspan="1" style="text-align:center;">V2</td> <td rowspan="1" colspan="1" style="text-align:center;">V3</td> <td rowspan="1" colspan="1" style="text-align:center;">V1-V2-V3</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">EDisper.-corr.</td> <td rowspan="1" colspan="1" style="text-align:center;">D1</td> <td rowspan="1" colspan="1" style="text-align:center;">D2</td> <td rowspan="1" colspan="1" style="text-align:center;">D3</td> <td rowspan="1" colspan="1" style="text-align:center;">D1-D2-D3</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">ECoulomb-SR</td> <td rowspan="1" colspan="1" style="text-align:center;">C1</td> <td rowspan="1" colspan="1" style="text-align:center;">C2</td> <td rowspan="1" colspan="1" style="text-align:center;">C3</td> <td rowspan="1" colspan="1" style="text-align:center;">C1-C2-C3</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">ECoul.-recip.</td> <td rowspan="1" colspan="1" style="text-align:center;">R1</td> <td rowspan="1" colspan="1" style="text-align:center;">R2</td> <td rowspan="1" colspan="1" style="text-align:center;">R3</td> <td rowspan="1" colspan="1" style="text-align:center;">R1-R2-R3</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">Etotal</td> <td rowspan="1" colspan="1" style="text-align:center;">V1+D1+C1+R1</td> <td rowspan="1" colspan="1" style="text-align:center;">V2+D2+C2+R2</td> <td rowspan="1" colspan="1" style="text-align:center;">V3+D3+C3+R3</td> <td rowspan="1" colspan="1" style="text-align:center;">(V1+D1+C1+R1)-(V2+D2+C2+R2)-(V3+D3+C3+R3)</td> </tr> <tfoot><tr><td colspan="5" style="text-align:left"> 注：AB为A和B整个体系的相互作用；A-B为A和B之间的相互作用.<br> </td></tr></tfoot> </table>下面是计算的两种染料分子(各27个)之间的相互作用, 在水溶液中, 染料分子均呈阴离子形式存在(每个染料分子均带两个负电荷).
<table id='tab-1'><caption>溶液中, 染料分子间的平均非键相互作用能(100 ns至200 ns) <input type='button' id='tab-1_tog' value='折叠表格' onclick="togtab('tab-1', this.value)"></caption><tr> <th rowspan="1" colspan="2" style="text-align:center;">分子间能量项(kJ/mol)</th> <th rowspan="1" colspan="1" style="text-align:center;">E(AB)</th> <th rowspan="1" colspan="1" style="text-align:center;">E(A)</th> <th rowspan="1" colspan="1" style="text-align:center;">E(B)</th> <th rowspan="1" colspan="2" style="text-align:center;">分子间相互作用能ΔE(A-B)</th> </tr> <tr> <td rowspan="2" colspan="1" style="text-align:center;">vdw</td> <td rowspan="1" colspan="1" style="text-align:center;">LJ(SR)</td> <td rowspan="1" colspan="1" style="text-align:center;">-2362.31</td> <td rowspan="1" colspan="1" style="text-align:center;">-1044.44</td> <td rowspan="1" colspan="1" style="text-align:center;">-647.835</td> <td rowspan="1" colspan="1" style="text-align:center;">-670.035</td> <td rowspan="2" colspan="1" style="text-align:center;">-672.487</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">Disper.-corr.</td> <td rowspan="1" colspan="1" style="text-align:center;">-4.28051</td> <td rowspan="1" colspan="1" style="text-align:center;">-1.02088</td> <td rowspan="1" colspan="1" style="text-align:center;">-0.80754</td> <td rowspan="1" colspan="1" style="text-align:center;">-2.452</td> </tr> <tr> <td rowspan="2" colspan="1" style="text-align:center;">coulomb</td> <td rowspan="1" colspan="1" style="text-align:center;">Coulomb(SR)</td> <td rowspan="1" colspan="1" style="text-align:center;">51184</td> <td rowspan="1" colspan="1" style="text-align:center;">28962.4</td> <td rowspan="1" colspan="1" style="text-align:center;">22365</td> <td rowspan="1" colspan="1" style="text-align:center;">-143.4</td> <td rowspan="2" colspan="1" style="text-align:center;">2743</td> </tr> <tr> <td rowspan="1" colspan="1" style="text-align:center;">Coul.-recip.</td> <td rowspan="1" colspan="1" style="text-align:center;">27988.2</td> <td rowspan="1" colspan="1" style="text-align:center;">12809.1</td> <td rowspan="1" colspan="1" style="text-align:center;">12292.7</td> <td rowspan="1" colspan="1" style="text-align:center;">2886.4</td> </tr> </table>由于计算的染料是都是阴离子, 相互之间形成聚集体, 虽然短程静电相互作用为负值, 但长程静电相互作用为正值, 总的静电相互作用能为正值, 并且总相互作用能也为正值.