小分子resp电荷拟合

介绍了几种拟合有机小分子的RESP电荷的途径。

R.E.D服务器计算拟合RESP电荷

原始文献为:RED Server: a web service for deriving RESP and ESP charges and building force field libraries for new molecules and molecular fragments,至今引用次数已达423
网站为:https://upjv.q4md-forcefieldtools.org/REDServer-Development/
可能该网站看上去花里胡哨的,让人摸不着头脑,但实际上该网站提供电荷的计算拟合,以及力场的开发等等。当然,我现在只是需要用到它来拟合RESP电荷,那就简单多了!

运行流程

其实当你把它理解以后,整个工作流程也就看上去简明了。主要是根据用户提供的计算设置参数以及小分子结构文件,提交给服务器,然后通过它的计算程序(包括Gaussian03,Gaussian09,Gaussian16,GAMESS-US或者Firefly各版本的计算程序)来完成后面的步骤。当然要想使用高斯程序,你必须先注册该网站的私人账号,要不然你就只能使用免费的计算程序。注册过程也十分简单,邮箱即可,用于通知。

实例

这里,我以三聚氰胺分子单体来介绍一下该如何使用R.E.D服务器。
为了使用该服务器,我们必须提供两类文件,将其打包成压缩包上传服务器。

  • 分子结构文件
    必须是PDB格式的文件,这一点不够人性化。另外单个分子原子数上限为225
  • 参数文件

    • Project.config
      该文件包括了一些项目配置参数。包括分子名称,分子总电荷,自旋多重度等等,以下是我的一些参数情况:
      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      MOLECULE1-TITLE = melamine 
      !! 分子名称
      MOLECULE1-TOTCHARGE = 0
      !! 分子总电荷为0
      MOLECULE1-SPINMULT = 1
      !! 分子自旋多重度1
      MOLECULE1-ATMREORDR = OFF
      !! 关闭原子顺序自动排列
      MOLECULE1-ATMCOR = ON1
      !! 只重命名多余的原子名
      MOLECULE1-CALCONECT = OFF
      !! 关闭自动确定原子间的成键,此时提供的PDB文件需要包含连接信息
      MOLECULE1-CALCHEMEQ = ON
      !! 使用等价
      MOLECULE1-RESMOD = ON
      !! 自动校正残基名称

    更多参数可以查看Project.config,以上基本可以了。

    • System.config
      该文件包括了一些计算参数,以及拟合电荷的方式选择等等,以下是我的情况:
      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      LANGUAGE = EN 
      !! PyRED log文件的语言,设置为英语即可
      MAXMEMVAL = 12288
      !! 最大可用内存,单位MB
      QMSOFT = GAUSSIAN
      !! 采用高斯计算
      NP = 8
      !! 并行计算处理器核心数为8
      CHR_TYP = RESP-X1
      !! 电荷拟合方式为 B3LYP/6-31G(d)
      SURFMK_MEPCALC = IOp(6/33=2,6/41=6,6/42=6)
      !! IOp关键字

    更多参数可以查看System.config

完成以后,将System.configProject.configmelamine.pdb三个文件一起打包成zip或者rar格式压缩包,然后一步步提交到R.E.D服务器,提交成功会有一封邮件信息发送到你的邮箱,包含你的任务号。并提供了任务进程面板,比如下图红框中为我的任务情况:

任务计算完成以后,也会有一封邮件通知。可以打包下载所有输出文件,其中...charge.txt文件为拟合的原子电荷。和自己用高斯在相似基组下拟合的电荷情况相比较,差别不算太大。

结论

R.E.D服务器拟合RESP电荷在一定程度上可行,具体参数设置还需要适当的调整。用该服务器的一个重要的好处就是电荷可以采用高斯计算,而发表文章时候引用该网站原始文献即可,不会涉及版权纠纷。

Gaussian+Ambertools计算拟合RESP电荷 <常用>

Ambertools16结合Gaussian 09 D.01/E.01,就用以下步骤就行了,不要做额外的修改(除了理论方法和基组)。6/50=1是不需要的。
(1)产生Gaussian计算MK电荷的输入文件:
antechamber -i lig.mol2 -fi mol2 -o lig.gjf -fo gcrt
(2)用Gaussian执行:
g09 < lig.gjf > lig.out
(3)做RESP电荷计算并产生分子的库文件(会从lig.out中读取分子附近格点上的静电势数据):
antechamber -i lig.out -fi gout -c resp -o lig.prepin -fo prepi

