VMD-L Mailing List

From: Irene Newhouse (einew_at_hotmail.com)
Date: Fri Nov 21 2008 - 12:57:03 CST

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 theligand atoms involved in the H-bond. As the original script was intended for usewith proteins, where you know which atom is involved, given a residue, the originalscript 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.
http://windowslive.com/Explore/Hotmail?ocid=TXT_TAGLM_WL_hotmail_acq_access_112008