仓库源文，站点原文

layout: post title: GROMACS QM/MM教程2：编译设置及简单运行 categories:

科 tags:
gmx

2017年02月02日 21:00:10

上一篇博文GROMACS QM/MM教程1: QM/MM方法的实现中, 整理了GROMACS QM/MM方法的基本实现原理, 略显理论化, 本教程中具体讲解如何操作. 以下说明基于Windows下的GROMACS 5.1.4版本.

编译设置

编译GROMACS时, 要打开QM/MM功能选项. 一般是使用高斯--with-qmmm-gaussian, 因为后面一般是使用自己的外挂脚本, 所以就选这个我比较熟悉的. 当然, 如果你更熟悉MOPAC, ORCA或GAMESS, 也可以使用这三个, 但使用MOPAC必须修改源码, 使用GAMESS太复杂, 因此建议选择高斯或ORCA之一.

编译完成后, 在mdp中使用QM/MM的相关选项, 运行grompp, 可能出现下面的错误

QMMM is currently only supported with cutoff-scheme=group

如果你的mdp中没有设置cutoff-scheme=group, 就会提示这个错误. 按要求使用此选项就可以.
QM/MM does not work in parallel, use a single rank instead

GROMACS的QM/MM模块无法并行, 只支持单核. 运行时只能使用 gmx mdrun -nt 1. 但根据我的经验, 有时没有加-nt 1选项, mdrun仍可以运行, 详细情况有待考察.
无法找到环境变量GMX_QM_GAUSS_DIR, GMX_QM_GAUSS_EXE, GMX_QM_MODIFIED_LINKS_DIR

这三个环境变量是GROMACS调用高斯时使用的, 必须设置, 否则就无法启动高斯.

GROMACS最终在命令行中调用高斯使用的命令为

【GMX_QM_GAUSS_DIR】/【GMX_QM_GAUSS_EXE】 < input.com > input.log

所以你必须设置好需要的环境变量, 保证在命令行中上面的命令能运行成功. 三个环境变量的具体说明和设置如下:
- GMX_QM_GAUSS_DIR: 高斯安装路径, 也即高斯的可执行程序所在的目录. 遵循Linux下路径规则, 使用斜线分隔符, 但注意最后无斜线. Windows下类似E:/G09w, Linux下类似/路径/g09.
- GMX_QM_GAUSS_EXE: 可执行程序/脚本的名称, 如g09.exe, 或者自己的脚本g09.bsh. 注意不能包含路径.
- GMX_QM_MODIFIED_LINKS_DIR: 修改后的高斯LINKS路径. 因为我们不打算修改LINKS, 故不使用, 置空或随意设置一个路径即可.
高斯无法运行

高斯运行时可能需要自己的环境变量 GAUSS_EXE_DIR, 需要的话, 按要求设置. 有时还需要将高斯路径加到path环境变量.

下面是不同系统的环境变量设置示例

Windows

Linux

<table class="highlighttable"><th colspan="2" style="text-align:left">~/.bashrc</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" style="background: #f8f8f8"><pre style="line-height: 125%">export GMX_QM_GAUSS_DIR=E:/G09w export GMX_QM_GAUSS_EXE=g09.bsh export GMX_QM_MODIFIED_LINKS_DIR='' </pre></div> </td></tr></table>

简单运行

上面的设置都成功后, 就可以拿一个简单的体系来测试运行了. 可惜的是, GROMACS的QM/MM模块已经很老, 且无人更新, 坐等废弃, 所以你直接拿来运行时不大会成功的. 可行的解决方案就是自己写点代码, 将高斯的运行过程包装一下, 来代替直接使用GROMACS运行高斯. 上篇博文中已经给了这么一个bash脚本gau, 但写得比较乱, 不易扩展修改. 我整理改进了一下, 并加以注释. 想了解具体细节, 遇到问题或进行修改时, 可参考注释.

