From: Jian Dai (jdai2_at_fsu.edu)
Date: Mon Dec 12 2011 - 10:13:16 CST
Dear all:
I'm trying to restrain a bond length between 2 atoms and a dihedral angle among 4 atoms. I expect tclforces and extrabonds will give close results (or not?) though fluctuation is expected.
By using extrabond feature I get expected results. But when using tclforces the result is not what is expected. Can anybody shoot me some suggestions please?
Another question is: from the tclforces example, the potential form seems to be V = 1/2*k*(x-x0)^2 for dihedral, but in the source code of NAMD (also in the manual), for example, in the function computeForce of ComputeBonds.C, the potential seems to be in V = k*(x-x0)^2. Which form should I use in the tclforces? Using the 1st form, the force will be -k*dx; while using the 2nd form the force will be -2*k*dx. Can anyone clear my confusion? Thanks.
The tclforce script is as follows.
# Definition for vecdot and veccross, copied from VMD source code
set DEG2RAD 0.01745
set k_dist 5
set init_dist 1.8
set dist_atoms {}
lappend dist_atoms [atomid A 37 HE2]
lappend dist_atoms [atomid B 37 ND1]
foreach atom [lsort -unique $dist_atoms] {
addatom $atom
}
set k_dihed 5
set init_dihed -20
set dihed_atoms {}
set A_CD2 [atomid A 37 CD2]
set A_NE2 [atomid A 37 NE2]
set B_CG [atomid B 37 CG]
set B_ND1 [atomid B 37 ND1]
lappend dihed_atoms $A_CD2
lappend dihed_atoms $A_NE2
lappend dihed_atoms $B_ND1
lappend dihed_atoms $B_CG
foreach atom [lsort -unique $dihed_atoms] {
addatom $atom
}
proc calcforces {} {
global DEG2RAD
global k_dist init_dist
global dist_atoms
global k_dihed init_dihed
global dihed_atoms
loadcoords crds
# Distance restraint
# Refer to NAMD source code ComputeBonded.C : compute_bonds()
foreach {atom1 atom2} $dist_atoms {
set r12 [vecsub $crds($atom2) $crds($atom1)]
set r [veclength $r12]
set diff [expr {$r - $init_dist}]
set f12 [vecscale $r12 [expr {$diff * $k_dist/$r}]]
addenergy [expr {0.5*$k_dist*$diff*$diff}]
addforce $atom1 $f12
addforce $atom2 [vecscale -1. $f12]
}
# Dihedral restraint. All angle unit is converted into radian
# Refer to NAMD tclforces example
foreach {atom1 atom2 atom3 atom4} $dihed_atoms {
set init_dihed [expr {$init_dihed*$DEG2RAD}]
set dihed [getdihedral $crds($atom1) $crds($atom2) $crds($atom3) $crds($atom4)]
set dihed [expr {$dihed*$DEG2RAD}]
set d_dihed [expr {$dihed - $init_dihed}]
addenergy [expr {0.5*$k_dihed*$d_dihed*$d_dihed}]
set df [expr {-1.*$k_dihed*$d_dihed}]
foreach {g1 g2 g3 g4} [dihedralgrad $crds($atom1) $crds($atom2) $crds($atom3) $crds($atom4)] {}
addforce $atom1 [vecscale $g1 $df]
addforce $atom2 [vecscale $g2 $df]
addforce $atom3 [vecscale $g3 $df]
addforce $atom4 [vecscale $g4 $df]
}
}
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:58:00 CST