NAMD计算自由能(2)

给出了一个等式，能精确描述自由能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 的文件中