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