博文
GROMACS QM/MM教程2：编译设置及简单运行
|
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 | exportGMX_QM_GAUSS_DIR=E:/G09wexportGMX_QM_GAUSS_EXE=g09.bshexportGMX_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对输入和输出文件进行转换##首先对GMXQM/MM模块给出的高斯输入文件进行处理,生成新的输入文件#然后调用系统已安装的高斯程序运行新的高斯输入文件#最后将高斯的输出文件转换为GMX需要的格式,供其调用##**注意**系统可能对长文件名支持不好,所以此脚本名称最好短一点#调用高斯程序的完整路径,根据自己的安装进行设置
GAUSSIAN=E:/G09w/g09.exe#清空前一步运行用到的临时文件
rm-rf_GMXtmp.*fort.7gxx.*#GMXQM/MM给出一个"正常"的高斯输入文件#里面包含用户在mdp文件中给出的一些设置#但这个文件格式古老且死板,不适用#因此你需要自己提供一个_Gkwd.gjf,放于工作目录下#里面指定计算时要使用的关键词(【标题】之前的部分)#这些关键词会覆盖GMX的默认设置##**注意**#1.脚本坚持在输入文件中使用埃为单位,所以你不能再使用Units=Bohr#2.脚本会根据情况自动添加Guess=Read,所以你不能再使用Guess=Read#3.脚本中必须使用NoSymmForcePunch=Derivatives三个关键词#4.如果体系中有MM电荷,必须使用ChargeProp=(Field,Read)关键词#5.其他关键词可根据自己的情况决定使用与否
[[!-e_Gkwd.gjf]]&&{echo'NO_Gkwd.gjf!!!';exit;}#使用自定义的关键词
awkNF_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.7fn",$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.7fn","",$1*Bhr,$2*Bhr,$3*Bhr}'input.com>>_GMXtmp.gjf#高斯输入文件准备完毕,调用高斯进行计算#如果运行时还需要其他文件(如态平均系数),可在此处添加$GAUSSIAN_GMXtmp.gjf_GMXtmp.out#处理得到的高斯输出文件,准备GMX需要的输入文件
awk'#获取MM电荷值/PointCharges:/{Nchg=0;getlinewhile($1=="XYZ="){Nchg++;Q[Nchg]=$6getline}}#获取MM电荷自能,SCF能量/Selfenergyofthecharges/{Eslf=$7}/SCFDone:/{Escf=$5}#获取每个MM电荷上的电场强度/--------ElectricField--------/{getline;getlineNpnt=0;getlinewhile(NF>4){if($2!="Atom"){Npnt++Ex[Npnt]=$3;Ey[Npnt]=$4;Ez[Npnt]=$5}getline}}#输出需要的能量和梯度END{printf"%19.12fn",Escf-Eslfsystem("catfort.7")for(i=1;i<=Nchg;i++)printf"%20.10e%20.10e%20.10en",-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#TRHF/STO-3GNoSymmForcePunch=Derivatives
ChargeProp=(Field,Read) |
这样就可以开始测试运行了.
简单示例两个水分子, 一个作为QM, 一个作为MM
CH4+与H2碰撞, 全部作为QM, 可能会反应
下载相关输入文件
说明 使用QM方法模拟反应时, 很难一次就能跑出一条反应轨迹. 这是正常的. 因为每个反应都有特定的发生几率, 不可能随便发生. 否则相应的反应速率非常大, 与事实不符. 如果你随便一跑就得到一条反应轨迹, 除了偶然之外, 更大的可能是你的设置有问题. 此外, 同样的初始条件(构型, 速度), 使用不同的QM方法进行模拟, 得到的轨迹也大不同. 一种QM方法下的反应轨迹, 换用另一种QM方法很可能就变成非反应轨迹了. 所以, 讨论反应速率的话, 需要对大量的轨迹进行统计, 单单几条轨迹意义不大. 当然, 根据几条反应轨迹说明反应的一些特点和机制, 还是可行的.
技术细节GROMACS QM/MM模块涉及的源代码位置.srcgromacsmdlibqmmm.h
.srcgromacsmdlibqmmm.cpp
.srcgromacsmdlibqm_orca.cpp
.srcgromacsmdlibqm_mopac.cpp
.srcgromacsmdlibqm_gamess.cpp
.srcgromacsmdlibqm_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水分子中每个原子上的力.
知晓了此文件的格式, 自己写其他量化程序的外挂脚本也就不困难了. 值得注意的是, 计算能量时一般要扣除MM电荷之间的自相互作用能.
◆本文地址: http://jerkwin.coding.me/2017/02/02/GMX_QMMM教程2-编译设置及简单运行/, 转载请注明◆
https://blog.sciencenet.cn/blog-548663-1031074.html
上一篇：GROMACS QM/MM教程1：QM/MM方法的实现
下一篇：自动调整VMD窗口的位置和大小
