# Re: tclforces and extrabonds

From: Aron Broom (broomsday_at_gmail.com)
Date: Mon Dec 12 2011 - 11:30:36 CST

Perhaps a pointless suggestion, but have you tried using collective
variables (colvars) for this? The harmonic function seems like it would be
ideal.

~Aron

On Mon, Dec 12, 2011 at 11:13 AM, Jian Dai <jdai2_at_fsu.edu> wrote:

> 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