umbrella sampling script

From: Hyundeok Song (songhk_at_email.uc.edu)
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(http://membrane.urmc.rochester.edu/Software/WHAM/WHAM.html) 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.

hyun.

- 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

# 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

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