From: John Stone (johns_at_ks.uiuc.edu)
Date: Sun Mar 18 2007 - 01:54:15 CDT

Sanjeev,
  I would add some additional 'puts' calls to determine how many
bonds you're getting from the first half of the code that calls
'measure hbonds'. If you're getting zero, then you'll want to
delve closer into this and determine what you're actually
getting from 'measure hbonds', and whether or not you've set
the selection criteria too narrowly. Just doing a quick read
through your script I didn't see anything obvious, so I suggest
that you print more information as the script runs to determine
which part of the code is actually giving you trouble.

  John Stone
  vmd_at_ks.uiuc.edu

On Sun, Mar 18, 2007 at 06:30:53AM -0000, Sanjeev wrote:
> Hello,
> I have written a script to find all h-bonds within a cutoff criteria. It is to run on all the frames that are loaded. But strangely it runs only on the first frame but does not process (skips) the rest of frames. Could someone hint at the bug please? The following is the script that is run (after verifying that all (two here) the snapshots are loaded).
>
> ------ > hb.tcl
>
> set numframes [molinfo top get numframes]
> set sel [atomselect top all]
>
> set PI 3.14159
>
> for {set i 0} {$i < $numframes} {incr i} {
> $sel frame $i
>
> puts "Frame : $i"
>
> foreach {d a h} [measure hbonds 3.2 30 $sel] {
> set numbonds 0
> foreach ctr $d {
> set d_ary($numbonds) $ctr
> set numbonds [incr numbonds]
> }
>
> set numbonds 0
> foreach ctr $h {
> set h_ary($numbonds) $ctr
> set numbonds [incr numbonds]
> }
>
> set numbonds 0
> foreach ctr $a {
> set a_ary($numbonds) $ctr
> set numbonds [incr numbonds]
> }
> }
>
> for {set j 0} {$j < $numbonds} {incr j} {
> set sel [atomselect top "index $d_ary($j)"]
> set d_a [$sel get name]
> set d_r [$sel get resname]
> set d_d [$sel get {x y z}]
> set d_i [$sel get resid]
> set D [lindex [$sel get { x y z }] $i]
>
>
> set sel [atomselect top "index $a_ary($j)"]
> set a_a [$sel get name]
> set a_r [$sel get resname]
> set a_d [$sel get {x y z}]
> set a_i [$sel get resid]
> set A [lindex [$sel get { x y z }] $i]
>
>
> set sel [atomselect top "index $h_ary($j)"]
> set h_a [$sel get name]
> set h_r [$sel get resname]
> set h_d [$sel get {x y z}]
> set h_i [$sel get resid]
> set H [lindex [$sel get { x y z }] $i]
>
>
> set dist [vecdist $D $A]
> set angl [expr 180 - acos([vecdot [vecsub $D $H] [vecsub $H $A]]/([vecdist $D $H] * [vecdist $H $A])) * 180 / $PI]
>
> set out_1 "$j $d_a $d_r $d_i $h_a $h_r $h_i $a_a $a_r $a_i"
> set out_2 [format "%s %5.2f %5.1f" $out_1 $dist $angl]
>
> puts $out_2
> }
>
> }
>
>
> -----< end file
>
>
> The o/p is:
>
> ------> output
>
> vmd > source hb.tcl
> No. of frames: 2
> Frame : 0
> 0 OH TYR 92 HH TYR 92 O LYS 37 2.70 157.8
> 0 N ARG 39 H ARG 39 O THR 36 3.79 125.8
> 0 N LYS 91 H LYS 91 O ASN 94 2.81 171.9
> 0 N GLY 68 H GLY 68 O CYX 65 2.88 147.5
> 0 NZ LYS 37 HZ3 LYS 37 OD2 ASP 38 3.62 123.1
> 0 ND2 ASN 67 HD21 ASN 67 OE1 GLN 69 2.88 162.0
> 0 NE2 GLN 69 HE22 GLN 69 OD1 ASN 71 3.23 168.0
> 0 N GLN 69 H GLN 69 OD1 ASN 67 3.93 161.4
> 0 OG SER 90 HG SER 90 OE1 GLU 86 2.71 122.7
> 0 N ALA 96 H ALA 96 OG SER 90 3.17 152.8
> 0 N SER 90 H SER 90 O THR 87 3.06 156.7
> 0 N LYS 66 H LYS 66 CG ASP 121 3.77 169.9
> 0 N LYS 66 H LYS 66 OD2 ASP 121 3.00 163.4
> 0 N LYS 37 H LYS 37 O ASN 34 3.11 156.3
> 0 N THR 36 H THR 36 O MET 30 3.54 126.8
> 0 ND2 ASN 71 HD22 ASN 71 O CYX 110 2.82 153.3
> 0 ND1 HIP 119 HD1 HIP 119 O1P RG3 126 2.77 93.8
> 0 O WAT 127 H2 WAT 127 O1P RG3 126 2.83 92.9
> 0 NE2 HIE 12 HE2 HIE 12 O WAT 127 3.12 109.1
> 0 CA HIP 119 HA HIP 119 O WAT 127 2.84 142.4
> 0 N PHE 120 H PHE 120 O WAT 127 2.35 141.1
> 0 NZ LYS 41 HZ3 LYS 41 OD1 ASN 44 3.06 139.3
> 0 NZ LYS 41 HZ1 LYS 41 O VAL 43 2.79 114.9
> 0 OH TYR 97 HH TYR 97 O LYS 41 2.73 136.4
> 0 NE2 HIP 119 HE2 HIP 119 OD1 ASP 121 2.81 162.3
> 0 OG1 THR 87 HG1 THR 87 O ALA 96 3.59 54.5
> 0 N CYX 65 H CYX 65 O GLN 69 2.90 166.4
> 0 NZ LYS 66 HZ3 LYS 66 O ASP 121 2.74 160.6
> 0 N GLU 86 H GLU 86 O PRO 42 2.89 143.3
> 0 N LYS 98 H LYS 98 O ARG 85 2.80 148.0
> 0 OG SER 80 HG SER 80 O SER 18 2.82 149.1
> 0 ND2 ASN 103 HD22 ASN 103 O THR 78 3.04 166.3
> Frame : 1
> vmd >
>
> ------< end output
>
>
> Thanks in advance,
> -Sanjeev

-- 
NIH Resource for Macromolecular Modeling and Bioinformatics
Beckman Institute for Advanced Science and Technology
University of Illinois, 405 N. Mathews Ave, Urbana, IL 61801
Email: johns_at_ks.uiuc.edu                 Phone: 217-244-3349
  WWW: http://www.ks.uiuc.edu/~johns/      Fax: 217-244-6078