From: Axel Kohlmeyer (
Date: Wed Mar 05 2008 - 16:34:14 CST

On Wed, 5 Mar 2008, Giovanni Bellesia wrote:

GB> Hi Axel,

hi giovanni,


GB> The idea is that these numbers can give me a measure of how much my small
GB> peptide is perturbating the water network.

ok, your script looks very reasonable. i'd actually increase the
"within" parameter for the hydrogens, because otherwise you'll
have water oxygens at the boundary of your area that are listed
as having no or fewer hydrogens in the neighborhood than they
actually have. secondly, i would then pick the water selection
to make a "hollow" sphere for the outer layers of water around
the peptide, so that you don't have the inner hydrogen bonding
averaged into the outer layers. something like:

name OH2 and (water and within 7.0 of protein) and (not (water and
within 5.0 of protein))

using some macros may make this easier.

other than that, you should actually get the kind of numbers
that you desire. please let me know how this works out. this
looks very interesting.


GB> set psf kffedimer-wb.psf
GB> set dcd T1.dcd
GB> mol load psf $psf dcd $dcd
GB> set outfile1 [open gdr-OH-50.dat w]
GB> set sel1 [atomselect top "name OH2 and water and same residue as
GB> within 5.0 of protein"]
GB> set sel2 [atomselect top "{name H1 or name H2} and water and same
GB> residue as within 5.0 of protein"]
GB> set gr0 [measure gofr $sel1 $sel2 delta 0.1 rmax 10.0 usepbc 1
GB> selupdate 1 first 0 last -1 step 1]
GB> set r [lindex $gr0 0]
GB> set gr [lindex $gr0 1]
GB> set igr [lindex $gr0 2]
GB> set isto [lindex $gr0 3]
GB> foreach j $r k $gr l $igr m $isto {
GB> puts $outfile1 [format "%.4f\t%.4f\t%.4f\t%.4f" $j $k $l $m]
GB> }
GB> close $outfile1
GB> exit

Axel Kohlmeyer
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
If you make something idiot-proof, the universe creates a better idiot.