NAMD计算自由能

Review NAMD计算自由能

bay__gulf618 (bay__gulf618@sina.com)

College of Chemistry, Chemical Engineering and Materirls Science, Soochow University 前言: 自由能的求算是分子模拟中最重要,也是最困难的工作之一. 重要是因为 . 本文介绍了4种常用的计算自由能的方法, 并详细讲述了各种方法的具体实现. 计算生化体系的自由能主要有如下4种方法，这些在namd 中都已得到实现 1 probability densities, 2 nonequilibrium work,

3 free energy perturbation (FEP), 4 thermodynamic integration (TI)

在namd2.7 版中对自由能的module进行了大幅度改进，可以使用复杂的collective variables 进行probability densities,计算, 从而实现在键角,二面角,RMSD等各种构象指标变化下的PMF。而在2.6版本中仅能实现距离变化时候的PMF, 本文只介绍2.6 中probability densities 的方法，就是ABF。

FEP 和TI 在思想和运行方式上极为相似，只是自由能的定义和处理方法略微不同。在namd ug 2.7 中，两种方法是放在一起介绍的，本文中也是放在一起介绍。

本文所用教程主要有 (1) NAMD ug 2.6/2.7 (2) Stretching 10Ala:

http://www.ks.uiuc.edu/Training/Tutorials/science/10Ala-tutorial/ (3) ABF-Mar2008

http://www.ks.uiuc.edu/Research/namd/tutorial/abf/ (4) C. Jarzynski, PRL, 78 (1997) 2690 (5) S. Park, JCP, 119 (2003) 3559

的还有伞形取样(umbrella sampling, US) 和WHAM. 该方法的基本公式为下式，此处P 即为在某点处发现例子的概率

A?????1?lnP??A0

这个方法涉及到当某点势垒过高式，取样很少甚至取不到样，可以采用给粒子加一个bias force，使之可以到达高势垒区。具体的处理方法本文不介绍，ABF方法在US基础上发展而来, 在namd计算ABF时候, US 并不是必选选项。

下面给出的例子是ABF-Mar2008 中， 10个ALA 的多肽（10Ala）从alpha 螺旋到拉伸成单链过程中的自由能变化，体系在气相中共104 个原子非常适合(现在)研究各种自由能方法。在tutorial中给出的那个conf 并不能直接使用，这里我们只要注意下面的内容就可以了

source ../abf-1.8/abf.tcl abf coordinate distance abf abf1 10 abf abf2 92 abf xiMin 12.0 abf xiMax 32.0 abf dxi 0.1 abf fullSamples 200 abf dSmooth 0.1 abf forceConst 10.0 abf writeXiFreq 100 abf outFile atom-based_01.abf

命令详解

1. ABF

ABF (average biasing force) 是namd 中实现probability densities 的一种方法. 同这个方法相关

source ../abf-1.8/abf.tcl

在2.6 版中，abf 是通过Jerome Henin等人编写的tclForce 的脚本实现的，这里要source 这个脚本。具体问题依你的情况而定，这个文件在namd 安装目录下的lib/abf 文件夹中，即时你用windows 版本，也要用/ 而不是\\。

在namd2.7 中ABF方法已经用c++写入内核, 使用方法参考ug.

abf coordinate distance

在ABF 的命令中，要以abf 开头，实际上这是在abf.tcl 中定义的一个函数

proc abf { keyword value } { set ::ABF::keyword $keyword set ::ABF::value $value … … }

此处的coordinate指的是reaction coordinate，RC，有原子距离，基团距离，在z 方向的距离等等但都是distance 而没有角度等等

abf abf1 10 abf abf2 92

前面说的距离，就是这两个之间的距离

数字是AtomID (index)，如果是基团就要用{} 都写起来

注意vmd 和namd 对AtomID 的起点有所不同，还有tcl 语法中{}和” ” 的相同与不同之处

