umbrella sampling script

From: Hyundeok Song (
Date: Tue May 05 2009 - 17:30:28 CDT

Dear NAMD users

I calculated 1 dimensional Potential of Mean Force(PMF) of one ion using umbrella sampling.
Alan Grossfield's WHAM( was used to get unbiased distribution. I used same spring constant that was used in NAMD simulation when I used WHAM to calculate 1D PMF. (10 kcal/mol/A^2)

1D PMF that I got is around 2 times higher than something I expected. I should check the reason. I guess there will many possible reasons. Here is my tclforce script that was used at umbrella sampling.

I put one ion at 10A from CM of proteins.
Reaction coordinate is the z position of the ion from the center of mass of proteins.

I want to make sure if my script is correct.
If someone's already done umbrella sampling, could you send the script so that I can compare mine with that?
Any comment is also welcome. Thanks.


- harmonic.tcl -

set dk 10 ; # force constant, unit of kcal/mol/A^2
set ref 10 ;# the distance of the ion from the CM of proteins, where reaction coordinate is z axis.

set K1 5366
addatom $K1

# addgroup: request center of mass coordinates of proteins, here 20 is just for examples. More than 500 were used at simulation.

set groupid [ addgroup {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20} ]

proc innerp {v1 v2} {
    set d 0.0
    foreach c1 $v1 c2 $v2 { set d [expr {$d + $c1*$c2} ] }
    return $d

proc disHarm { a1 } {
     global dk groupid ref

     # load the coordinates of the ion
     loadcoords x

     set r1 $x($a1)
     set z1 [lindex $r1 2]

     # load the coordinates of the CM of proteins
     set r2 $x($groupid)
     set x2 [lindex $r2 0]
     set y2 [lindex $r2 1]
     set z2 [lindex $r2 2]

     # z22: reference point(z component) of the ion, 10A from CM
     set z22 [ expr {$z2 + $ref} ]

     # ref12: reference position of the ion.
     set ref12 {}
     lappend ref12 $x2 $y2 $z22

     set vect [vecsub $x($a1) $ref12 ]
     set d2 [ innerp $vect $vect ]
     set d [ expr {sqrt ($d2) }]

     puts "zpos:[expr $z1 - $z2]" ;# z position of the ion from C.M. of protein
     addenergy [ expr {0.5 * $dk * $d * $d }]
     addforce $a1 [vecscale -$dk $vect ]


proc calcforces {} {
     global K1

     disHarm $K1

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:52:44 CST