VMD-L Mailing List

From: alexandra.marques_at_fc.up.pt
Date: Sat Feb 17 2007 - 17:52:45 CST

Hi

I have a script to calculate the hydrogen bonds between some waters of the
active site and the protein residues.However this script seems to work with
trajectories but not with pdb files. The other problem is that I cannot obtain
"the number of h-bonds in frame"; instead I obtain an output like this:
number of h-bonds in frame 0: {19711 19711 19711}
I am not seeing where the errors are in the script. Could someone please help
me?
Here is the script:

set outfile [open "hbonds.out" w]
set pep [atomselect top "protein"]
set wat [atomselect top "water and ((within 4.0 of residue 0 to 254) and (within
4.0 of residue 255))"]
set nf [molinfo top get numframes]
for {set i 0} {$i < $nf} {incr i} {
  $pep frame $i
  $wat frame $i
  $pep update
  $wat update
  set nhb [measure hbonds 3.5 120 $wat $pep]
  puts ${outfile} "number of h-bonds in frame $i: $nhb"
  foreach donor [lindex $nhb 0] acceptor [lindex $nhb 1] {
   set seldonor [atomselect top "index $donor" frame $i]
   set donorResidue [$seldonor get residue]
   $seldonor delete
   set selacceptor [atomselect top "index $acceptor" frame $i]
   set acceptorResidue [$selacceptor get residue]
   $selacceptor delete
   puts "The donor has index $donor, belongs to residue $donorResidue"
   puts "The acceptor has index $acceptor, belongs to residue $acceptorResidue"
   puts ${outfile} "frame $i: donor mask :$donorResidue"
   puts ${outfile} "frame $i: acceptor mask :$acceptorResidue"
 }
}
close $outfile

Thanks
Alexandra

-------------------------------------------------------------
A FCUP utiliza o sistema de webmail Horde/IMP (www.horde.org)

Visite: http://www.fc.up.pt/