GROMACS QM/MM教程2:编译设置及简单运行

类别:    标签: gmx   阅读次数:   版权: (CC) BY-NC-SA

上一篇博文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, 可能出现下面的错误

  1. QMMM is currently only supported with cutoff-scheme=group

    如果你的mdp中没有设置cutoff-scheme=group, 就会提示这个错误. 按要求使用此选项就可以.

  2. QM/MM does not work in parallel, use a single rank instead

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

  3. 无法找到环境变量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, 故不使用, 置空或随意设置一个路径即可.
  4. 高斯无法运行

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

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

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

技术细节

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

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

其他高斯环境变量

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

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

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

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

  1. 优化后的坐标(非必需, 当工作类型为优化时才需要)
  2. 能量
  3. QM原子上的梯度
  4. 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电荷之间的自相互作用能.

◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: GROMACS QM/MM教程1:QM/MM方法的实现
后一篇: 自动调整VMD窗口的位置和大小

访问人次(2015年7月 9日起): | 最后更新: 2024-01-20 10:40:28 UTC | 版权所有 © 2008 - 2024 Jerkwin