VMD-L Mailing List

From: Dong Luo (us917_at_yahoo.com)
Date: Wed Apr 21 2010 - 10:38:30 CDT

Hi Yongping,

I bet it's not available as an option in VMD yet. The only thing you can do is to reorgnize the water atoms that you are interested to make them all in one pbc box. I think a tcl script can do the job.

It can be something like this:

set all [atomselect top all]
set center [measure center $all] #get the center of the water box
$all delete

pbc set {9.89 9.89 9.89} -all
set se1 [atomselect top "name O and pbwithin 3.0 of name O"] # considering the period boundary condition
$sel writepdb sel.pdb
$sel delete

#now to move the atoms outside box back to box
#here I assume 9.89x9.89x9.89 is the size of the water box
mol load pdb sel.pdb
set atomX1 [atomselect top "x < [expr {[lindex $center 0] - (9.89 / 2.0)}]"]
$atomX1 moveby {9.89 0 0}
set atomY1 [atomselect top "y < [expr {[lindex $center 1] - (9.89 / 2.0)}]"]
$atomY1 moveby {0 9.89 0}
set atomZ1 [atomselect top "z < [expr {[lindex $center 2] - (9.89 / 2.0)}]"]
$atomZ1 moveby {0 0 9.89}
set atomX2 [atomselect top "x > [expr {[lindex $center 0] + (9.89 / 2.0)}]"]
$atomX2 moveby {-9.89 0 0}
set atomY2 [atomselect top "y > [expr {[lindex $center 1] + (9.89 / 2.0)}]"]
$atomY2 moveby {0 -9.89 0}
set atomZ2 [atomselect top "z > [expr {[lindex $center 2] + (9.89 / 2.0)}]"]
$atomZ2 moveby {0 0 -9.89}
$atomX1 delete
$atomY1 delete
$atomZ1 delete
$atomX2 delete
$atomY2 delete
$atomZ2 delete

set sel [atomselect top all]

hbond 3.5 30.0 $sel none hb.dat
$sel delete

Dong

________________________________
From: yongpingzeng <yongpingzeng_at_163.com>
To: vmd-l_at_ks.uiuc.edu
Cc: vmd-l_at_ks.uiuc.edu
Sent: Wed, April 21, 2010 12:16:04 AM
Subject: vmd-l: Re:pbc and measure hbond

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]
        &nbs! p; }
        ! ; & nbsp;
           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
>
>
>
>
>