VMD-L Mailing List
From: Sanjeev (snjv_bs_at_rediffmail.com)
Date: Sun Mar 18 2007 - 01:30:53 CDT
- Next message: John Stone: "Re: Scripting problem with h-bond listing"
- Previous message: Axel Kohlmeyer: "Re: question about gofr"
- Next in thread: John Stone: "Re: Scripting problem with h-bond listing"
- Reply: John Stone: "Re: Scripting problem with h-bond listing"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]
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
- Next message: John Stone: "Re: Scripting problem with h-bond listing"
- Previous message: Axel Kohlmeyer: "Re: question about gofr"
- Next in thread: John Stone: "Re: Scripting problem with h-bond listing"
- Reply: John Stone: "Re: Scripting problem with h-bond listing"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]