小分子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电荷在一定程度上可行,具体参数设置还需要适当的调整。用该服务器的一个重要的好处就是电荷可以采用高斯计算,而发表文章时候引用该网站原始文献即可,不会涉及版权纠纷。

Ambertools+Gaussian计算拟合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版本的ambertools和gaussian09W可以在参考资料中找到。具体还需要在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表示该分子结构不完整或者其他问题。完整输出如下图所示:

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原子电荷的超级懒人脚本(一行命令就算出结果)