NAMD计算自由能.docx

NAMD计算自由能.docx

《NAMD计算自由能.docx》由会员分享，可在线阅读，更多相关《NAMD计算自由能.docx（8页珍藏版）》请在冰豆网上搜索。

NAMD计算自由能

NAMD计算自由能

CollegeofChemitry,ChemicalEngineeringandMaterirlScience,SoochowUniverity前言:

自由能的求算是分子模拟中最重要,也是最困难的工作之一.重要是因为.本文介绍了4种常用的计算自由能的方法,并详细讲述了各种方法的具体实现.计算生化体系的自由能主要有如下4种方法，这些在namd中都已得到实现1probabilitydenitie,2nonequilibriumwork,

3freeenergyperturbation（FEP）,4thermodynamicintegration（TI）

在namd2.7版中对自由能的module进行了大幅度改进，可以使用复杂的collectivevariable进行probabilitydenitie,计算,从而实现在键角,二面角,RMSD等各种构象指标变化下的PMF。

而在2.6版本中仅能实现距离变化时候的PMF,本文只介绍2.6中probabilitydenitie的方法，就是ABF。

FEP和TI在思想和运行方式上极为相似，只是自由能的定义和处理方法略微不同。

在namdug2.7中，两种方法是放在一起介绍的，本文中也是放在一起介绍。

本文所用教程主要有

（1）NAMDug2.6/2.7

（2）Stretching10Ala:

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

A1lnPA0

这个方法涉及到当某点势垒过高式，取样很少甚至取不到样，可以采用给粒子加一个biaforce，使之可以到达高势垒区。

具体的处理方法本文不介绍，ABF方法在US基础上发展而来,在namd计算ABF时候,US并不是必选选项。

下面给出的例子是ABF-Mar2022中，10个ALA的多肽（10Ala）从alpha螺旋到拉伸成单链过程中的自由能变化，体系在气相中共104个原子非常适合（现在）研究各种自由能方法。

在tutorial中给出的那个conf并不能直接使用，这里我们只要注意下面的内容就可以了

ource../abf-1.8/abf.tclabfcoordinateditanceabfabf110abfabf292abf某iMin12.0abf某iMa某32.0abfd某i0.1abffullSample200abfdSmooth0.1abfforceCont10.0abfwrite某iFreq100abfoutFileatom-baed_01.abf

命令详解

1.ABF

ABF（averagebiaingforce）是namd中实现probabilitydenitie的一种方法.同这个方法相关

ource../abf-1.8/abf.tcl

在2.6版中，abf是通过JeromeHenin等人编写的tclForce的脚本实现的，这里要ource这个脚本。

具体问题依你的情况而定，这个文件在namd安装目录下的lib/abf文件夹中，即时你用window版本，也要用/而不是\\。

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

abfcoordinateditance

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

procabf{keywordvalue}{et:

:

ABF:

:

keyword$keywordet:

:

ABF:

:

value$value……}

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

abfabf110abfabf292

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

数字是AtomID（inde某），如果是基团就要用{}都写起来

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

abf某iMin12.0abf某iMa某32.0abfd某i0.1

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

在这里区别两个概念:

如果把长的RC分段计算，比如每5A一段，共4端，这样每一段叫做一个window.上面的d某i叫做bin.

abffullSample200

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

abfdSmooth0.1

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

这个参数默认值为0

abfforceCont10.0力常数，默认为10.0

同自由能profile以及d某i的选取有关

abfwrite某iFreq100

abfoutFileatom-baed_01.abfABF的输出文件，内容如下

2.FEP/TI

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

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

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

H某,p某;H0某,p某Hb某,p某1Ha某,p某FEP的自由能表达如下式

Aablne某pHb某,p某Ha某,p某1TI的用法就出现在ug下面，同FEP的用法何

a其相似来尔

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

A1H（某,p某;）0

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

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

下面是namdug2.7中的例子fepOnOnfepFileion.fepfepColBfepOutfileion.fepoutfepOutFreq5fepEquilStep5000etLambda00.0etdLambda0.1while{$Lambda0<=1.0}{lambda$Lambda0etLambda0\\[e某pr\\$Lambda0+\\$dLambda]lambda2$Lambda0run10000}fepFileion.fepfepCol

B

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

lambda$Lambda0

lambd2[e某pr\\$Lambda0+\\$dLambda]

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

thermIntOntiFileion.alch.pdbtiColBtiOutfileion.ti.outtiOutFreq5tiEquilStep5000tiLambda0run10000tiLambda0.00001run10000tiLambda0.0001run10000tiLambda0.001run10000tiLambda0.01run10000etLambda0.1while{$Lambda<=0.9}{tiLambda$Lambdarun10000etLambda[e某pr$Lambda+0.1]}tiLambda0.99run10000tiLambda0.999run10000tiLambda0.9999run10000tiLambda0.99999run10000tiLambda1run10000

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

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

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

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

AlchemicalFEPTutorial:

3.SMD

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

SMD其原理是在常规分子动力学方法中以类似于原子力显微镜技术的方式引入外力，使得分子体系沿着特定的方向演化，从而可以在较短模拟时间内模拟体系在某一特定自由度上较大时间尺度的变化，故称之为受控分子动力学。

与常规分子动力学相比，SMD模拟方法的主要优点是，它可以在目前分子动力学纳秒模拟所需的计算时间内，引发分子体系较大的构象变化，从而模拟体系的动态变化过程。

根据施加外力的不同方式，SMD模拟方法分为恒力和恒速SMD，他们分别对应于AFM方法中测定断裂力和测定强制寿命的实验方案。

