Gromacs结合自由能计算教程（二）
以下文章来源于网络，如有侵权,请联系管理员删除!
如何计算自由能？
续上回讲到的，上一章简要介绍了自由能计算的相关背景与炼金术热力学途径。这一节将会讲解计算炼金术自由能中每一步热力学过程所涉及的自由能计算方法。
首先回顾一下炼金术热力学的途径，我们会涉及到一个讲研究对象从其环境中解耦和的过程，在这个过程中我们会先后且逐步地关闭与环境的静电作用与范德华作用，保证相邻模拟间尽可能大的相空间重叠，降低自由能计算的不准确性。再次强调这个过程是分开进行的，并且静电作用总是先于范德华力解耦以避免正负电荷过于近的不合理接触（范德华力的交换互斥项可以避免这种不合理接触）。
另外，正如你此时可能会感到疑问的一点是，我在这里所述的范德华力与我们大学课本中所讲述的似乎略有出入。
在大学课本中的范德华力是基本上是这么说的：
“范德华力的主要来源有三种机制：取向力：极性分子之间的相互作用力诱导力：极性分子极化非极性分子，产生诱导偶极矩，并相互作用色散力：一对非极性分子由于本身电子的概率运动，可以相互配合产生一对方向相反的瞬时偶极矩，这一对瞬时偶极矩相互作用。这种机制是非极性分子中范德华力的主要来源，1930年由F.W.伦敦首先根据量子力学原理给出解释“ ——from Wikipedia
到目前为止我所说的范德华力实际指的大约（说是大约是因为这样说也不准确）只有第三项——色散力，严谨的表达是分子动力学软件用来描述原子间作用的所谓L-J势，包含了近距离的交换互斥项与远距离的吸引项。搞清楚分子动力学中的范德华力这个基础概念十分重要，我会在之后的文章中单独论述。事实上许多中文教材在范德华力的描述相当old school，存在着自相矛盾与纠缠不清的地方。
回到炼金术自由能计算。我们在这里引入了一个变量λ来线性地描述解耦过程，Gromacs提供了一系列关键字如“couple-lambda”等来实现相关操作（前提是使用关键字“free-energy = yes”）。在热力学积分法中，我们计算自由能的方法实际上是使用∂U/∂λ对λ进行积分，λ实际取的是离散值，所以我们得到的自由能实际也是近似值；原则上如果λ无限细分，所得的自由能就可以无限逼近真实值。
如上图就是先后解耦库伦作用与范德华力整个过程的热力学积分示意。
事实上按照这个炼金术途径处理自由能的方法有很多种，在之后的分析结果会进行进一步的讨论。值得一提的是，为了保证∂A/∂λ函数的连续性，在关闭范德华力时我们总是会使用“软核势能”来描述原子，着可以实现真正意义上的“逐步关闭原子相互作用”。
Gromacs中关于自由能的参数设置
自由能计算在Gromacs中主要涉及以下额外关键词
free-energy = yes
couple-moltype = ligand
couple-lambda0 = none
couple-lambda1 = vdw-q
couple-intramol = no
separate-dhdl-file = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
init-lambda-state = 0
coul-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.00 0.00 0.0 0.1 0.2 0.4 0.5 0.6 0.8 1.0
vdw-lambdas = 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.85 0.9 0.95 0.97 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
nstdhdl = 100
calc-lambda-neighbors = -1
Couple-moltype是你要耦合（或解耦）的分子（基团）名称
Couple-lambda0是你期望的起点状态，Couple-lambda1是终点状态；在给出的例子中，Couple-lambda0=none，Couple-lambda1=vdw-q意思是在初始状态不与环境相互作用，在最终状态完全打开范德华力与静电作用，所以着定义的是一个耦合的过程。如果是解耦过程应该反过来写。
couple-intramol =no 指的是在这个过程中不处理分子内去耦合/耦合的问题
sc-alpha, sc-power, sc-sigma三项是关于软核作用势的设定，需要参考软核作用势的相关文献。不知道软核作用势的情况下就用这套值也没太大问题。
init-lambda-state 指的是当前模拟的λ值。事实上在一个计算中会涉及十几个或更多λ下的计算。可以使用脚本mk_mdp.sh 一键生成不同lambda的mdp文件，给不同过程调用。
vdw-lambdas与coul-lambda对应每个λ下的范德华力与静电力的耦合量，0代表没有相互作用，1代表完全打开。这些值的设定有一定技巧性，需要根据模拟体系进行调整，之后会进一步讨论。在这里是一个耦合过程，所以先打开了范德华作用，再打开静电作用。
nstdhdl = 100指的是每100帧保存一次∂H/∂λ数据
calc-lambda-neighbors = -1指的是计算并输出所有lambda值的ΔH供MBAR分析。通常情况下，在设置了init-lambda-state时，calc-lambda-neighbors = 1指的是计算当前lambda下前后一个lambda值对应的ΔH，对于普通BAR分析已经足够了，而对于MBAR分析，需要使用-1，表示计算所有值。关于MBAR分析之后会提到。
