杂谈自由能计算，PMF，伞形抽样，WHAM

以两个分子结合过程为例，那么结合能与结合自由能有什么差别？“结合能”很好理解，真空中2个分子离的远远的，干材烈火地碰到一起，体系总能量降低，释放的能量就是“结合能”，它是分子间吸引作用的表现（分子间永久偶极或瞬时偶极引起）。结合能就是内能的改变量。但它显然不是“结合自由能”。那什么是“结合自由能”？学过物化的立刻会想到自由能的定义里有“熵”这个东东——无序度。是的，结合前，2个分子爱去哪去哪儿，自由度多大阿，现在搞在了一起，去哪都得一起，无序度显然降低了啊。如果反应在溶液中进行，不光参与结合的分子的无序度减低，还会影响周围的环境分子。让我把这两个分子仍到水里，离开一定距离。假设这俩分子是苯，它们是疏水的，所以两个人在水里也会很快干材烈火地搞在一起。这个过程显然是个自发过程，前后能量有改变，体系总能量仍然降低。虽然这个能量包括了环境能量的变化，那么这个能量变化是不是结合自由能？不是，它仍然是内能的变化，能量减低以热的形式释放。但是，由于结合前后，这俩搞基的苯却对周围水分子的排列产生了影响，引起了熵变。自由能变化=内能变化+熵能变化， ΔA=ΔU-TΔS 。为什么这里是负号呢？苯苯结合前，黑白参杂，熵肯定大；现在黑是黑，白是白，熵显然降低。上帝老人家显然是不喜欢有序的，至少对平衡态是这个态度，所以熵降低，对应的熵能却要增加，很不利啊。所以公式里TΔS是负的. 可是你会问，熵S这么一个玄密的东西，实验怎么测，模拟怎么算？先别忙，下面再详述。这里先总结一下结合自由能ΔA和结合能ΔU的区别：ΔA=ΔU-TΔS。即，结合能ΔU是结合自由能ΔA的一部分，另一部分是熵能TΔS变化。

上面的叙述是对NVT体系来说的，因为体积不变，不存在反抗外压做功PΔV。如果是NPT体系，那么盒子体积要变化，这时的自由能称为吉布斯自由能，多出一项对外做功PΔV，即ΔG=ΔU+PΔV-TΔUS。如果把前面两项变化和起来，称为焓变ΔH=ΔU+PΔV。这是题外话。那么前面NVT体系的自由能变化称为什么呢？对，亥姆霍兹自由能。下面对这两项不加区别，统称自由能，具体是什么自由能取决于你模拟使用什么系综，NVT还是NPT.

在继续往下说之前，量化版跑gauss的兄弟问了：那前面真空中两分子结合自由能怎么算？其实，你那计算都是0K下的，热力学第三定律说0K下熵为0，不存在熵能变化，因此，你那结合自由能就等于结合能。这里还要补充一点，Gauss中还可以计算2个分子非0k下的结合自由能，那是怎么算的呢？那是要把分子的振动，转动，平动等对熵有贡献的项，按照二次项展开，计算它的S与温度的关系，从而可以求出一定温度下的自由能变化。扯远了，回来接着说我们最关心的问题：怎么模拟一个过程的自由能变化。

但是在开始之前，我必须先举一个例子，这个例子能帮我们时刻清醒地认识体系放热/做功什么的并不等于全部的自由能变化：想象把一块盐扔到一杯水里，这时候盐很快溶解了。有的盐溶解要放热，有的盐溶解要吸热。吸热意味着体系的能量要升高。请问，为什么溶解要吸热的过程也能自发进行呢？

答案在前面已经说的很明白了：能够自发进行的过程，说明体系的总自由能是降低的；而自由能分为2部分，内能和熵能；虽然吸热导致内能升高，但只要熵能降低的足够多，最终自由能仍然是降低的，所以能够自发进行。那这个熵能降低的来源是什么呢？很明显，原先阴阳离子结合在一起，只有一个状态；溶解后能够游来游去，活动空间大大增加，意味着总状态数大大增加，熵也大大增加啊。。。但是要注意的是，虽然离子的熵是增加的，水分子呢？由于一部分水分子被离子束缚在它的壳层内，活动能力反而下降，因此，部分水分子的熵能稍有升高。但是，不影响大局。

现在你应该已经彻底明白自由能，内能，熵能这些概念的区别了。下面切入到主题，怎么用分子模拟来模拟一个过程的自由能变化？比如，一个底物分子与蛋白分子的结合？或一个H2与沸石的吸附？或者一个小分子穿过双层膜？这些都是变化过程。凡是变化过程，都能够用一个叫反应坐标的东西，把这个过程从头到尾串起来。举例来说，比如小分子穿过某一通道，那么通道的中轴可以定义为反应坐标。这个例子比较简单，对于某些复杂过程很难定义反应坐标。比如，你给很多个H2在沸石表面的吸附定一个反应坐标试试？你可以只关注一个H2，但你很难确定它要吸附的方位点。

