tclforces and extrabonds

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 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] {
}

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] {
}

proc calcforces {} {
global k_dist init_dist
global dist_atoms
global k_dihed init_dihed
global dihed_atoms

# 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}]]
}

# Dihedral restraint. All angle unit is converted into radian
# Refer to NAMD tclforces example
foreach {atom1 atom2 atom3 atom4} \$dihed_atoms {
set dihed [getdihedral \$crds(\$atom1) \$crds(\$atom2) \$crds(\$atom3) \$crds(\$atom4)]