分子动力学中结合自由能的认识

第一章 自由能与 “能量”的差别

       分子动力学的模拟体系，不论是NPT，还是NVT，都有一个总能量内能。内能分为动能和势能两项，势能项分解成很多子项，如键能、非键作用能等。模拟是否平衡，总能量是否守恒，需要绘制一下总能可以得知。

      关于自由能：在热力学第二定律中自由能是自发过程的判据，即一个过程的自由能降低就能自发进行；能自发进行的过程，其自由能必然降低，即小于0。

1.1 关于结合能

        以两个分子的结合过程为例，进一步讨论结合能与结合自由能的差别。

        真空中两个分子离的很远，结合在一起，体系总能量将降低，释放的能量就是“结合能”，它是分子间吸引作用的表现（分子间永久偶极或瞬时偶极引起）。结合能就是内能的改变量。在溶液中能自由结合的两个分子，是自发过程，前后能量发生改变，体系总能量降低。尽管包括了环境能量的变化，但其能量变化不是结合自由能，仍然是内能的变化，能量减低以热的形式释放。

1.2 关于熵

        “熵”：是无序度，体系的混乱程度。结合前，分子自由度大，结合后无序度将会降低。在溶液中进行时，参与结合的分子的无序度降低，周围的环境分子也会受影响。两个分子在结合前后，对周围溶液分子的排列产生影响，这将引起了熵变。  

1.3 结合自由能。

            由上可得：结合自由能变化=结合能(内能)变化+熵能变化；例如ΔA=ΔU-TΔS。

           备注：-TΔS，分子结合后，熵降低，即混乱程度降低，对应的熵能增加，所以公式里TΔS是负的。结合能ΔU是结合自由能ΔA的一部分，另一部分是熵能TΔS变化。

        帮助理解部分：能量越低越稳定，稳定状态下的能量会较低。自然界的万物能量的变化是趋于一个自然状态而言的。比如两个分子应该结合，结合后对其描述的总能将较低，释放了结合能。对熵能而言，熵是描述混乱程度，自然界趋于一个无序的状态，结合后无序度降低，不是其稳定状态，其能量自然增加。ΔU是内能，分子结合内能降低，熵能增加，为了保持能量平衡，因此其应该是一个负号来让其保持在一个能量平衡的状态。

1.4 讨论

      以上叙述是对NVT体系讨论，NVT体系的自由能变化称为亥姆霍兹结合自由能A，因为体积不变，不存在反抗外压做功PΔV。若是NPT体系，盒子体积要变化，这时的自由能称为吉布斯结合自由能G，多出一项对外做功PΔV，即ΔG=ΔU+PΔV-TΔUS。如果把前面两项变化和起来，称为焓变ΔH=ΔU+PΔV。

     具体是什么自由能取决于你模拟使用什么系综，NVT还是NPT。G=A+PV(P为压力，V为体积)。在生物的动力学反应中，因为Δ(PV)可以忽略不计，所以两者是相同的，这点需要理解。

     体系放热/做功并不等于全部的结合自由能变化。把一块盐扔到一杯水里，盐溶解，有的盐溶解要放热，有的盐溶解要吸热。吸热意味着体系的能量要升高。溶解要吸热的过程也能自发进行是因为体系的总自由能是降低的；自由能分为2部分，内能和熵能；吸热导致内能升高，但只要熵能降低的足够多，最终自由能仍然是降低的，所以能够自发进行。原先阴阳离子结合在一起，只有一个状态，比较有序；溶解后游来游去，活动空间大大增加，混乱度增加，熵大大增加，熵能即降低，因为稳定态趋于的是无序态。离子的熵增加，其中一部分水分子被离子束缚在它的壳层内，活动能力反而下降，因此，部分水分子的熵能稍有升高。

第二章 分子动力学模拟与结合自由能