另外,我也提供了一个能在Windows平台计算的脚本,和acpype在线网站差不多,但是电荷方式不同。主要是利用ambertools的antechamberacpypegaussian09W在Windows10上计算resp电荷,脚本点击下载。计算的是B3LYP/6-311G**级别下的电荷。Windows版本的ambertoolsgaussian09W可以在参考资料中找到。具体还需要在Windows10环境变量中添加:

1
2
3
4
5
ABERHOME=E:/User-software/amber18
Path=E:\User-software\amber18\bin
GAUSS_EXEDIR=D:\G09W
GAUSS_SCRDIR=D:\G09W\Scratch
Path=D:\G09W

此处用的Cygwin终端,其他支持bash的终端也是可以的,比如Cmder。使用方式也很简单,只需要将小分子pdb结构文件 (不能用中文名字)和resp.bsh脚本放在同一文件夹下,命令运行bash resp.bsh等待任务完成即可。此脚本支持一次性读取多个单独pdb结构文件,按顺序完成所有计算工作。注意检查终端输出内容,特别是Errors字段,不为0表示该分子结构不完整或者其他问题。完整输出如下图所示:

GAMESS量化软件+resp工具(Ambertools自带)

也可以通过GAMESS软件计算,然后通过resp工具拟合电荷,具体步骤如下(以水分子为例):

  • GAMESS输入文件water.inp如下:
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
     $CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE DFTTYP=B3LYP ICHARG=0 MULT=1 $END
    $SYSTEM MEMORY=100000000 $END
    $STATPT NSTEP=500 $END
    $PCM SOLVNT=WATER $END
    $BASIS GBASIS=N311 NGAUSS=6 NDFUNC=1 NPFUNC=1 DIFFSP=.TRUE. DIFFS=.TRUE. $END
    $ELPOT IEPOT=1 WHERE=PDC OUTPUT=PUNCH $END
    $PDC PTSEL=CONNOLLY LAYER=4 VDWINC=0.2 MAXPDC=100000 $END
    $DATA
    Water
    cnv 2

    O 8
    H 1 .000000 .926627 -.239987
    $END

此输入文件表明是在B3LYP/6-311++G**进行计算,水分子的H采用对称性。

  • GAMESS计算。在Windows下执行:
    1
    rungms.bat water.inp 2016-pgi-linux-mkl 4

