- 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
~/.bashrc | |
---|---|
1 2 3 | export GMX_QM_GAUSS_DIR=E:/G09w
export GMX_QM_GAUSS_EXE=g09.bsh
export GMX_QM_MODIFIED_LINKS_DIR=''
|
简单运行
上面的设置都成功后, 就可以拿一个简单的体系来测试运行了. 可惜的是, GROMACS的QM/MM模块已经很老, 且无人更新, 坐等废弃, 所以你直接拿来运行时不大会成功的. 可行的解决方案就是自己写点代码, 将高斯的运行过程包装一下, 来代替直接使用GROMACS运行高斯. 上篇博文中已经给了这么一个bash脚本gau
, 但写得比较乱, 不易扩展修改. 我整理改进了一下, 并加以注释. 想了解具体细节, 遇到问题或进行修改时, 可参考注释.
g09.bsh | |
---|---|
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 | #/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
|
将此脚本g09.bsh
放于高斯可执行程序所在目录下, 再将环境变量GMX_QM_GAUSS_EXE
设为g09.bsh
, 然后准备一个高斯关键词文件_Gkwd.gjf
, 大致如下
_Gkwd.gjf | |
---|---|
1 2 3 4 | %chk=input
%mem=1GB
#T RHF/STO-3G NoSymm Force Punch=Derivatives
Charge Prop=(Field, Read)
|
这样就可以开始测试运行了.
简单示例
-
两个水分子, 一个作为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电荷之间的自相互作用能.