From: alexandra.marques_at_fc.up.pt
Date: Mon Mar 26 2007 - 15:49:19 CDT

Hi

I am trying to calculate the number of hydrogen bonds of a protein during an MD
simulation but my script don’t work properly. Could someone please help me?
 Here is the script:

set outfile [open "hbond.out" w]
set nf [molinfo top get numframes]
set allsel [atomselect top "all and name CA"]
set residlist [lsort -unique -integer [$allsel get residue]]

for {set i 0} {$i < $nf} {incr i} {

       foreach res $residlist {
       set sel1 [atomselect top "residue $res and name \"N.*\" \"O.*\" \"S.*\""]
       set sel2 [atomselect top "protein and not residue $res and name \"N.*\"
\"O.*\" \"S.*\""]
       set hbonds [measure hbonds 3.5 60 $sel1 $sel2]
       $sel1 frame $i
       $sel2 frame $i
       $sel1 update
       $sel2 update
       $sel1 delete
       $sel2 delete
       puts $outfile "number of h-bonds in frame $i: $hbonds"
      }
     }
close $outfile

Thanks
Alexandra

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

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