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 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 : Mon Dec 31 2012 - 23:21:03 CST