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 的每个小区间，每 …… 此处隐藏：4984字，全部文档内容请下载后查看。喜欢就下载吧 ……