结合自由能计算：MD vs. QM

nachtmusik

本文首发于我的知乎专栏https://zhuanlan.zhihu.com/siaolang

结合自由能的计算可以使用基于量子化学的方法，也可以使用基于分子动力学模拟（MD）的方法，然而两者的原理大为不同。量子化学（QM）的方法考虑了实际电子结构，原则上说在计算相互作用焓的时候应该比做了大量近似处理的分子力学精确的多，但缺陷在于其巨大的计算成本；只能对于某一个或者几个期待的构象进行静态的研究，难以对构象变化带来的熵进行计算——量化中熵的计算仅能考虑振转熵贡献，基于hessian矩阵的计算得到。

相比之下分子力学能够对体系的构象变化进行充分的采样，并且可以处理大如蛋白的体系，使用显式溶剂模型——显式模型下才能考虑到某些复杂的能量变化，例如疏水空腔的“高能水”释放，离散溶质效应……

相信很多人和我一样一直抱有疑问，既然MD在底层计算上精度不够，QM又难以进行充分的统计学采样，二者在实际计算中，到底那种更准确？

这里就要讲到一篇文献Hyrophobicchallenge, 文中比较了几种QM方法和两种MD方法对于预测疏水中性客体与葫芦脲的结合自由能与实验值的差值。

文中所使用的主客体如下图，其中DAPI为实验测量键合常数所用的荧光探针
对于MD部分，使用牵拉法计算自由能

MD1：使用OPC水模型，从R.E.D网站获得的RESP电荷结合GAFF2力场生成拓扑。

MD2：使用TIP3P水模型，使用AM1-BCC电荷结合GAFF力场生成拓扑。

原则上MD1的精度应该是要高于MD2的，无论是水模型还是电荷、力场大概都要优于MD2。


对于QM部分，自由能为电子能（内能）、振转熵贡献、溶剂化自由能三步独立计算之和。也许是为了比较与MD的水和能计算差异，QM是在真空下基于真空优化的结构做单点和振转分析的，然后再单独计算水和能的贡献（除了QM4考虑了溶剂对构象的影响）。这种做法我个人并不推荐在实践中这么做，真空中构象和溶剂中的构象往往会有很大差异，尤其是对于带电分子这么做就真是很不明智了；正如之后结果也会提到的，QM4因为考虑了溶剂对构象的影响，实际结果更为准确。如果条件允许，构象优化与结合能、振转分析都在SMD水溶剂模型下做，不用特意拆开来去算溶剂化自由能，可能会更准确。

几何构型优化与振转熵贡献均在PBEh-3c下完成。作者采用了四种不同的方法计算结合能：

QM1：色散矫正meta-GGATPSS-D3ATM/def2-QZVP，水和能在COSMO-RS 12fine模型下
QM2：色散矫正meta-hybridPW6B95-D3ATM/def2-QZVP'（对于重原子舍去f轨道，对于氢原子舍去g轨道），水和能在COSMO-RS 13模型下
QM3：色散矫正DLPNO-CCSD(T)/CBS*，水和能在COSMO-RS 13模型下
QM4：色散矫正DLPNO-CCSD(T)/CBS*，为了考虑溶剂对最优构象的影响，使用XTB+GBSA预优化，然后在PBEh-3c+DCOSMO-RS进一步优化的结构下做单点，水和能在COSMO-RS 13模型下
特别地，对于QM1, QM2，为了改善QM方法对构象采样不足的弱点，还使用了QMDFF（Prof. Grimme自己开发的基于量化计算hessian矩阵生成的分子力场）在优化构象进行模拟退火，得到了几种能量最低构象，然后在PBEh-3c级别再优化， 一并进行相同级别QM计算，最后进行玻尔兹曼加权，以考虑构象熵贡献。对于QM3、QM4，由于量化方法的高精度耗时极高，只考虑单个构象能量。

最后的结果如下，纵轴代表理论计算结果，横轴为实验结果，理想状态下所有点都应该落在对角线上。

总体而言，几种QM方法计算的结合自由能与实验数据有着更小的误差，但是和实验数据的相关性较差；而MD方法总体与实验数据有着较大的误差，对于结合能大的体系误差更大，但与实验数据有良好的相关性。

这个结果可以说在意料之中，它们的误差是由其固有的缺陷造成的——如前面所述，MD的底层理论简陋，QM无法进行充分采样。

MD1与MD2的差距十分有限，这也可以理解，下面会谈到溶剂化能其实各种方法都算的挺准的，而GAFF2对GAFF的改进有限，AM1-BCC虽然略显粗糙但也有不错的准确度，综合下来看就是预测上基本没什么不同。

再来看QM的数据，几种方法共同的一点是，对于简单碳氢化合物1~5预测基本都在回归线上，但对于19、20、21、17，前两者为六元环，后两者为较长链烷烃，此时计算值就与回归线有较大偏差 （同为长链的26和六元环的25误差也较大，但只在QM4中出现），这是可能是因为在这几种分子实际结合时的构象较为复杂（苯环看似简单，实则存在经式纬式构象，并且容易旋转），量化方法难以充分采样的缺点在此暴露了出来。另外一点可以作证这个观点的证据是，采用了构象加权矫正的QM1、QM2，虽然与实验数据偏差较大，但与实验数据的相关性好了一些，这和MD的情况类似。