2.1.分子动力学模拟时，为什么经常说结合自由能里面的熵不能模拟？

      熵 “一个体系所有的可及微观状态数Ω的对数”， S=k*lnΩ。

      分子模拟跑一步算一个采样，假如一个采样对应一个状态的话，等于说跑几步就有几个状态。所以，分子模拟几乎不能模拟所有的可微状态，因此可及微观状态数Ω没法确定。“熵不是体系微观状态的热力学平均量”，或者“熵是与相空间的体积有关的量”。只有那些有对应微观热力学量的宏观性质才能模拟平均，比如体积，比如内能，比如压强。

2.2 能模拟结合能（内能），却不能模拟熵，怎么计算结合自由能？

                     以下内容来自网络许多资源统计，很久之前摘抄。

        实验通过测浓度变化测定一个反应前后的结合自由能差。通过反应平衡后的浓度分布，可得平衡常数（产物浓度^m/反应物浓度^n），有了平衡常数K，然后∆G＝-RT ln K。关键在这个平衡常数，它体现了体系中的分子的选择：是选择反应物态，还是产物状态。如果向往产物态，K值很大；如果向往反应物态，K值较小。这个向往，体现了分子本能抉择，以及“哪自由能低哪儿就去哪儿”的自由意志。所以最终浓度反应的是一个由于自由能落差而产生的在两种状态之间的分布状态不同，即概率密度！因此分子模拟要模拟自由能，就从分布概率入手。 

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

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

     理解：（1）反应坐标。比如，一个底物分子与蛋白分子的结合；或一个H2与沸石的吸附；或者一个小分子穿过双层膜，这些都是变化过程。凡是变化过程，都能够用反应坐标把这个过程从头到尾串起来。举例来说，比如小分子穿过某一通道，那么通道的中轴可以定义为反应坐标。这个例子比较简单，对于某些复杂过程很难定义反应坐标。比如，很多个H2在沸石表面的吸附定一个反应坐标时可以只关注一个H2，但又很难确定要吸附的方位点。

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

       W(r)就是结合自由能沿反应坐标，即二者之间距离r变化的能量曲线。通常你会在一定距离找到一个最小点，那个点就是他俩最乐观的距离。自由能是跑出来的概率。这个概率与反应坐标势能面U(r)有关。因为模拟体系访问一个点，是按照一定概率访问的，这个概率取决于这个点的势能U，势能越高，访问概率越小。

     Boltzmann分布：P = k*exp(-U/kT)。这个概率P是整个体系访问某个能量状态的概率，我们要求的概率，是在某个反应坐标上的概率，已经把体系分为了反应部分和环境部分，然后关心的是反应部分不同态间的相对概率。即，假设用r表示反应坐标P(r) = k* exp(-U(r)/kT)。通过这个公式，把内能势能面转化得到概率P(r)，概率再转化为自由能势能面W(r)，如此这样，熵就包含在了其中！关键就在于Boltzmann分布。分布函数是势能面和自由能势能面之间的纽带。Boltzmann公式中的kT， kT以热驱动的形式，将分子从它喜欢的低势能状态推向它不喜欢的高势能状态，造成相空间中的可及状态数增加(相空间体积增加)；温度越高，kT推动下的状态数增加就越多。熵就是这么包括进来的。

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

     甲烷的例子中，是放开来跑的。但由于分子模拟本身的缺陷，对于大部分变化过程，“放开跑”的方法都不适用。“放开跑”的方法，要求体系对反应坐标上的每个点都采样的很好。模拟体系不会自动访问所有的点。因为正如Boltzmann公式所示，模拟体系访问一个点的概率，与这个点的体系总能量E的指数次方成反比。模拟体系，绝大部分时间是在能量低的构象里呆着。结果采样不足，无法统计。尽管模拟体系爬不上去，但真实的体系能爬上去，一个构象能量高，不是说永不可及，而是概率低，但只要时间足够长，体系就能爬上去在该构象停留一段时间，这就是为什么有的反应要1-2天，就是中间要穿越很高的能垒。穿越概率较低，依赖于小小的热运动kT=0.6kcal/mol，升高温度就能增加穿越概率。但模拟时不可能模拟一两天，只能跑几个ns的模拟，以至体系根本无法访问反应坐标上的高能垒部分。因此需要增强采样方法。例如伞形采样窗口。伞形抽样代表了结合自由能计算的一个大的分类。

        另外一个大的分类，叫热力学积分，包括微扰，缓慢生长。

        热力学积分的思想，与有偏采样有很大的区别。有偏采样基于“概率密度-->自由能”这样一个关系，而热力学积分，则是对势能函数的微扰。它设计一个逐步微扰的过程，使得通过参数λ从0变到1，则体系从初态变到终态。然后∆A=A(λ=1)-A(λ=0)=∫d。在每个λ上模拟，并通过前后λ步模拟可以得到∂U(λ)/∂λ的数值一阶导，然后积分。本质上讲，通过引入参数λ，将自由能差这个体系的不可测微观量转变了∂U(λ)/∂λ这个可观测量进行积分。然而，由于基于热力学积分的模拟方法需要设计一个转变路径，要求体系前后两个状态之间差别不大，比如甲烷缓慢长成乙烷以求得二者溶剂化自由能差。然而，你用甲烷长成高分子那就不行了。所以热力学积分法的应用局限于非常特定的例子。热力学积分引入参数λ，进行逐步平衡态模拟，与前面说的将反应坐标分段非常类似。然而，二者有本质不同：自由能与∂U(λ)/∂λ存在积分关系，然而自由能对反应坐标没有任何的直接关系。这也正是两大类方法的分水岭。

