VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Tue Jul 29 2008 - 20:40:30 CDT

On Tue, 29 Jul 2008, Christopher Gillespie wrote:

CG> Hello again,
CG> Thank you Eduard and John for the help.
CG>
CG> Using Eduard's suggestion for surface selection I modified my script to set
CG> only protein atoms within 2.8A of water, updated at each frame. The problem
CG> is that it takes until ~40A to get g(r) to 1, specifically g(r) is always <
CG> 1 until it nears 40A. I believe this is due to the averaging including the
CG> shell area that is blocked by the rest of the protein atoms, i.e. the shel
CG> volume should be less then 4*pi*r^2*dr to remove the volume occupied by the
CG> protein. Has someone already tackled this issue.

the problem is that there is no straightforward way to determine the
volume that a _single_ molecule (e.g. a protein) occupies. you can only
determine the _total_ density and based on that an _average_ volume.
but that only works if you have only a single atom or molecule type.
as soon as you mix items of different size you have all kinds of
non-linear behavior, e.g. the mixture requiring less total volume
since the smaller particles can go into places that would not be
accessable to larger ones. how do you want to account for that in
a systematical way?

[...]

CG> set dr 0.002 ;# bin size (Ang.) Shortened due to Axel Kohlmeyer's
CG>
CG> set dim 40 ;# rdf Range (Ang.)

the bugs related to normalization should be fixed in
the latest alpha release.

depending on how large the excluded volume from the protein is
relative to the rest, you should get a behavior of the g(r) that
in some mid-rang g(r) will be > 1.0 and in the long range g(r)
will be < 1.0. you can rationalize this with statistical arguments.
you don't have any water molecules where the protein is, but you
are computing your reference density as if you had a homogeneous
distribution. yet you only compute g(r) for molecules where you
have water and it is more probably to have waters in the neighborhood
than what you would compute as total density, whereas the void would
be farthest away from the waters in the center of the water area.
at least this is the behavior if you have significant density
fluctuations (e.g. in supercritical water). not sure how much it
shows in your case.

hope this helps,
    axel.

CG>
CG>
CG>
CG> set out1 [open gofrSurfWat.dat w] ;# Output file for water rdf
CG>
CG> set out2 [open gofrSurfNa.dat w] ;# Output file for sodium rdf
CG>
CG> set out3 [open gofrSurfCl.dat w] ;# Output file for chloride rdf
CG>
CG> set frm [molinfo top get numframes] ;# Determine number of loaded frames
CG>
CG> set i [expr $frm -1]
CG>
CG> puts " Number of frames to process : $i"
CG>
CG> # Sodium
CG>
CG> puts ""
CG>
CG> puts "Starting with Sodium"
CG>
CG> set grN [measure gofr $surf $cat delta $dr rmax $dim usepbc 1 selupdate 1
CG> first 1 last -1 step 1]
CG>
CG>
CG>
CG> set rN [lindex $grN 0]
CG>
CG> set gN [lindex $grN 1]
CG>
CG>
CG>
CG> foreach j $rN k $gN {
CG>
CG> puts $out2 "$j $k"
CG>
CG> }
CG>
CG>
CG>
CG> # Chloride
CG>
CG> puts ""
CG>
CG> puts "Now Chloride"
CG>
CG> set grC [measure gofr $surf $ani delta $dr rmax $dim usepbc 1 selupdate 1
CG> first 1 last -1 step 1]
CG>
CG>
CG>
CG> set rC [lindex $grC 0]
CG>
CG> set gC [lindex $grC 1]
CG>
CG>
CG>
CG> foreach j $rC k $gC {
CG>
CG> puts $out3 "$j $k"
CG>
CG> }
CG>
CG>
CG>
CG> # Water
CG>
CG> puts ""
CG>
CG> puts "Finally Water"
CG>
CG> set grO [measure gofr $surf $sol delta $dr rmax $dim usepbc 1 selupdate 1
CG> first 1 last -1 step 1]
CG>
CG>
CG>
CG> set rO [lindex $grO 0]
CG>
CG> set gO [lindex $grO 1]
CG>
CG>
CG>
CG> foreach j $rO k $gO {
CG>
CG> puts $out1 "$j $k"
CG>
CG> }
CG>
CG>
CG> puts " "
CG>
CG> puts "All Done"
CG>
CG> puts " "
CG>
CG>
CG>
CG>
CG> On Mon, Jul 21, 2008 at 11:22 AM, John Stone <johns_at_ks.uiuc.edu> wrote:
CG>
CG> >
CG> > Beyond the reply that Eduard already sent, I just wanted to mention
CG> > that you can use a script to compute SASA per residue, assign the
CG> > per-residue SASA to "beta" or "user" (Just as an example. You could
CG> > just as easily use any of the per-atom fields) and then use a selection
CG> > like:
CG> > beta > 0
CG> >
CG> > Cheers,
CG> > John Stone
CG> > vmd_at_ks.uiuc.edu
CG> >
CG> > On Sun, Jul 20, 2008 at 05:28:04PM -0400, Christopher Gillespie wrote:
CG> > > Hello,
CG> > > I am trying to generate some rdfs for systems that I have already run. I
CG> > > want to do these using the protein surface as one selection. I know that
CG> > > using the selection "Protein and Surface" only implies "Not Buried" which
CG> > > removes hydrophobic residues from the selection set, but I want to see
CG> > the
CG> > > actual atoms, including any hydrophobic ones that have at least some
CG> > SASA,
CG> > > while also not including polar atoms that are actually buried, e.g. at
CG> > the
CG> > > dimer interface. I know it may not be exactly correct since I should
CG> > have
CG> > > to determine the SASA for each frame, but to make this doable in a
CG> > > reasonable amount of time and since my protein has an average RMSD of
CG> > ~1.2A
CG> > > the atoms that make up the set from frame 0 should not deviate that much,
CG> > I
CG> > > assume.
CG> > >
CG> > > So does anyone know how to define a selection of SASA > 0 atoms for use
CG> > in
CG> > > gofr?
CG> > >
CG> > > Thanks
CG> > >
CG> > > Chris
CG> >
CG> > --
CG> > NIH Resource for Macromolecular Modeling and Bioinformatics
CG> > Beckman Institute for Advanced Science and Technology
CG> > University of Illinois, 405 N. Mathews Ave, Urbana, IL 61801
CG> > Email: johns_at_ks.uiuc.edu Phone: 217-244-3349
CG> > WWW: http://www.ks.uiuc.edu/~johns/ Fax: 217-244-6078
CG> >
CG> 

-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   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.