至于水和能的计算方面，COSMO-RS13（红线）与MD（蓝线）方法都与实验吻合得很好，看来主要的误差是在主客体结合能的计算与熵贡献方面。

总体看来，这一轮较量之中MD是处于下风。但MD就真的不适合计算自由能了吗？

其实不然，首先要考虑到在这组研究中所使用的主体分子是葫芦脲，一种较为刚性的分子，它本身的构象变化较少，因此本身比较利于量化计算；并且例子中太多分子都是结构比较简单的分子，但即便如此还是有个别分子存在构象问题造成的偏差。对于实践中常常处理的蛋白质空腔以及其客体分子可能会有更大的构象复杂度，这时候量化方法的偏差可能会更大。

另一方面，MD的好处还是在于其较小的计算开销，能够处理显式溶剂、大体系（蛋白质）、长时间采样，在显卡加速和多核计算加持下计算速度远远优于高精度QM，尤其是对于蛋白质体系，QM几乎是无法处理的，虽然也有多精度组合方法的QM，但是操作过于复杂，怎样划分不同精度处理的界限也是难题。所以我个人认为，对于蛋白质与配体的结合自由能处理应该首先考虑MD方法。

文中还有更多有意思的内容，对此感兴趣的请务必读一下原文：
https://pubs.acs.org.ccindex.cn/ ... .jpcb.7b09175#t1fn5

naoki

正有此疑惑，学习了，谢谢~

fhh2626

MD作者用的是attach-pull-release方法，不是伞状采样

APR理论上很严谨，但是主客体相互作用计算的话alchemistry才是标准方法，APR能不能代表MD还有待讨论

k64_cc

fhh2626 发表于 2019-5-25 11:40
MD作者用的是attach-pull-release方法，不是伞状采样

APR理论上很严谨，但是主客体相互作用计算的话alch ...

万一数据量堆上来的话应该也差不太多…？

fhh2626

k64_cc 发表于 2019-5-25 11:55
万一数据量堆上来的话应该也差不太多…？

其实我觉得这个例子的话应该差不多。。。但是如果蛋白质结构比较复杂，配体又插得比较深（好污）的话，pull那一步有点unphysical

nachtmusik

fhh2626 发表于 2019-5-25 11:40
MD作者用的是attach-pull-release方法，不是伞状采样

APR理论上很严谨，但是主客体相互作用计算的话alch ...

已修改
另：我个人做过几例类似主客体体系，使用alchemical方法，也存在对于强键合客体有高估结合自由能的趋势（弱键合客体相对准确）；我觉得这应该是力场问题，尤其是静电能计算方面高估得比例比较大。

nachtmusik

fhh2626 发表于 2019-5-25 12:18
其实我觉得这个例子的话应该差不多。。。但是如果蛋白质结构比较复杂，配体又插得比较深（好污）的话，pu ...

没错，这是APR的缺陷。
一方面需要设定一个足够大的模拟盒子，耗时会增大很多；另一方面，尤其是对于处理蛋白质-配体系统时，牵拉路径未必是简单的直线，而可能要穿过曲折的通道，这种情况下牵拉路径是极难定义的。

fhh2626

nachtmusik 发表于 2019-5-25 12:31
已修改
另：我个人做过几例类似主客体体系，使用alchemical方法，也存在对于强键合客体有高估结合自由能 ...

这个就是经典力场最大的问题，感觉目前也没有什么好办法解决

tianyongpan

fhh2626 发表于 2019-5-25 11:40
MD作者用的是attach-pull-release方法，不是伞状采样

APR理论上很严谨，但是主客体相互作用计算的话alch ...

“主客体相互作用计算的话alchemistry才是标准方法”，请问这种方法有文章讲述嘛，刚接触这一块内容

fhh2626

tianyongpan 发表于 2019-5-26 10:43
“主客体相互作用计算的话alchemistry才是标准方法”，请问这种方法有文章讲述嘛，刚接触这一块内容

自由能微扰（FEP）或者热力学积分（TI），你随便一搜就一大堆

nachtmusik

tianyongpan 发表于 2019-5-26 10:43
“主客体相互作用计算的话alchemistry才是标准方法”，请问这种方法有文章讲述嘛，刚接触这一块内容

我的专栏有介绍这种方法
https://zhuanlan.zhihu.com/p/60963446
里面提到了几个教程和分析方法，你可以顺着这些链接看下去。
alchemistry只是一种方便计算的热力学途径，分析最后的结合自由能有又很多种方法，包括自由能微扰，热力学积分，MBAR等等。

JCenter

如果是跑BOMD 100ps左右。精度相比于DFT静态计算和经典力场MD下的效果，会更好吧。这种情况下为了节省算力可以使用（QM/MM方法，水分子使用经典力场，主客体还是考虑电子，从而降低待计算的电子数目）

CrysW555

感谢楼主分享，原帖文献打不开，是否是 https://pubs.acs.org/doi/full/10.1021/acs.jpcb.7b09175 这个?