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
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