proc hbond_occupancy { {dis 3} {ang 20} {sel1 protein} {sel2 none} {outfile stdout} {mol top} } {
  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 [atomselect $mol "index [lindex $hbondallframes $i 0]"]
      set acceptor [atomselect $mol "index [lindex $hbondallframes $i 1]"]
      set occupancy [expr { 100*[lindex $hbondcount $i]/($numberofframes+0.0) } ]
      if { ([$donor get name] == "N") || ([$donor get name] == "CA") || ([$donor get name] == "C") || ([$donor get name] == "O") } {
         set donortype "Main"
      } else {
         set donortype "Side"
      } 
      if { ([$acceptor get name] == "N") || ([$acceptor get name] == "CA") || ([$acceptor get name] == "C") || ([$acceptor get name] == "O") } {
         set acceptortype "Main"
      } else {
         set acceptortype "Side"
      } 

      if {$i == 0 } {
         puts $outfile "donor \t\t acceptor \t occupancy"
      }
      puts $outfile [format "%s%i-%s \t %s%i-%s \t %.2f%%" [$donor get resname] [$donor get resid] $donortype [$acceptor get resname] [$acceptor get resid] $acceptortype $occupancy]
  }

  if {[string compare $outfile stdout]} {
     close $outfile
  }

  molinfo $mol set frame $framenumberbackup
}