<table class="highlighttable"><th colspan="2" style="text-align:left">g09.bsh</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 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">#/bin/bash # # 请仔细阅读此脚本, 了解具体的处理过程 # 因为你可能需要对其修改才能适用于自己的体系 # # 此脚本对 高斯GAUSSIAN 程序进行封装 # 主要使用awk对输入和输出文件进行转换 # # 首先对GMX QM/MM模块给出的高斯输入文件进行处理, 生成新的输入文件 # 然后调用系统已安装的高斯程序运行新的高斯输入文件 # 最后将高斯的输出文件转换为GMX需要的格式, 供其调用 # # **注意** 系统可能对长文件名支持不好, 所以此脚本名称最好短一点 # 调用高斯程序的完整路径, 根据自己的安装进行设置 GAUSSIAN=E:/G09w/g09.exe # 清空前一步运行用到的临时文件 rm -rf _GMXtmp.* fort.7 gxx.* # GMX QM/MM给出一个"正常"的高斯输入文件 # 里面包含用户在mdp文件中给出的一些设置 # 但这个文件格式古老且死板, 不适用 # 因此你需要自己提供一个 _Gkwd.gjf, 放于工作目录下 # 里面指定计算时要使用的关键词(【标题】之前的部分) # 这些关键词会覆盖GMX的默认设置 # # **注意** # 1. 脚本坚持在输入文件中使用 埃 为单位, 所以你 不能 再使用 Units=Bohr # 2. 脚本会根据情况自动添加 Guess=Read, 所以你 不能 再使用 Guess=Read # 3. 脚本中 必须 使用 NoSymm Force Punch=Derivatives 三个关键词 # 4. 如果体系中有MM电荷, 必须 使用 Charge Prop=(Field, Read) 关键词 # 5. 其他关键词可根据自己的情况决定使用与否 [[ ! -e _Gkwd.gjf ]] && { echo 'NO _Gkwd.gjf !!!'; exit; } # 使用自定义的关键词 awk NF _Gkwd.gjf >> _GMXtmp.gjf awk '/guess=read/{print "Guess=Read"}' input.com >> _GMXtmp.gjf echo "" >> _GMXtmp.gjf # 使用GMX给出的【标题】【电荷 多重度】【QM坐标】【MM坐标 电荷大小】 awk ' BEGIN { n=0; Bhr=0.52917721092 } NF==0 { n++ } n==1 && NF >0 { print } n==2 && NF==4 { printf"%4d %16.7f %16.7f %16.7f\n", \ $1, $2*Bhr, $3*Bhr, $4*Bhr } n >2 && NF==4 { printf"%4s %16.7f %16.7f %16.7f %10.6f\n", \ "", $1*Bhr, $2*Bhr, $3*Bhr, $4 } n >1 && NF!=4 { print } ' input.com >> _GMXtmp.gjf # 高斯无法直接给出MM电荷上的梯度 # 为得到MM电荷上的梯度, 需要指定高斯计算MM位置上的静电场 # 然后根据MM电荷的大小和静电场计算出所需要的梯度 # 为此, 需要再次提供MM电荷位置的坐标, 并将其写到文件的最后 # 这些坐标仍然由GMX提供的文件获得 # # **注意** 下面两点脚本会自动处理 # 1. MM电荷部分与坐标部分之间 必须 有空行 # 2. MM电荷部分 必须 以埃为单位 awk ' BEGIN { n=0; Bhr=0.52917721067 } NF==0 { n++; next } n > 2 { printf"%4s %16.7f %16.7f %16.7f\n", \ "", $1*Bhr, $2*Bhr, $3*Bhr} ' input.com >> _GMXtmp.gjf # 高斯输入文件准备完毕, 调用高斯进行计算 # 如果运行时还需要其他文件(如态平均系数), 可在此处添加 $GAUSSIAN _GMXtmp.gjf _GMXtmp.out # 处理得到的高斯输出文件, 准备GMX需要的输入文件 awk ' # 获取MM电荷值 /Point Charges:/ { Nchg=0; getline while($1=="XYZ=") { Nchg++; Q[Nchg]=$6 getline } } # 获取MM电荷自能, SCF能量 /Self energy of the charges/ { Eslf=$7 } /SCF Done:/{ Escf=$5 } # 获取每个MM电荷上的电场强度 /-------- Electric Field --------/ { getline; getline Npnt=0; getline while(NF>4) { if($2!="Atom") { Npnt++ Ex[Npnt]=$3; Ey[Npnt]=$4; Ez[Npnt]=$5 } getline } } # 输出需要的能量和梯度 END { printf"%19.12f\n", Escf-Eslf system("cat fort.7") for(i=1; i<=Nchg; i++) printf"%20.10e%20.10e%20.10e\n", \ -Q[i]*Ex[i], -Q[i]*Ey[i], -Q[i]*Ez[i] } ' _GMXtmp.out > _GMXtmp.7 # 高斯使用 D 而不是 E 作为科学计数法的标识, 改过来, 否则GMX无法识别 sed 's/D/E/g' _GMXtmp.7 > fort.7 </pre></div> </td></tr></table>