前面说了，自由能里面有个叫熵的东西，这个东西能模拟吗？不能。熵的定义，叫“一个体系所有的可及微观状态数Ω的对数”， S=k*lnΩ， Ω这个东西，你没法数。另外我们都知道，分子模拟一步算一个采样，假如一个采样对应一个状态的话，等于说你让跑几步它就有几个状态。。。所以，分子模拟不能模拟熵。用严格的话来讲，应该说“熵不是体系微观状态的热力学平均量”，或者“熵是与相空间的体积有关的量”。只有那些有对应微观热力学量的宏观性质，才能模拟之，平均之，比如体积，比如内能，比如压强。

好了，你能模拟内能，你却不能模拟熵，那你怎么计算自由能？？？

兄弟，在说这个问题之前，我们来看看实验怎么测定一个反应前后的自由能差的。怎么测的，还记得吗？。。测浓度变化。。有了反应平衡后的浓度分布，可得平衡常数（产物浓度^m/反应物浓度^n, m,n是配平系数，你懂的），有了平衡常数K，然后∆G°＝-RT ln K。。。可是这个魔术是怎么变的呢？关键在这个平衡常数，它体现了体系中的分子的选择：是选择反应物态，还是产物状态。。如果向往产物态，K值很大；如果向往反应物态，K值较小。。。这个向往，体现了分子“哪儿凉快哪儿呆着去”的本能抉择，以及“哪自由能低哪儿就去哪儿”的自由意志。。。所以最终浓度，反应的是一个由于自由能落差而产生的在两种状态之间的分布状态不同，即概率密度！！如果你看懂了，恭喜你，往下理解就简单了。分子模拟要模拟自由能，正要从分布概率入手！！！

分子模拟怎么模拟自由能变化？简单说，假如体系有A,B 2个状态，你让体系自己去跑。。。跑跑跑。。。跑完一统计，在A态待了20%时间，在B态待了P2=80%时间，概率分别为PA=0.2, PB=0.8, 好，平衡常数K=P2/P1=4, AB两态的自由能差是∆G°＝-RT lnK=-0.6*ln4 ≅ 0.8 kcal/mol (T~300K, 一个kT大约0.6kcal/mol, 你懂的) 这就是分子模拟模拟自由能变化的原理。

有2点需要补充。

首先真实体系不可能只有A, B两个状态。沿着反应坐标，比如两分子相互靠近过程中，有无数个状态。如果你把反应坐标离散化，切成一小段一小段，那么状态数取决与你怎么切。你可以放开体系让它自己在每一段上自由地跳来跳去，最后你统计它落在每一段上的舞步数目，最终可以得到一个概率密度分布，最后用公式转换成每一小段（状态）上的自由能，连成一条起伏的曲线，叫PMF，全称叫Potential of mean force。它，就是传说中的自由能势能面。对，能量变化的曲线或曲面，叫势能面；自由能势能面，叫PMF。为什么取这么个令人无限纠结的名字呢？因为这两个分子靠近的过程中，整个体系其它所有分子对他俩都有影响（作用力），最后得到的是个平均效应，所以叫Potential of mean force。反正wiki上是真么解释的。把这个东西限定在两分子相互靠近的过程。但我觉的这背后的思想可以应用于任何变化过程，只要你能定义反应坐标。然后坐标上每个点的概率，也体现了所有分子对这一点的平均作用。对吧？

其次，概率密度怎么转化成PMF数值呢？回想一个经典的例子，那就是模拟两个甲烷分子在水里面的结合自由能势能面：放开这俩分子让它们在水里跑，最后计算它们之间的径向分布函数g(r)，又叫RDF, 而RDF恰好就是两分子相距不同距离的概率密度。然后按照下式计算每一点的PMF值：

W(r)=-kTln

W(r)就是结合自由能沿反应坐标，即二者之间距离r变化的优美曲线。通常你会在一定距离找到一个最小点，那个点就是他俩最happy的距离。上面这个公式就是所谓的Kirkwood公式。现在你应该已经隐约明白怎么回事了：自由能是跑出来的概率。那这个概率与什么有关呢？反应坐标势能面U(r)。为什么呢？因为模拟体系访问一个点，是按照一定概率访问的，这个概率取决于这个点的势能U，势能越高，访问概率越小。你可能已经立刻想到了Boltzmann分布这个东东，对的，

