From: Alberto Sergio Garay (sgaray_at_fbcb.unl.edu.ar)
Date: Tue Jun 30 2009 - 08:02:57 CDT

Hi all

I have a funny problem. I'm running a script which selects some
particular residue and wraps the rest of the atoms around it.
After that the script selects some residues and some water molecules
to look for hBonds between them. The script informs that there were
not any hbonds between this two selections. I checked visually that
the selected residues were corrects i.e. all the atoms selected were
close each other. The strange is when I select the residues
(previously wrote down by my tcl/tk script) in the tcl/tk window and
run hbonds measure I can find the hbonds, which I hoped.

Could any one give any clue of what I'am doing wrong?
Thank you in advance.

Below I paste part of my script.....

#Please, enter the residue to use for the centre of the calculation
set resid 60
#Enter the radius around the residue to make the calculations (in Angstroms)
set radius 13
#set radius for water
set radius2 [expr $radius + 5]

##################################################################################
#Enter the name for the output file
set filename1 "Water_Bridges_dmpc60_13A.csv"

# set outputPath {}; # Output path for results files. If it is empty,
output is written to the same location as an input file
# set outputPath {/Temp/ceramide}

#Set H-bond definition criteria
set cutoff 3.51
set angle 30.1

        #Extract frames from file
        set num_steps [molinfo top get numframes]

        #Open files for writing

        set fid1 [open $filename1 w]
        set fid1Log [open $filename1.log w]

        # Create file header
        puts $fid1
"frame_write\tNlipids\tWaters\twater_donors\twater_acceptors\twater_hBonds\twater_1hBonds\twater_2hBonds\twater_3hBonds\twater_4hBonds\twater_m4hBonds\twater_avghBonds\tself_hBridges\tinter_hBridges\tinter_1hBridges\tinter_2hBridges\tinter_3hBridges\tinter_4hBridges\tinter_m4hBridges\tinter_avghBridges\tinter_dahBonds\tinter_ddhBonds\tinter_aahBonds"
# puts $fid1Log "Output file $filename1 initiated"

######################################################
#Wraping the trayectory around the selected residue ##
######################################################

        #Selecting the residue to use as a centre for the calculation
        #It can be changed to make focus on some part of the residue, e.g.:
only the lipid's head
        puts "###################################"
        puts "WRAPING TRAYECTORY .........."
        #Wrapping the all molecules around the selected residue
        pbc wrap -centersel "resid $resid" -center com -all

#####################################################################################
#####################################################################################

        for {set frame 1} {$frame < $num_steps} {incr frame} {

        puts "Starting calculus on frame $frame"

        #Making a list of water and lipids molecules inside a sphere defined
around a residue

set water [atomselect top "same residue as (resname SOL) and within
$radius2 of (resid $resid and name CD)" frame $frame]
#selecting some particular atoms.....
set lipids [atomselect top "(resname DMP) and (name CD and within
$radius of (resid $resid and name CD))" frame $frame]

#write the number of residues found
puts "water: [llength [lsort -integer -unique [$water get residue]]]"
puts "water: [lsort -integer -unique [$water get residue]]"
puts "lipids: [llength [lsort -integer -unique [$lipids get residue]]]"
puts "lipids: [lsort -integer -unique [$lipids get residue]]"

set Nwaters [llength [lsort -integer -unique [$water get residue]]]
set Nlipids [llength [lsort -integer -unique [$lipids get residue]]]
puts "Nwaters: $Nwaters"
puts "Nlipids: $Nlipids"

        set LipidsMoleculesNumber [llength [lsort -integer -unique [$lipids
get residue]]]

                #Count H-bonds in selection
                set lipid_water_hBonds [measure hBonds $cutoff $angle $lipids $water]
                set water_lipid_hBonds [measure hBonds $cutoff $angle $water $lipids]

        $water delete
        $lipids delete

        #writing the number of ceramide_water_hBonds found
        set n 0
        foreach a [lindex $ceramide_water_hBonds 0] {
        set n [expr $n + 1]
        }
        puts "lipids_water_hBonds found:# $n #"

        #writing the number of water_ceramide_hBonds found
        set n 0
        foreach a [lindex $water_ceramide_hBonds 0] {
        set n [expr $n + 1]
        }
        puts "water_lipids_hBonds found:# $n #"

-- 
Dr. Sergio Garay
Facultad de Bioquimica y Cs. Biológicas
Universidad Nacional del Litoral
Santa Fe - Argentina
C.C. 242 - Ciudad Universitaria - C.P. S3000ZAA
Argentina
Ph. +54 (342) 4575-213
Fax. +54 (342) 4575-221