计算结束后会产生water.dat输出文件,里面包含了拟合电荷所需数据。

  • water.dat文件中提取相关数据,用作resp工具输入文件。原csh脚本。在这里我提供一个我修改后的bash脚本如下:

    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
    107
    #!/bin/bash
    # a script to generate resp input file from gamess dat file

    ########################
    #Modified by liuyujie #
    #Date 2019.01.09 #
    ########################

    if [ $# -lt 1 ];then
    echo "Usage:$0 gamess.dat"
    exit 1
    fi

    if [ ! -e $1 ];then
    echo "Input file $1 not exist!"
    exit 2
    fi

    # judge whether dat file from gamess is right for fitting resp charges
    opt='$ELPOT IEPOT=1 WHERE=PDC OUTPUT=PUNCH $END'

    # judge whether dat file includes esp data
    if [[ `grep "POTENTIAL COMPUTED FOR" $1`a == a ]];then
    echo "Your $1 does not appear to have ESP data"
    # judge file suffix
    file=$1
    if [[ "${file##*.}" == "dat" || "${file##*.}" == "DAT" ]];then
    echo "Add the line below to ${file%.*}.inp and rerun"
    echo "$opt"
    else
    echo "Perhaps, you are specifying a GAMESS output (not punch) file?"
    fi
    exit 3
    fi

    # define output file
    file=$1
    resp=${file%.*}.in
    esp=${file%.*}.esp

    # create .in file for resp
    cat > $resp << EOF
    TITLE
    &cntrl nmol=1, ihfree=1
    &end
    1.0
    EOF

    # find molecule name from dat file
    n=`grep -n "DATA" $1 |head -1 |cut -d: -f1`
    let n++
    sed -n -e "${n}p" $1 >> $resp

    # find molecule net charge and the number of atom
    na=`grep "POTENTIAL COMPUTED FOR" $1 | tail -1`
    ch=$(echo $na |cut -d' ' -f9)
    na=$(echo $na |cut -d' ' -f5)
    echo "$ch $na" |awk '{printf "%5d%5d\n", $1, $2}' >> $resp

    # find atomic number of molecule and write into .in file
    n=`grep -n "POTENTIAL COMPUTED FOR" $1 | tail -1 | cut -d: -f1`
    let n0=$n+$na
    head -n $n0 $1 |tail -n $na |awk '{printf "%5d%5d\n", $2, 0}' >> $resp
    echo 0 |awk '{printf "\n\n\n\n"}' >> $resp

    # create .esp file for resp
    # find points number and write data into .esp file
    npt=`grep "POINTS NPT=" $1 |tail -1`
    npt=$(echo $npt |cut -d' ' -f7)
    echo "$na $npt 0" | awk '{ printf "%5d%5d%5d\n",$1,$2,$3 }' > $esp
    head -n $n0 $1 | tail -n $na | awk -v s=\ '{ printf "%16c%16.7e%16.7e%16.7e\n",s,$3,$4,$5 }' >>$esp

    # write ELPOTT, x, y, z into .esp file
    n=`grep -n "POINTS NPT=" $1 | tail -1 | cut -d: -f1`
    let n1=$n+$npt
    head -n $n1 $1 | tail -n $npt | awk '{ printf "%16.7e%16.7e%16.7e%16.7e\n",$5,$2,$3,$4}' >>$esp

    # finish
    echo "RESP files were created"
    echo "Example of usage:"
    echo -e "For first fitting without equivalence : \n\t resp -i $resp -e $esp -o ${file%.*}.out -p ${file%.*}.pch -t ${file%.*}.chg"

    # Whether do second fitting
    file2=$2
    file2suffix=${file2##*.}
    if [[ $# == 2 && ${file2suffix} == "txt" ]];then
    head -6 $resp > ${file%.*}2.in

    # modify .in file for second fitting
    # add iqopt, that is, read in new initial charges from -q unit
    sed -i 's/ihfree=1/&, iqopt=2/' ${file%.*}2.in

    # atom number
    head -n $n0 $1 |tail -n $na |awk '{printf "%5d\n", $2}' > tmp.in

    # custom ivary for atom equivalence
    cat *.txt |awk '{printf "%5d\n", $1}' > tmp2.in

    # merge atom number and ivary and append *2.in file
    paste tmp.in tmp2.in |awk '{printf "%5d%5d\n", $1, $2}' >> ${file%.*}2.in
    echo 0 |awk '{printf "\n\n\n\n"}' >> ${file%.*}2.in
    rm -rf tmp.in tmp2.in

    # finish
    echo ""
    echo -e "For second fitting with equivalence : \n\t resp -i ${file%.*}2.in -e $esp -q ${file%.*}.chg -o ${file%.*}2.out -p ${file%.*}2.pch -t ${file%.*}2.chg"
    fi

    使用方式为:

    1
    bash makeresp.bsh water.dat [*.txt]

会得到water.espwater.in(如果提供等价性.txt文件,则还会生成wate2.in用于第二步拟合)输出文件,这些文件作为后续的输入文件。

  • resp工具拟合。第一次不使用等价性:
    1
    resp -i water.in -e water.esp -o water.out -p water.pch -t water.chg

会得到water.outwater.chgwater.pch三个文件,其中后两个文件中就包含resp电荷值。

  • 如果你要使用原子等价,自己另外需提供一个后缀名为.txt的文件,比如以乙醇分子为例,内容如下:
    1
    2
    3
    4
    5
    6
    7
    8
    9
    0
    0
    2
    2
    0
    0
    6
    -1
    -1

0表示不使用等价,2表示该行原子与第二个原子等价,6同理;-1或者-99表示第二步拟合任然使用第一步的电荷。
然后运行:

1
resp -i water2.in -e water.esp -q water.chg -o water2.out -p water2.pch -t water2.chg

即可得到二次拟合以后具备等价性的RESP电荷了。

说明:这里其实你没必要为了使用resp单独去安装Ambertools,你可以直接用Fortran编译器编译这个源文件(里面包含一个编译好的Windows版本)即可,Windows下相关命令为:

1
ifort /fpp /DCYGWIN resp.F lapack.F

Multiwfn计算拟合

.fch+Multiwfn

利用gaussian计算产生的chk另存为fch文件,然后利用Multiwfn中的RESP计算模块拟合原子电荷。

Multiwfn自己算

Multiwfn自身有计算静电势的代码,计算小体系耗时不太高,但是体系大了的话耗时会比较高。另外,Multiwfn作者在新版本中也提供了一次性计算RESP电荷的脚本,见参考资料4。
注意 pop=mk或者CHELPG IOP等关键词不要和opt一起用,这样Multiwfn可能识别错误。

不同计算方式计算拟合的电荷有些差距。

参考资料

1.论坛-sobereva
2.RESP拟合静电势电荷的原理以及在Multiwfn中的计算
3.李继存老师编译的Win版本的AmberTool:
(1)AmberTools-win下载
(2)win下的环境变量设置
4.计算RESP原子电荷的超级懒人脚本(一行命令就算出结果)