P = k*exp(-U/kT)

注意，这个概率P是整个体系访问某个能量状态的概率，我们要求的概率，是在某个反应坐标上的概率，已经把体系分为了反应部分和环境部分，然后关心的是反应部分不同态间的相对概率。即，假设用r表示反应坐标

P(r) = k* exp(-U(r)/kT)

通过这个公式，我们现在已经成功地把内能势能面转化U(r)为了概率P(r)，概率再转化为自由能势能面W(r)，倒来倒去，熵就包含在了其中！！！神奇吧？关键就在于Boltzmann分布。分布函数是势能面和自由能势能面之间的纽带。『补注：关键是Boltzmann公式中的kT， kT以热驱动的形式，将分子从它喜欢的低势能状态推向它不喜欢的高势能状态，造成相空间中的可及状态数增加(相空间体积增加)；温度越高，kT推动下的状态数增加就越多。熵就是这么包括进来的』

那好，你会说，现在简单了，只要我们知道沿反应坐标的分布概率，就能求出目标反应体系的自由能势能面。 那怎么求反应坐标分布概率呢？一个最直接的想法就是把反应坐标分成很多点，然后在每一点上模拟，然后求每一点上的概率密度Pi，最后转成自由能势能面。现在，关键是求每一点上的分布概率。

在前面两甲烷的例子中，我们是放开来跑的。现在我要告诉你一个不幸的消息，由于分子模拟本身的缺陷，对于大部分变化过程，这种“放开跑”的方法都不适用。

范伟：怎么个情况？这种“放开跑”的方法，要求体系对反应坐标上的每个点都sample的很好。那我问你，你的模拟会自动访问所有的点吗？一般不会。为什么？正如Boltzmann公式所示，模拟体系访问一个点的概率，与这个点的体系总能量E的指数次方成反比!!! 伤不起啊。我们的模拟体系，绝大部分时间老老实实呆在能量低的构象里不思进取。结果采样不足，这可咋做统计啊亲们？这时有些兄弟疑惑了，能量高的构象，模拟体系爬不上去，那真实的体系就能爬上去吗？。。。能，一个构象能量高，不是说永不可及，而是概率低，但只要时间足够长，体系就能爬上去呆一会儿，这就是为什么有的反应要1-2天，就是中间要穿越很高的能垒，以至穿越概率较低。。。升高温度能增加穿越概率。。。一切动力，都依赖于小小的热运动kT=0.6kcal/mol去创造鲤鱼跳龙门的奇迹啊亲们。。

你能模拟1-2天的时间吗？不能，可悲的研究僧们只能跑几个ns的模拟，以至体系根本无法访问反应坐标上的高能垒部分，他们该肿么毕业呢? 这时候伞形抽样横空出世了。。。苦逼青年有福了。。。

Umbrella Sampling, firstly proposed by G.M. Torrie, J.P. Valleau from Canada at 1977, 提供了一个巧妙的方法: 你这个体系不愿意待在我反应坐标的高能量的段段上不是？哥我用绳子把你栓在那里。具体做法是用一个弹簧势w=k*(x-x0)^2施加在反应物上，一旦反应物坐标x不在我设定的x0上，势能就会升高，你就得乖乖回来。哥把反应坐标切成很多小段(学名叫window)，在每一小段上都这样跑上一会儿，这叫有偏采样， biased sampling，暴力方式强迫体系沿反应坐标前进。可是，我们能接受一个扭曲的结果吗？不能啊。有偏采样的结果，怎么才能转换回无偏采样的概率密度呢？这需要一个“三步上篮”的过程，叫WHAM。

WHAM,  Weighted Histogram Analysis Method, 加权柱状图分析法，目的是把有偏采样的统计结果转换为无偏采样的统计结果，基本思路是这样的：

第一，首先按照下式依次求第i个window上的无偏概率Pi(x)，得到N个window区间上的无偏概率密度。

               exp
Pi(x)  =  ---------------- * Pib(x) 
           

 Pi(x), 无偏概率，即w=0时真实概率密度; i 是window的序号, x是反应坐标
 Pib(x), 有偏概率, 即w=wi时的概率密度; b-biased
 wi(x), 加在第i个窗口上的有偏势大小：

第二，总体概率密度P(x)并非是各window上概率密度Pi(x)的简单相加，而是线性组合P(x)=sum{ci*Pi(x)}

第三，根据P(x), 按照A(x)=-kT*ln求各点自由能，最终得到PMF。

上面的关键是第一步那个公式的理解，和第二步系数ci的求解。
详细公式看这里：www.cse.illinois.edu/~ericcyr/uiuc/wham.pdf
我不打算鏖战于数学细节，而只想说说它背后的思路:

