VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Sun Mar 18 2007 - 15:03:09 CDT

On 18 Mar 2007, Sanjeev wrote:

SB> Dear John, I have tried to fix it with '$sel update' and the likes,
SB> but it did not help. I have now nine frames and there are valid
SB> h-bonds in all of them. The only interesting aberration I noticed is
SB> that the atomselect_n is increasing by 40 each time I run, which I
SB> do not consider to be unusual either, except that the first frame
SB> has 13 bonds/13*3=39 sized array (for acceptors, donors and
SB> hydrogens). Beyond that it beats me. I suppose the aberration is

please re-read your script and watch the output. you're constantly
overwriting $sel, no wonder it does not work. you have to use other
variables names in the innerloops.

cheers,
   axel.

SB> reproducible as I have run it with different PDBs. The following is
SB> the script as of now:

--- hb.tcl ---

set sel [atomselect top all]
set n_frames [molinfo top get numframes]

set PI 3.14159

puts "No. of frames: $n_frames"

for {set i 0} {$i < $n_frames} {incr i} {
   puts "Processing frame: $i"

   $sel frame $i
   $sel update

   foreach {d a h} [measure hbonds 2.8 30 $sel] {
      puts "Sel is: $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]
      }

      puts "Numbonds : $numbonds"
   }

   for {set j 0} {$j < $numbonds} {incr j} {
      set sel [atomselect top "index $d_ary($j)" frame $i]
      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)" frame $i]
      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)" frame $i]
      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 "$i $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
   }

}

unset sel
unset d_ary a_ary h_ary
unset n_frames

--- hb.tcl ---

The output is below:

--- output ---

vmd > source hb.tcl
No. of frames: 9
Processing frame: 0
Sel is: atomselect0
Numbonds : 13
0 OG1 THR 67 HG1 THR 67 OD1 ASN 65 2.79 162.0
0 O WAT 4626 H2 WAT 4626 O1P RA3 135 2.67 176.1
0 OH TYR 98 HH TYR 98 O LYS 38 2.73 171.6
0 N ILE 86 H ILE 86 O ARG 97 2.76 169.0
0 NH1 ARG 1 1HH1 ARG 1 O PHE 5 2.73 166.4
0 NE2 GLN 14 2HE2 GLN 14 O TRP 10 2.75 162.9
0 NZ LYS 38 HZ3 LYS 38 OD1 ASN 41 2.80 163.6
0 N ALA 99 H ALA 99 O ASP 84 2.76 153.4
0 N TRP 10 H TRP 10 O THR 6 2.79 164.8
0 N CYX 23 H CYX 23 OD2 ASP 100 2.80 165.1
0 OG SER 74 HG SER 74 O PHE 76 2.70 176.8
0 NH1 ARG 77 1HH1 ARG 77 O ARG 104 2.77 157.6
0 N THR 46 H THR 46 O LEU 44 2.76 152.5
Processing frame: 1
Sel is: atomselect39
Numbonds : 0
Processing frame: 2
Sel is: atomselect39
Numbonds : 0
Processing frame: 3
Sel is: atomselect39
Numbonds : 0
Processing frame: 4
Sel is: atomselect39
Numbonds : 0
Processing frame: 5
Sel is: atomselect39
Numbonds : 0
Processing frame: 6
Sel is: atomselect39
Numbonds : 0
Processing frame: 7
Sel is: atomselect39
Numbonds : 0
Processing frame: 8
Sel is: atomselect39
Numbonds : 0
vmd > source hb.tcl
No. of frames: 9
Processing frame: 0
Sel is: atomselect40
Numbonds : 13
0 OG1 THR 67 HG1 THR 67 OD1 ASN 65 2.79 162.0
0 O WAT 4626 H2 WAT 4626 O1P RA3 135 2.67 176.1
0 OH TYR 98 HH TYR 98 O LYS 38 2.73 171.6
0 N ILE 86 H ILE 86 O ARG 97 2.76 169.0
0 NH1 ARG 1 1HH1 ARG 1 O PHE 5 2.73 166.4
0 NE2 GLN 14 2HE2 GLN 14 O TRP 10 2.75 162.9
0 NZ LYS 38 HZ3 LYS 38 OD1 ASN 41 2.80 163.6
0 N ALA 99 H ALA 99 O ASP 84 2.76 153.4
0 N TRP 10 H TRP 10 O THR 6 2.79 164.8
0 N CYX 23 H CYX 23 OD2 ASP 100 2.80 165.1
0 OG SER 74 HG SER 74 O PHE 76 2.70 176.8
0 NH1 ARG 77 1HH1 ARG 77 O ARG 104 2.77 157.6
0 N THR 46 H THR 46 O LEU 44 2.76 152.5
Processing frame: 1
Sel is: atomselect79
Numbonds : 0
Processing frame: 2
Sel is: atomselect79
Numbonds : 0
Processing frame: 3
Sel is: atomselect79
Numbonds : 0
Processing frame: 4
Sel is: atomselect79
Numbonds : 0
Processing frame: 5
Sel is: atomselect79
Numbonds : 0
Processing frame: 6
Sel is: atomselect79
Numbonds : 0
Processing frame: 7
Sel is: atomselect79
Numbonds : 0
Processing frame: 8
Sel is: atomselect79
Numbonds : 0
vmd > 

---
It baffles me as to where the problem is. Especially now that I have use update frame as well.
-Sanjeev
-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.