VMD-L Mailing List

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Fri Nov 21 2008 - 13:22:04 CST

Hi Irene,
I'm not sure offhand why this doesn't work, but... could you verify that
you're using the double quote character " (ascii 34) to open and close
your atomselection text? This may just be a font issue, but your
addition prints differently, and I'm wondering if this is due to the
mail clients only or shows up in your raw tcl.
Also, I'd note that the use of anonymous atom selections (those which
are never assigned names) should be considered harmful, as it leaks
memory. This script ought in principle to be rewritten to instead make
one atom selection for each set of atoms of interest, do the puts
statement using those atomselections, and then delete the atom
selections. I know this isn't your fault, but figured I'd mention it
since I've seen this script pop up before.
Best,
Peter

Irene Newhouse wrote:
>
> I tried posting this 2 days ago, but never got an echo of it, so I
> don't think it made it.
>
>
>
> I've been using the hydrogen bond occupancy script sent by Dong Luo [
> http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l/573.57.html]
>
> I'm working with a glycan ligand, though, and there's lots & lots of
> choices for the
> ligand atoms involved in the H-bond. As the original script was
> intended for use
> with proteins, where you know which atom is involved, given a residue,
> the original
> script identifies only the resname & resid in the output.
>
> The original statement that does this reads:
>
> puts $outfile [format "%s%i \t %s%i \t\t %.2f%%" [[atomselect $mol
> "index $donor"] get resname] [[atomselect $mol "index $donor"] get
> resid] [[atomselect $mol "index $acceptor"] get resname] [[atomselect
> $mol "index $acceptor"] get resid] $occupancy]
>
> I tried to use the same pattern to add the atom name to the printed
> string, eg [[atomselect $mol “index $acceptor”] get name]:
>
> puts $outfile [format "%s%i%s \t %s%i%s \t\t %.2f%%" [[atomselect $mol
> "index $donor"] get resname] [[atomselect $mol "index $donor"] get
> resid] [[atomselect $mol “index $donor”] get name] [[atomselect $mol
> "index $acceptor"] get resname] [[atomselect $mol "index $acceptor"]
> get resid] [[atomselect $mol “index $acceptor”] get name] $occupancy]
>
> I also added another %s to the end of the format specifiers for the
> donor & acceptor atoms. But when I did this, I get a message giving me
> references to the usage of atomselect.
>
> Can anyone please help me out with this? I append the entire script,
> and how I'm calling it, below
>
> Thanks,
> Irene Newhouse
>
>
> proc hbond_occupancy { {dis 3.0} {ang 61} {sel1 protein} {sel2 none}
> {outfile stdout} {mol top} } {
>
> # found on VMD mailer; may not work
> http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l/573.57.html
>
> if {! [string compare $dis help]} {
>
> puts "Usage: hbond_occupancy distance angle selection1 selection2
> mol"
>
> return
>
> }
>
>
>
> if {! [string compare $mol top]} {
>
> set mol [molinfo top]
>
> }
>
>
>
> if {[string compare $outfile stdout]} {
>
> set outfile [open $outfile w];
>
> }
>
>
>
> set hbondallframes {}
>
> set hbondcount {}
>
>
>
> if { ! [string compare $sel1 protein]} {
>
> set sel1 [atomselect $mol $sel1]
>
> }
>
>
>
> set framenumberbackup [molinfo $mol get frame]
>
> set numberofframes [molinfo $mol get numframes]
>
>
>
> for { set i 0 } { $i < $numberofframes } { incr i } {
>
> molinfo $mol set frame $i
>
>
>
> if {[string compare $sel2 none]} {
>
> set hbondsingleframe [measure hbonds $dis $ang $sel1 $sel2]
>
> } else {
>
> set hbondsingleframe [measure hbonds $dis $ang $sel1]
>
> }
>
>
>
> for { set j 0 } { $j < [llength [lindex $hbondsingleframe 0] ] }
> { incr j } {
>
> set newhbond {}
>
> lappend newhbond [lindex $hbondsingleframe 0 $j ]
>
> lappend newhbond [lindex $hbondsingleframe 1 $j ]
>
> lappend newhbond [lindex $hbondsingleframe 2 $j ]
>
> set hbondexist [lsearch $hbondallframes $newhbond]
>
> if { $hbondexist == -1 } {
>
> lappend hbondallframes $newhbond
>
> lappend hbondcount 1
>
> } else {
>
> lset hbondcount $hbondexist [expr { [lindex $hbondcount
> $hbondexist] + 1 } ]
>
> }
>
> }
>
> }
>
>
>
> for { set i 0 } { $i < [llength $hbondallframes] } { incr i } {
>
> set donor [lindex $hbondallframes $i 0]
>
> set acceptor [lindex $hbondallframes $i 1]
>
> set occupancy [expr { 100*[lindex $hbondcount
> $i]/($numberofframes+0.0) } ]
>
>
>
> if {$i == 0 } {
>
> puts $outfile "donor \t acceptor \t occupancy"
>
> }
>
> # puts $outfile [format "%s%i \t %s%i \t\t %.2f%%" [[atomselect
> $mol "index $donor"] get resname] [[atomselect $mol "index $donor"]
> get resid] [[atomselect $mol "index $acceptor"] get resname]
> [[atomselect $mol "index $acceptor"] get resid] $occupancy]
>
> puts $outfile [format "%s%i%s \t %s%i%s \t\t %.2f%%"
> [[atomselect $mol "index $donor"] get resname] [[atomselect $mol
> "index $donor"] get resid] [[atomselect $mol “index $donor”] get name]
> [[atomselect $mol "index $acceptor"] get resname] [[atomselect $mol
> "index $acceptor"] get resid] [[atomselect $mol “index $acceptor”] get
> name] $occupancy]
>
> }
>
> if {[string compare $outfile stdout]} {
>
> close $outfile
>
> }
>
>
>
> molinfo $mol set frame $framenumberbackup
>
> }
>
> # use this thing
>
>
>
> #gly-gly
>
> set mymol [mol new h5anw.pdb type pdb waitfor all]
>
> mol addfile h5a100nw.dcd type dcd first 1 last 379 waitfor all molid
> $mymol
>
> set gly [atomselect top "resname 0SA 3LA 3LB 6LA 6LB 3YB 4YB 4GB ROH"]
>
> set pr [atomselect top "protein"]
>
> hbond_occupancy 3.5 61 $gly none hbondh5agl.dat $mymol
>
> mol delete $mymol
>
>
>
>
> ------------------------------------------------------------------------
> Access your email online and on the go with Windows Live Hotmail. Sign
> up today.
> <http://windowslive.com/Explore/Hotmail?ocid=TXT_TAGLM_WL_hotmail_acq_access_112008>