VMD-L Mailing List

From: yongpingzeng (yongpingzeng_at_163.com)
Date: Tue Apr 20 2010 - 23:16:04 CDT

Hi,
I would like to count all hbonds for pure water(48 water molecules in 9.89A cubic box). i am testing my script as below.
pbc set {9.89 9.89 9.89} -all
set sel1 [atomselect top "name O and pbwithin 3.0 of name O"] # considering the period boundary condition
hbond 3.5 30.0 $sel none hb.dat
But still it dosnt work and only count partial hydrogen bonds. Thanks in advance!
Best regards
Yongping Zeng
Yangzhou University
proc hbond { {dis 3} {ang 20} {sel1 protein} {sel2 none} {outfile stdout} {mol top} } {
   #usage: set sel1 [atomselect top "name O"]
   #usage: source hbondtest.tcl
   #usage: hbond_occupancey 3.0 20.0 $sel1 none hbond.dat #3.0 is distance of O-O, 20.0 is the angle of O-O-H
  
     set mol [molinfo top]

  if {[string compare $outfile stdout]} {
     set outfile [open $outfile w];
  }
 
  if { ! [string compare $sel1 protein]} {
     set sel1 [atomselect $mol $sel1]
  }
  set framenumberbackup [molinfo $mol get frame]
  set numberofframes [molinfo $mol get numframes]
# set framenumberbackup 10
# set numberofframes 10
       
  for { set i 0 } { $i < $numberofframes } { incr i } {
      molinfo $mol set frame $i
       set hbondcount0 0
       set hbondcount1 0
       set hbondcount2 0
       set hbondcount3 0
       set hbondcount4 0
      
         set hbondsingleframe [measure hbonds $dis $ang $sel1]
   
        set newhbond {}
      for { set j 0 } { $j < [llength [lindex $hbondsingleframe 0] ] } { incr j } {
             
          lappend newhbond [lindex $hbondsingleframe 0 $j ]
         
          lappend newhbond [lindex $hbondsingleframe 1 $j ]
      }
       
       for { set k 0 } { $k < 32 } { incr k } {
       set hbondexist [lsearch -all $newhbond $k]
  # set hbondexist1 [list [lsearch -all $newhbond $k]]
       set len [llength $hbondexist]
      
      if {$len==0} {
        set hbondcount0 [expr $hbondcount0 + 1]
          }
 
          if {$len==1} {
        set hbondcount1 [expr $hbondcount1 + 1]
          }
          
          if {$len==2} {
           set hbondcount2 [expr $hbondcount2 + 1]
           }
           
           if {$len==3} {
          set hbondcount3 [expr $hbondcount3 + 1]
           }
           
           if {$len==4} {
           set hbondcount4 [expr $hbondcount4 + 1]
           }
        
          
  
  }
   puts $outfile "$i $hbondcount0 $hbondcount1 $hbondcount2 $hbondcount3 $hbondcount4 "
  }
 
}

Dear VMD users:
 I encounter with the issue of handling periodic boundary conditions while detecting Hydrogen bonds for pure water. For instance, the command 'measure hbonds' is not aware of PBC, and will miss whatever bonds may extend across the edges of the periodic box.
Jonh Stone have given some hints:
In the post in 2005
> existing VMD functions that do:
> - 'within' and 'exwithin' atom selections
> - 'measure contacts'
> - 'measure hbonds'
> - automatic bond search (allow bonds across PBC cell boundaries)
> - draw bond components of the Bonds, CPK, and Licorice representations
> - draw dynamic bonds representation
>
> Adding PBC-awareness to the drawing code is the most risky in terms
> of hurting performance, but this may be something that could be dealt with
> by adding a checkbox so that VMD only spends effort checking for this when
> the user actually needs it.
>
> Before embarking on implementing all of this, it wvmd-l_at_ks.uiuc.eduould be nice to know
> which of these users would find most beneficial. We could work on the
> most-wanted features first...
>
> John Stone

I dont know if the problem have been resolved? Could anyone like to give me some help?
 
Best
Yongping Zeng
Yangzhou University