请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 920|回复: 4

[Multiwfn资源与经验] 使用Multiwfn做基于分子力场的能量分解分析

[复制链接]

1万

帖子

25

威望

1万

eV
积分
34922

管理员

公社社长

发表于 2018-9-19 17:23:29 | 显示全部楼层 |阅读模式
使用Multiwfn做基于分子力场的能量分解分析

文/Sobereva @北京科音
First release: 2018-Sep-19  Last update: 2018-Sep-21


摘要:本文介绍量子化学波函数分析程序Multiwfn (http://sobereva.com/multiwfn)中的基于分子力场的能量分解的原理和使用。此方法耗时极低,结果物理意义明确,不仅在一些场合可以替代较昂贵的基于波函数的能量分解方法,而且还具有一些独特优势,诸如结果可以可视化、可方便地考察分子内弱相互作用等。在Multiwfn使这种分析变得非常简单后,在未来必将得到大量应用。对Multiwfn不了解的话建议参看《Multiwfn入门tips》(http://sobereva.com/167)、《Multiwfn波函数分析程序的意义、功能与用途》(http://sobereva.com/184)。


1 基于力场的能量分解的原理

能量分解是量子化学分析方法中的一个重要组成,它可以把片段间总相互作用能分解为有物理意义的能量项以便于考察相互作用的本质,这在《Multiwfn支持的弱相互作用的分析方法概览》http://sobereva.com/252)一文中有简要介绍。实际上,基于形式非常简单的分子力场(forcefield),也可以对弱相互作用的成份进行拆解。虽然分子力场已被广泛用于分子动力学等领域,但是基于分子力场做能量分解研究的文章极少。考虑到通过力场做能量分解的潜在应用价值,笔者将之写入了Multiwfn程序中作为主功能21的子功能1,经实测发现效果不错。

分子力场里所谓的非键相互作用包括静电作用(electrostatic)和范德华作用(van der Waals),后者可划分为起排斥作用的“交换互斥项”(repulsion)和起吸引作用的“色散项”(dispersion)。如果你对原子间相互作用本质缺乏了解,很建议读一下《谈谈“计算时是否需要加DFT-D3色散校正?”》(http://sobereva.com/413)中的相关介绍。大多数分子力场采用对势形式来计算原子间非键相互作用的静电(ele)、交换互斥(rep)和色散吸引(disp)部分,公式如下所示
1.png
其中A,B是原子标号,q是原子电荷,r是原子间距离,ε是范德华作用势阱深度,R0是原子间非键距离,当r=R0时原子间范德华作用能恰等于势阱深度。

ε和R0是在分子力场里定义的。不同力场里根据原子所处化学环境的不同定义了不同原子类型,多数力场是根据原子类型来定义的范德华参数,不同原子类型的范德华参数有的相同有的不同。实际计算原子间范德华作用时用的参数一般是基于原子的德华参数通过混合规则来产生的,比如UFF力场用的是下图第一行的规则,不同原子类型间的ε和R0取两个原子类型的ε和R0的几何平均。而对于AMBER和GAFF力场,规则如下图第二行所示,ε混合规则与UFF相同,但是力场对每个原子定义的是非键半径参数R*,原子间的非键距离参数通过相应原子的非键半径加和得到。
2.png
如上可见,计算原子间非键相互作用能非常简单。要想计算特定片段间的非键相互作用,就把片段间每一对原子的非键作用加和即可。这就是基于力场做能量分解的原理,很容易理解。Multiwfn目前支持流行的AMBER99、GAFF和UFF力场做能量分解,这三种力场的函数形式都是相同的,都是最简单形式的分子力场,静电作用仅通过原子电荷来体现,并没有像一些更复杂力场那样考虑原子多极矩或者可极化效应。(以后Multiwfn也可能会支持这些更复杂的力场形式)

Multiwfn中基于力场的能量分解功能相对于Morokuma、SAPT等主流的基于波函数的能量分解方法有以下优势
·速度非常快,对很大的体系也只要一瞬间就能给出结果。而基于波函数的能量分解方法很难用到一百多个原子的体系(有的只能用到几十原子体系)
·片段可以非常自由地划分,想考察哪两个区域的相互作用就相应地输入原子编号来定义即可
·分子内弱相互作用可以方便地考察。而对于基于波函数的能量分解方法,考察分子内弱相互作用需要把分子拆成片段,其中牵扯很多麻烦的事情,具体细节上也很具有任意性。

然而基于力场的能量分解也有很多局限性:
·受制于分子力场简单的函数形式,能量分解精度和基于波函数的方法还有差距,不能在定量精度上有太高要求
·不同力场支持的元素不同,要考察的片段里面有力场不支持的元素就没法用
·计算原子电荷的方法存在一定任意性,而且有的弱相互作用注定没法基于原子电荷来定性正确表现,比如卤键是由于卤原子有sigma-hole特征所产生的,这需要考虑原子四极矩或者非原子中心的点电荷才能描述
·给出的结果没法体现片段间相互作用所产生的电子云的极化、电荷转移等因素的影响。
·只能考察弱相互作用,而没法考察化学键这样的强相互作用(这方面即便是SAPT等很多基于波函数的能量分解方法也没法考察)

做这种能量分解所用的原子电荷必须对分子范德华表面附近的静电势有很好的重现性,否则不可能靠原子电荷较可靠地计算片段间静电相互作用。因此使用的原子电荷首选是拟合静电势电荷,其中最常见的是MK和CHELPG电荷。如果对拟合静电势电荷不了解,参看《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》(http://sobereva.com/441)和《原子电荷计算方法的对比》(http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml)。所用原子电荷对静电势重现性哪怕达到“不错”的档次都太不够,比如用ADCH、CM5、AM1-BCC等电荷,虽然对静电势重现性也不错,仅次于拟合静电势电荷,但是算出来的静电相互作用能和拟合静电势电荷还是有一定差距的。

关于能量分解时力场的选择,如果是有机体系,就用AMBER或者GAFF就挺好。这两个力场的原子类型的名字虽然不同(一个很大差别是AMBER的是大写,GAFF的是小写),但实际上GAFF的范德华参数是从AMBER力场继承过来的,主要差别在成键参数上,因此基于哪个力场做能量分解都可以。虽然Multiwfn的能量分解也支持UFF力场,但一般不建议用UFF做能量分解,如后文的例子所示,对于哪怕合理的结构,UFF也往往严重高估互斥作用,甚至令本身稳定结合的分子间相互作用为正。不过将结构用UFF力场优化后可以消除这个问题,但由于UFF力场优化出的弱相互作用体系的结构往往也不太理想,所以能量分解结果还是不如AMBER和GAFF好。UFF力场最大好处是支持元素非常广泛,几乎覆盖了整个周期表。

用于做能量分解的几何结构应当是通过靠谱的级别进行优化过的,如果是来自于X光衍射晶体结构,应当冻结重原子的位置而对氢原子进行优化,因为氢原子的位置准确度很差,往往只是经验性地确定的。



2 基于力场的能量分解在Multiwfn中的使用

其实基于力场做能量分解,利用诸如GROMACS等分子动力学就能做,但是真要用起来,过程很麻烦。Multiwfn专门设计了一个模块来做这种能量分解,在设计上已经尽可能把分析流程简化了(但注定不可能做到完全傻瓜化)。基本使用流程如下:

(1)对体系中每类分子都创建一个文本文件(分子类型文件),里面包含分子中各个原子的原子类型和原子电荷。然后再写一个文本文件(分子列表文件),里面是分子类型文件的路径,还有相应类型分子的数目,详见后文。
(2)启动Multiwfn,读取一个含有当前体系结构信息的Multiwfn可以支持的输入文件。哪些文件被Multiwfn支持且带有结构信息,在《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》。(http://sobereva.com/379)里有明示,比如你用.fch、.mol、.pdb、.xyz、.wfn、.molden等等格式都可以
(3)进入主功能21(能量分解主功能)中的子功能1(基于力场的能量分解模块)。
(4)用选项3定义读取分子列表文件,从而给体系中各个原子分配原子类型和电荷。如果你要用的力场不是默认的AMBER & GAFF,则读取前先选选项-1切换成要用的力场。
(5)用选项2来定义片段,先设置要定义的片段数,然后输入各个片段包含的原子序号,输入格式很灵活,比如3,5-9,12-20,69。第(4)步和第(5)的顺序无所谓,可以互换。
(6)最后,选择选项1开始进行分析。各个片段间的相互作用能及物理成分,以及片段里各个原子对片段间相互作用能的贡献都会输出到屏幕上。

在写分子列表文件和分子类型文件的时候一定要*万分注意*,把分子列表文件和其中分子类型文件按顺序完全展开之后,原子的顺序必须和Multiwfn启动时载入的文件里的原子顺序完全一致