abf xiMin 12.0 abf xiMax 32.0 abf dxi 0.1

这是所要考察RC 的起始点，以及bins 的长度 如果体系的自由能profile较复杂，dxi 要小一些，但会增加计算量

在这里区别两个概念: 如果把长的RC分段计算，比如每5A一段，共4端，这样每一段叫做一个window. 上面的dxi 叫做bin.

abf fullSamples 200

在每个bin中200 步之后取样，用以消除nonequilibrium effects.

abf dSmooth 0.1

平滑处理，消除曲线的波动，但会使精细的profile变化消失

这个参数默认值为0

abf forceConst 10.0 力常数，默认为10.0

同自由能profile 以及dxi 的选取有关

abf writeXiFreq 100

在stdout 中每个100步输出类似下面的行 TCL: ABF> Xi at timestep 48200 : 18.7445838067

abf outFile atom-based_01.abf ABF 的输出文件，内容如下

# xi A(xi) av_force n_samples 10.100 0.0000 0.9421 1026 10.300 -0.1954 1.0121 1398 … … A(xi) 就是在这点处的A 值，av_force 是A 对xi的导数

2.FEP/TI

因为nonequilibrium work 内容较多，这一节先介绍FEP 和TI. 前面讲过，只是自由能的定义和处理方法不同，思想和和在namd 中的运行方式是相同的如果说ABF 是对位置进行限制，那么FEP 和TI 就是对状态进行限制. FEP/TI 计算的是某个原子或基团逐渐grow 和vanish 过程中自由能的变化，在user guide 中放在alchemical free energy中进行介绍的。

FEP/TI 可以计算蛋白质mutate 产生的自由能变，可以通过设计热力学循环计算某两个状态间的能差，也可以，举个例子，在水中的一个乙醇渐渐消失，在膜中另一个乙醇渐渐出现，计算乙醇在水和膜中的自由能差。

使用一个≤1的参数lambda λ来定义粒子存在的状态，0为初始态a，1为终态b. 体系哈密顿量用公式表示

Hx,px;???H0?x,px???Hb?x,px???1???Ha?x,px?FEP 的自由能表达如下式

?Aa?b??lnexp????Hb?x,px??Ha?x,px????1??TI 的用法就出现在ug下面，同FEP 的用法何

a其相似来尔

TI 的自由能表达如图3 中所示的形式

?A??1?H(x,px;?)0??

?这里只是让大家看一下公式长得什么样子，计算时候用不到处理这些。

同样使用ABF-Mar2008 中的例子，就是那个甲烷从真空中到水中的能变. 教程中设计了一个热力学循环，不细讲这个循环，总之就是那个物理过程的自由能变等于水中出现一个甲烷的alchemy 过程的自由能变。

下面是namd ug2.7 中的例子 fepOn On fepFile ion.fep fepCol B fepOutfile ion.fepout fepOutFreq 5 fepEquilSteps 5000 set Lambda0 0.0 set dLambda 0.1 while {$Lambda0 <= 1.0} { lambda $Lambda0 set Lambda0 \\ [expr \\$Lambda0 + \\$dLambda] lambda2 $Lambda0 run 10000 } fepFile ion.fep fepCol

B

指定哪些原子出现或消失，当为-1时原子逐渐消失，为1时逐渐出现，为0时不变

lambda $Lambda0

lambd2 [expr \\$Lambda0 + \\$dLambda]

每次计算一个window 内的能变，这个过程在一个循环中完成，相当于取了从0.0到1.0 的每个小区间，每次run 10000步，共10万步

thermInt On tiFile ion.alch.pdb tiCol B tiOutfile ion.ti.out tiOutFreq 5 tiEquilSteps 5000 tiLambda 0 run 10000 tiLambda 0.00001 run 10000 tiLambda 0.0001 run 10000 tiLambda 0.001 run 10000 tiLambda 0.01 run 10000 set Lambda 0.1 while {$Lambda <= 0.9} { tiLambda $Lambda run 10000 set Lambda [expr $Lambda + 0.1] } tiLambda 0.99 run 10000 tiLambda 0.999 run 10000 tiLambda 0.9999 run 10000 tiLambda 0.99999 run 10000 tiLambda 1 run 10000

