From: Giovanni Bellesia (
Date: Wed Mar 05 2008 - 21:26:22 CST

Hi Axel,
> 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.
yes, this is a very good call ! I totally missed the fact that
PBC do not apply when I'm using the "limited selections". Thanks !
As a matter of fact, increasing the "within" parameter for hydrogens
gives me a more reasonable decay of the perturbation (actually the decay
is ~ 2 times faster than the one I had before, and this makes sense to me).
The number integrals are "stable" over a reasonable range of "upper
cutoffs" for the hydrogens (I'm using 2.5 since I know that this is some
sort of critical distance for the gdr(OH) )
> 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))
About this second point, I agree with you that going for it leads to a
more precise definition of the perturbation (especially for larger
distances). However, neglecting the lower cutoff maybe better for a
qualitative illustration of the results.
BTW, the difference is not so dramatic after all (considering that I'm
measuring an average number of Hbonds), so I'll try to plot both data
and see what I like better :) Thanks for this one too.

Finally, I was actually thinking of using this idea of partitioning the
space in different "spherical slices" to see if there is some gradient
in the number density of water in the vicinity of the peptides. I don't
have a clear idea about how to do this at the moment ... any suggestions ?


> 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.
> cheers,
> axel.
> GB>
> GB>
> 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>
> 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>
> 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>
> GB> exit
> GB>
> GB>
> GB>