1. 有偏势w对结果概率密度造成的影响，与有偏势的大小有关：越大，偏的越多；等于0时得到真实概率密度；因此，w=0时的概率密度，与w=w的概率密度，他们之间呈比例关系，比例系数则与有偏势能w的Boltzmann分布权重exp(wi)有关；可见学好概率是多玛的重要；

2. 有偏势w对结果概率密度造成的影响，还与i所处的位置有关，能量低的地方，收到的影响就小，能量高的地方，受到的影响就大，因此第二步ci的计算将考虑这个差别。

好了，伞形有偏抽样，以及将有偏采样结果转换成无偏概率密度并得到自由能势能面就是这样的过程。用Gromacs，用NAMD，用CHARMM， AMBER都能跑这样的计算。

伞形抽样代表了自由能计算的一个大的分类。另外一个大的分类，叫热力学积分，包括什么微扰啦，缓慢生长啦。。。热力学积分的思想，与有偏采样有很大的区别。有偏采样上面说啦，基于“概率密度-->自由能”这样一个关系，而热力学积分，则是对势能函数的微扰。它设计一个逐步微扰的过程，使得通过参数λ从0变到1，则体系从初态变到终态。然后:

∆A=A(λ=1)-A(λ=0)=∫<∂U(λ)/∂λ>dλ

稍具数学基础的可以看出，在每个λ上模拟，并通过前后λ步模拟可以得到∂U(λ)/∂λ的数值一阶导，然后积分，pretty easy。

本质上讲，通过引入参数λ，将自由能差这个体系的不可测微观量转变了∂U(λ)/∂λ这个可观测量进行积分。然而，由于基于热力学积分的模拟方法需要设计一个转变路径，要求体系前后两个状态之间差别不大，比如甲烷缓慢长成乙烷以求得二者溶剂化自由能差。然而，你用甲烷长成高分子那就不行了。所以热力学积分法的应用局限于非常特定的例子。

热力学积分引入参数λ，进行逐步平衡态模拟，与前面说的将反应坐标分段非常类似。然而，二者有本质不同：自由能与∂U(λ)/∂λ存在积分关系，然而自由能对反应坐标没有任何的直接关系。这也正是两大类方法的分水岭。
下面关于Umbrella sampling的叙述并不清晰，下面是更新：

原内容：
"Umbrella Sampling, firstly proposed by G.M. Torrie, J.P. Valleau from Canada at 1977, 提供了一个巧妙的方法: 你这个体系不愿意待在我反应坐标的高能量的段段上不是？哥我用绳子把你栓在那里。具体做法是用一个弹簧势w=k*(x-x0)^2施加在反应物上，一旦反应物坐标x不在我设定的x0上，势能就会升高，你就得乖乖回来。哥把反应坐标切成很多小段(学名叫window)，在每一小段上都这样跑上一会儿，这叫有偏采样， biased sampling，暴力方式强迫体系沿反应坐标前进。可是，我们能接受一个扭曲的结果吗？不能啊。有偏采样的结果，怎么才能转换回无偏采样的概率密度呢？这需要一个“三步上篮”的过程，叫WHAM。"

补充更正：

其实，不只是在势能高的地方加势能w，在反应坐标的每个window上，我们都加势能w。但是，原本势能低的构象，
反应坐标x就在此widow的x0附近就很舒服，晃动不那么大，那么计算出来的w就很小；在势能高的构象上，它老是想
往势能低的地方滑，不过由于弹簧拴着，所以跑不走，但弹簧给拉的很长，即x-x0很大，那么计算得到的w当然也就
很大了。最后各window点的概率嘛，按照该第i个window上的平均势能Ui计算，Pi(biased)=<exp(-Ui/kT)>。这个Ui除
了包括这个体系的原子间作用势能之和之外，还有那个弹簧势wi。这肯定不行阿，按上面说的话，是扭曲了正确的概
率分布，或者使得正确的概率密度曲线变形了。肿么办呢？你上面不是有每次采样的拉力w吗，你可以统计得到在这个
window上采样的平均<wi>，这个就是把概率拉变形的罪魁势能，把这个<wi>从平均Ui里面减去就行了阿。
http://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=2&ved=0CCoQFjAB&url=http%3A%2F%2Fwww.pci.tu-bs.de%2Fagelstner%2Flehre%2Fthc4-08%2Fumbintro.
pdf&ei=eGPsTo7-DcXV0QHd6syzCQ&usg=AFQjCNG5QWlhezBkYStVk9Si5svofMSMYA&sig2=AScf_ONUy_E-d0QHQAqRSQ
http://emuch.net/html/201110/3672663.html

欢迎关注微信公众账号