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。

在10Ala 中采用了tclBC的方式实现cvSMD，为实现拉伸10Ala 一端的N，并限制 (constraint )另一端的N原子，实现拉伸。程序每隔一段时间往一个文本中输出当前时刻，针尖位置，和作用在针尖上力等信息。拉伸方向为z 的正方向，但肽链上两个施力点并不平行于z，拉伸时候会拉过去，不需要调整。

SI中有一个smd.tcl 文件，在conf 文件中放到tclForcesScript 的后面. 这个代码很容易读懂，建议大家读一下，有必要时可以更改里面的部分内容 tclBC较tclForce 有很多优点，此处用tclBC 可能更好一些.

教程中使用了3种拉伸速度，拉伸20A 总耗时分别为20ps，2ns 和200ns，第三种为可逆拉伸，第二种近似可逆拉伸，用于自由能计算，第一种为不可逆拉伸。需要再修改outputFreq，保证既能几下足够多的数据，又不使文件太大

当拉伸速度足够慢时候，拉伸力所做的功就是体系的自由能变, 功的积分公式如下

这个脚本前已定义了拉伸速度v (A/dt),积分时间dt (存储两组数据的间隔)，以及各个存储时刻的针尖受力。 set w {} set fsum 0 foreach ftemp $f { set fsum \\ [expr $fsum + $ftemp * $v * $dt] lappend w $fsum }

10Ala教程给出的各个tcl 脚本都比较简单，但需要多个脚本合作才能完成一项具体操作。通过运行这个脚本，可以计算单次模拟过程中，针尖所做的功，根据热力学第二定律，W ≥ ΔA 仅当完全可逆拉伸时候等号成立. 这一步详细脚本从略，但需要注意v，dt 等参数，这个v 同smd.tcl 定义的v 都是拉伸速度，但并不一样

下面重磅推出本文的明星方程，Jarzynski 方程

W(t)?v?f(t')dt'

0t'在算法实现时候，可以把积分化为求和。 下面是calcwork1.tcl 文件中的内容，程序在source

热II定律给出了拉伸功同自由能之间的不等式，

仅仅能用来估计自由能变的上限而Jarzynski 方程

给出了一个等式，能精确描述自由能profile, 也就是PMF. S. Park, JCP, 119 (2003) 3559 中提出了J方程的累级近似，并发现用J 方程的二级近似比原方程更精确些。10Ala 例子中就同时计算了J 方程，一级和二级近似得到的free energy profile，并同可逆拉伸功做了比较。

教程中，以2 ns 的时间拉伸10Ala 20A 的距离，平行10个样(这里的<>才算真正意义上的系宗平均)，得到的位置和力的文件存储在一个文件中。本文单独存储10个文件，这都是由smd.tcl 生成的。

首先说一下各个量的单位，在这里的smd.tcl 中，

set Fexp {} Set F1 {} set F2 {} set Tclfreq set v set T set dt set v

50 ;# every 0.1ps ;# A/timestep

0.00002

在计算pmf 的pmf.tcl 文件中

0.6 0.1 0.01

;# kT, T=300k, Kcal/mol ;# 存储间隔，ps ;# 0.01A/dt

当使用calcwork2.tcl (功能同calcwork1.tcl 几乎一样，只是能计算多条轨迹) 计算出来各条轨迹的w 后，就可以通过J 方程计算自由能了。那个方程也不深奥，也就仅仅是一个对10 条轨迹的平均问题， 详细算法么，下面的代码就像流水账一样，比文字还好懂

for {set i 0} {$i < 20001} {incr i 1} { set texp 0 set t1 0 set t2 0 foreach l [array names w] { set e [lindex $w($l) $i] set texp [expr $texp + exp([expr - $e / $T]) ] set t1 [expr $t1 + $e] set t2 [expr $t2 + $e * $e] } lappend Fexp [expr - $T * log([expr $texp / 10])] lappend F1 [expr $t1 / 10] lappend F2 [expr $t1 / 10 - $t2 / 12 + $t1 * $t1 / 120 ] }

后面几行中有除以10 的代码，这个10 是指10条轨迹，那个12 是什么呢在编程上这种来历不明的数叫做幻数，是一种极坏的编程风格但这里这么写是提醒大家注意，10Ala 教程中给出的是10 和100 ，而不是12 和120认为这里写错了，给个提示，那个T=0.6，大家想一下，如果我错了还请指正。

下面的PMF 结果是我求得的，发现同S. Park 的结论有点差别，那位同学说二级近似更精确些，我

的结果好像是不近似更精确。毕竟人家Jarzynski 明明推出来一个等式。

附上一些代码, 部分内容同10Ala教程中名同菜不同 calcwork2.tcl set w($l) {} set fsum 0 foreach ftemp $f($l) { set fsum [expr $fsum + $ftemp * $v * $dt] lappend w($l) $fsum } } cumulants.tcl set Fexp {} set F1 {} set F2 {} for {set i 0} {$i < 20001} {incr i 1} { set texp 0 set t1 0 set t2 0 foreach l [array names w] { set e [lindex $w($l) $i] set texp [expr $texp + exp([expr - $e / $T]) ] set t1 [expr $t1 + $e] set t2 [expr $t2 + $e * $e] } lappend Fexp [expr - $T * log([expr $texp / 10])] lappend F1 [expr $t1 / 10] lappend F2 [expr $t1 / 10 - $t2 / 12 + $t1 * $t1 / 120 ] } set pmff [open pmf.out w] puts $pmff \ Fexp F1 F2\set i 0 while {$i < 20001} { set ci [expr 13 + $v * $i /10] set fei [lindex $Fexp $i] set f1i [lindex $F1 $i] set f2i [lindex $F2 $i] puts $pmff \ incr i } pmf.tcl

相当于c语言中的main

计算PMF

计算各条轨迹的功

foreach l [array names f] { source load-traj.tcl set T 0.6 set dt 0.1 set v 0.01 source calcwork2.tcl source cumulants.tcl exit 使用方法，

vmd -dispdev text -e pmf.tcl 结果在一个较pmf.out 的文件中

本文来源：https://www.bwwdw.com/article/18yp.html