第三章.初步总结蛋白-配体结合自由能的计算方法

基于几何约束的蛋白质-配体准确结合自由能计算 -付浩浩

      通过以上文献进行总结。谈到结合自由能，首先分类为准确结合自由能计算和近似结合自由能计算，准确即标准。

                                                                准确计算

      “准确”一词有两个方面的含义: 一是经典力场或打分函数可以准确地描述微观粒子的相互作用, 二是模拟过程完整地捕获了蛋白质-配体结合过程各个自由度的贡献, 并准确地描述统计力学行为. 其中前者属于力场拟合范畴的问题, 近年来发展的极化力场对非键相互作用有着较好的描述。从统计力学的角度看, 蛋白质-配体的标准结合自由能可以通过下式计算:其中ΔG为标准结合自由能。KB为玻尔兹曼常数。T为温度。K为标准平衡常数。p（bound）为结合态出现的概率，p（unbound）为非结合状态。

        然而, 直接通过式(1)计算结合自由能在实际中十分困难, 因为在可及的分子模拟时间尺度内, 很难看到结合态和非结合态间的自发转化。

        第一类：应用重要性采样方法(importance sampling), 如metadynamics或自适应偏置力(adaptive biasing force, ABF)可以加速蛋白质-配体的结合-分离过程, 实现在分子模拟的时间尺度内观察到该过程. 然而由于在结合-分离过程中, 蛋白质-配体间还伴随着复杂的取向和构象变化, 而在现有的计算条件下, 重要性采样模拟不可能同时捕获沿平动、转动和构象变化自由度方向上的运动, 使得直接使用重要性采样方法计算蛋白质-配体的结合自由能亦十分困难。

        第二类：另一种标准结合自由能计算的思路为使用Alchemistry 类方法(炼金术方法)。在炼金术方法中，通过定义一系列虚拟的中间态缓慢地将一个配体转换成另一个配体，或将一个配体转换成一组没有相互作用的虚拟粒子。因此，炼金术方法不仅可以计算蛋白与不同配体的相对结合自由能，也可以计算蛋白与配体结合的绝对自由能。如自由能微扰(free-energy perturbation,FEP)和热力学积分(thermodynamics integration,TI)，通过在体系中“消失”或“生成”一部分原子, 收集模拟过程中体系的内能变化, 进一步计算“消失”或“生长”过程对应的自由能变化. 在实际应用中,常通过设计热力学循环, 进行多次alchemistry 计算, 得到目标结合自由能。在该方法中，随着配体的逐渐“消失”, 配体相对于蛋白质会出现复杂的平动、转动和构象变化等过程. 而在分子模拟的时间尺度内难以捕获这些过程, 从而阻碍Alchemistry 类方法计算的收敛性。但是这类方法在计算中仍然是值得学习和关注的。

        第三类：为解决重要性采样或Alchemistry 算法在计算蛋白质-配体结合自由能过程中的收敛性难题, 可以采用几何约束减少需要采样的构象空间, 并通过后处理准确地扣除约束对计算结果的影响。近年来陆续提出了漏斗状(funnel)约束, 球状(spherical)约束, 六/七自由度(six/seven degree of freedom)约束等.这些约束分别成功地与伞状采样(umbrella sampling, US) 、metadynamics、自适应偏置力、自由能微扰、热力学积分等重要性采样或alchemistry 算法联用, 准确计算了大量蛋白质-配体的结合自由能。

                                                                 近似计算

       进行蛋白质-配体的准确结合自由能计算复杂且昂贵. 在许多情况下并不需要得到十分准确的结合自由能, 仅需定性比较结构类似的配体与蛋白质作用能力的相对强弱。分子对接方法采用完全经验的打分函数和全局搜索算法, 可以极其快速的定性判断结合能力强弱。在需要半定量的场合, 采用基于终点(end-point)模拟的方法可以十分便捷地给出近似的结合能力。另外, 近年提出的CL-FEP 方法将中心极限定理与自由能微扰算法相结合, 使自由能微扰公式适用于终点模拟, 是一种非常有潜力的近似结合自由能计算方法。

      端点法（endpoint methods）最粗糙但是计算量最小的方法。由于只考虑蛋白和配体结合与解离这两种末端状态，计算效率显著高于其他三种方法，适合于大通量筛选。但是计算结果往往是这四种方法里最不可靠的，尤其是对于绝对自由能的计算，往往给出非常糟糕的结果。但是对于相同类配体的比较计算却很有用。

      常用的端点法有：Linear interaction energy (LIE)、Molecular mechanics/generalized Born surface area (MM/GBSA)、Molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA)。