将此脚本g09.bsh放于高斯可执行程序所在目录下, 再将环境变量GMX_QM_GAUSS_EXE设为g09.bsh, 然后准备一个高斯关键词文件_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" style="background: #f8f8f8"><pre style="line-height: 125%">%chk=input %mem=1GB #T RHF/STO-3G NoSymm Force Punch=Derivatives Charge Prop=(Field, Read) </pre></div> </td></tr></table>

这样就可以开始测试运行了.

简单示例

两个水分子, 一个作为QM, 一个作为MM
CH4+与H2碰撞, 全部作为QM, 可能会反应

下载相关输入文件

说明使用QM方法模拟反应时, 很难一次就能跑出一条反应轨迹. 这是正常的. 因为每个反应都有特定的发生几率, 不可能随便发生. 否则相应的反应速率非常大, 与事实不符. 如果你随便一跑就得到一条反应轨迹, 除了偶然之外, 更大的可能是你的设置有问题. 此外, 同样的初始条件(构型, 速度), 使用不同的QM方法进行模拟, 得到的轨迹也大不同. 一种QM方法下的反应轨迹, 换用另一种QM方法很可能就变成非反应轨迹了. 所以, 讨论反应速率的话, 需要对大量的轨迹进行统计, 单单几条轨迹意义不大. 当然, 根据几条反应轨迹说明反应的一些特点和机制, 还是可行的.

技术细节

GROMACS QM/MM模块涉及的源代码位置

.\src\gromacs\mdlib\qmmm.h
.\src\gromacs\mdlib\qmmm.cpp
.\src\gromacs\mdlib\qm_orca.cpp
.\src\gromacs\mdlib\qm_mopac.cpp
.\src\gromacs\mdlib\qm_gamess.cpp
.\src\gromacs\mdlib\qm_gaussian.cpp

各部分的功能文件名称揭示得很清楚.

其他高斯环境变量

实际上GROMACS与高斯联用时, 根据源代码还有其他一些非必要的环境变量

GMX_QM_CPMCSCF
GMX_QM_SA_STEP
GMX_QM_ACCURACY
GMX_QM_GAUSSIAN_NCPUS
GMX_QM_GAUSSIAN_MEMORY

实际上这些设置直接在高斯的输入文件中设置更方便, 使用环境变量反而罗嗦.

GROMACS支持的高斯输出文件格式

GROMACS与高斯联用时, 读取的高斯输出文件为fort.7. 这是一个文本文件, 里面的数据为

优化后的坐标(非必需, 当工作类型为优化时才需要)
能量
QM原子上的梯度
MM原子上的梯度

只要此文件满足上面的数据格式, 且数据数量与体系需要的一致, GROMACS就可以读取. 下面是两个水分子的例子, 第一行为体系能量, 下面三行为QM水分子中每个原子上的力, 最后三行为MM水分子中每个原子上的力.

   -74.959189786700
    0.8122428486E-01    0.4936266263E-02    0.2405930866E-01
   -0.6954208610E-01    0.1296978682E-02   -0.7381477271E-02
   -0.1168492817E-01   -0.6264009111E-02   -0.1670717046E-01
    4.5870000000e-05   -3.5862000000e-05   -1.0508400000e-04
   -2.0433000000e-05    3.7947000000e-05    7.2558000000e-05
   -2.2935000000e-05    2.8773000000e-05    6.1716000000e-05

知晓了此文件的格式, 自己写其他量化程序的外挂脚本也就不困难了. 值得注意的是, 计算能量时一般要扣除MM电荷之间的自相互作用能.