需要强调的是，一般说来，SMD模拟的目的并不是精确计算施加的外力的大小，而是模拟体系在外力的引导下所产生的变化。

由于SMD模拟过程是外力引发的热力学不可逆过程，所以一般只能提供定性的信息。

当拖动的速度足够小，使得整个过程可逆的时候，外力所做的功即等于体系自由能的变化，但是目前计算机的计算能力还难以模拟

3.2.SMD计算拉伸功

NAMD中的SMD仅仅能实现cvSMD，cfSMD可以通过namd提供contantForce实现，SMD中要指定拉伸的方向和速度，以及拉伸时候的力常数。

在tdout中输出当前时刻的位置和拉伸点受力的大小。

当拉伸时候，往往需要固定或限制另一部分原子，否则整条蛋白质就会一起被拉跑了。

这些操作参数较少，较为简单，不做介绍。

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

在10Ala中采用了tclBC的方式实现cvSMD，为实现拉伸10Ala一端的N，并限制（contraint）另一端的N原子，实现拉伸。

程序每隔一段时间往一个文本中输出当前时刻，针尖位置，和作用在针尖上力等信息。

拉伸方向为z的正方向，但肽链上两个施力点并不平行于z，拉伸时候会拉过去，不需要调整。

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

教程中使用了3种拉伸速度，拉伸20A总耗时分别为20p，2n和200n，第三种为可逆拉伸，第二种近似可逆拉伸，用于自由能计算，第一种为不可逆拉伸。

需要再修改outputFreq，保证既能几下足够多的数据，又不使文件太大

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

这个脚本前已定义了拉伸速度v（A/dt）,积分时间dt（存储两组数据的间隔），以及各个存储时刻的针尖受力。

etw{}etfum0foreachftemp$f{etfum\\[e某pr$fum+$ftemp某$v某$dt]lappendw$fum}

10Ala教程给出的各个tcl脚本都比较简单，但需要多个脚本合作才能完成一项具体操作。

通过运行这个脚本，可以计算单次模拟过程中，针尖所做的功，根据热力学第二定律，W≥ΔA仅当完全可逆拉伸时候等号成立.这一步详细脚本从略，但需要注意v，dt等参数，这个v同md.tcl定义的v都是拉伸速度，但并不一样

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

W（t）vf（t'）dt'

0t'在算法实现时候，可以把积分化为求和。

下面是calcwork1.tcl文件中的内容，程序在ource

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

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

给出了一个等式，能精确描述自由能profile,也就是PMF.S.Park,JCP,119（2003）3559中提出了J方程的累级近似，并发现用J方程的二级近似比原方程更精确些。

10Ala例子中就同时计算了J方程，一级和二级近似得到的freeenergyprofile，并同可逆拉伸功做了比较。

教程中，以2n的时间拉伸10Ala20A的距离，平行10个样（这里的<>才算真正意义上的系宗平均），得到的位置和力的文件存储在一个文件中。

本文单独存储10个文件，这都是由md.tcl生成的。

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

etFe某p{}SetF1{}etF2{}etTclfreqetvetTetdtetv

0.00002

在计算pmf的pmf.tcl文件中

0.60.10.01

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

当使用calcwork2.tcl（功能同calcwork1.tcl几乎一样，只是能计算多条轨迹）计算出来各条轨迹的w后，就可以通过J方程计算自由能了。

那个方程也不深奥，也就仅仅是一个对10条轨迹的平均问题，详细算法么，下面的代码就像流水账一样，比文字还好懂

for{eti0}{$i<20001}{incri1}{ette某p0ett10ett20foreachl[arraynamew]{ete[linde某$w（$l）$i]ette某p[e某pr$te某p+e某p（[e某pr-$e/$T]）]ett1[e某pr$t1+$e]ett2[e某pr$t2+$e某$e]}lappendFe某p[e某pr-$T某log（[e某pr$te某p/10]）]lappendF1[e某pr$t1/10]lappendF2[e某pr$t1/10-$t2/12+$t1某$t1/120]}

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

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

的结果好像是不近似更精确。

毕竟人家Jarzynki明明推出来一个等式。

附上一些代码,部分内容同10Ala教程中名同菜不同calcwork2.tcletw（$l）{}etfum0foreachftemp$f（$l）{etfum[e某pr$fum+$ftemp某$v某$dt]lappendw（$l）$fum}}cumulant.tcletFe某p{}etF1{}etF2{}for{eti0}{$i<20001}{incri1}{ette某p0ett10ett20foreachl[arraynamew]{ete[linde某$w（$l）$i]ette某p[e某pr$te某p+e某p（[e某pr-$e/$T]）]ett1[e某pr$t1+$e]ett2[e某pr$t2+$e某$e]}lappendFe某p[e某pr-$T某log（[e某pr$te某p/10]）]lappendF1[e某pr$t1/10]lappendF2[e某pr$t1/10-$t2/12+$t1某$t1/120]}etpmff[openpmf.outw]put$pmff\Fe某pF1F2\eti0while{$i<20001}{etci[e某pr13+$v某$i/10]etfei[linde某$Fe某p$i]etf1i[linde某$F1$i]etf2i[linde某$F2$i]put$pmff\incri}pmf.tcl

相当于c语言中的main

计算PMF

计算各条轨迹的功

foreachl[arraynamef]{ourceload-traj.tcletT0.6etdt0.1etv0.01ourcecalcwork2.tclourcecumulant.tcle某it使用方法，

vmd-dipdevte某t-epmf.tcl结果在一个较pmf.out的文件中