VMD-L Mailing List
From: Sanjeev (snjv_bs_at_rediffmail.com)
Date: Sun Mar 18 2007 - 14:16:42 CDT
- Next message: Axel Kohlmeyer: "Re: Re: Scripting problem with h-bond listing"
- Previous message: John Stone: "Re: Scripting problem with h-bond listing"
- Maybe in reply to: John Stone: "Re: Scripting problem with h-bond listing"
- Next in thread: Axel Kohlmeyer: "Re: Re: Scripting problem with h-bond listing"
- Reply: Axel Kohlmeyer: "Re: Re: Scripting problem with h-bond listing"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]
Dear John,
I have tried to fix it with '$sel update' and the likes, but it did not help. I have now nine frames and there are valid h-bonds in all of them. The only interesting aberration I noticed is that the atomselect_n is increasing by 40 each time I run, which I do not consider to be unusual either, except that the first frame has 13 bonds/13*3=39 sized array (for acceptors, donors and hydrogens). Beyond that it beats me. I suppose the aberration is reproducible as I have run it with different PDBs. The following is 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
- Next message: Axel Kohlmeyer: "Re: Re: Scripting problem with h-bond listing"
- Previous message: John Stone: "Re: Scripting problem with h-bond listing"
- Maybe in reply to: John Stone: "Re: Scripting problem with h-bond listing"
- Next in thread: Axel Kohlmeyer: "Re: Re: Scripting problem with h-bond listing"
- Reply: Axel Kohlmeyer: "Re: Re: Scripting problem with h-bond listing"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]