第四章 结合自由能与药物的关系

      基于结构的药物设计方法有多种多样，包括分子对接、数据库搜索以及从头设计等等，但是它们最终往往都面临一个同样的问题就是药物分子活性的评价。许多药物和其它生物分子的活性都是通过与受体大分子之间的相互作用表现出来的，所以受体和配体之间结合自由能(binding afinity)评价是基于结构的计算机辅助药物分子设计的核心问题。

     精确的自由能预测方法能够大大提高药物设计的效率。在过去的20年中，随着受体-配体相互作用的理论研究以及计算机辅助药物分子设计方法的快速发展，自由能预测方法的研究受到了越来越多的关注。

     在配体(包括药物)和受体的相互作用过程中，牵涉到两类相互作用，即非键相互作用和共价相互作用。非键相互作用包括静电相互作用、范德华相互作用以及氢键相互作用等等，这些相互作用都可以通过力场计算进行比较简单的表达;而共价相互作用则牵涉到化学键的断裂和生成，需要用量子力学的方法进行考察。

   在药物分子和受体的相互作用时，一般只牵涉到非键相互作用，而共价相互作用只存在少数体系中。

第五章 结合自由能计算的发展与困难

      准确计算蛋白质-配体的结合自由能被认为是药物设计领域的圣杯。近二十年来，计算化学领域的研究人员提出了不同的方法，并成功应用于不同领域。

     然而，进行准确结合自由能计算过程十分复杂，且实际操作中有可能会遇到各种技术问题（如质子化态的选择、自由能计算的收敛问题、数据的记录与后处理、结合过程中结合位点与外界水交换的处理等），因此，对于非该领域的专业人士而言，实现蛋白质-配体标准结合自由能计算依旧十分困难。

     之后将为自由能专门做一个part系列。此外，之后的这些part系列，我会根据学习的进度和知识的深入，不断改进和提升。望关注。