TI 只需要定义一个tiLambda 即可，同FEP 有两个lambda 不同

两种方法都有一个问题就是当lambda 接近0或1的时候，会出现很大的斜率. 需要采用这种越来越精细的划分，对FEP 亦当如此. 在tutorial 中有一

个fep.tcl 文件，可以简化这种操作，作为support information 刊出, 用法大家可以参考tutorial 中的conf 文件自己捉摸.

计算蛋白质mutate 过程时候，需要制作dual 的psf 文件，这个在psfgen 中处理不了. 制作方法参考: Alchemical FEP Tutorial:

http://www.ks.uiuc.edu/Research/namd/tutorial/fep/ 这里不做介绍

如此缓慢的拖动速度。 这个技术的英文名称自身就极为混乱，很多人把所有这些受力模拟称为\。这个称谓主要来自美国某位比较知名教授，本人觉得莫名其妙，同时十分难翻译。已有的译法有受控分子动力学，拉伸分子动力学，牵引分子动力 学……。另外的一种称谓与AFM联系起来，分别称常速拉伸与常力拉伸为\MD simulation (FPMD)\和\。这种跟模拟原理相结合的称为比较合理，可以翻译为\力学探针分子动力学模拟\和\力学探钳分子动力学模拟\。也有人建议称之为\in silico\silico 是拉丁语， 意为\硅\，因此一般写成斜体,In silico是指\在硅之中\，也就是\进行于电脑中，或是经由电脑模拟\之意，衍生自另外两个在生物学上常用的术语：in vivo及in vitro。这个称谓尚未见中文译名。

3. SMD

3.1. 先介绍一点SMD 的基本概念

SMD其原理是在常规分子动力学方法中以类似于原子力显微镜技术的方式引入外力，使得分子体系沿着特定的方向演化，从而可以在较短模拟时间内模拟体系在某 一特定自由度上较大时间尺度的变化，故称之为受控分子动力学。与常规分子动力学相比，SMD模拟方法的主要优点是，它可以在目前分子动力学纳秒模拟所需的 计算时间内，引发分子体系较大的构象变化，从而模拟体系的动态变化过程。根据施加外力的不同方式，SMD模拟方法分为恒力和恒速SMD，他们分别对应于 AFM方法中测定断裂力和测定强制寿命的实验方案。需要强调的是，一般说来，SMD模拟的目的并不是精确计算施加的外力的大小，而是模拟体系在外力的引导 下所产生的变化。由于SMD模拟过程是外力引发的热力学不可逆过程，所以一般只能提供定性的信息。当拖动的速度足够小，使得整个过程可逆的时候，外力所做 的功即等于体系自由能的变化，但是目前计算机的计算能力还难以模拟

3.2. SMD 计算拉伸功

NAMD 中的SMD仅仅能实现cvSMD，cfSMD 可以通过namd 提供constantForce 实现，SMD 中要指定拉伸的方向和速度，以及拉伸时候的力常数。在stdout 中输出当前时刻的位置和拉伸点受力的大小。

当拉伸时候，往往需要固定或限制另一部分原子，否则整条蛋白质就会一起被拉跑了。这些操作参数较少，较为简单，不做介绍。

cfSMD 用处不大，Liu et al. J. Phys. Chem. B 110 (2006) 12789 中有一个精妙例子，通过cfSMD 拉伸在跨膜短杆菌肽A 中Na+ ，使Na+ 在通过这个通道时候在自由能的谷底持续停留而在峰处一下子就跳过去了，可以快速估计粒子在管中的自由能变，但不能精确得到profile。