VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Fri Sep 01 2006 - 15:01:50 CDT

On 9/1/06, alexandra.marques_at_fc.up.pt <alexandra.marques_at_fc.up.pt> wrote:
>
> Hi
>
> Thanks. My protein is complexed with a ligand and I am interested in identify
> waters with long residence times, particularly in the active site, near the
> ligand. I thought that one way to do this would be to identify the
> waters frame
> by frame around the ligand. Do you know if there is an easy or more "correct"
> way (or a script...) to identify waters with long residence times across a
> trajectory?

well, you could do this from _within_ VMD. i would do something along the lines
of the following.

initialize an array for counting the residence time with zeroes. one entry for
each water in the system.

set allres <total number of residues>
for {set i 0} {$i < $allres} {incr i} {
   set hist($i) 0
}

then build the required selection(s) and '$sel get residue' the indices of the
waters around the ligand. compare this list to the list of the previous
frame (initialized to {} ) and now if an index is new or present in the
previous list increment the entry in the array: incr $hist($res)
if an entry from the previous list is missing in the new the water has
left, so you store the info via: lappend restime $hist($res)
in a list and reset the counter via. set hist($res) 0
now with the list in $restime, you can build a histogram or
compute the average. voila.

you can find something similar in the mk3drama script (just for
a 2-d histogram) at:
http://www.ks.uiuc.edu/Research/vmd/script_library/scripts/mk3drama/

cheers,
    axel.

>
>
> Many thanks
> Alex
>
>
>
>
> Quoting Axel Kohlmeyer <akohlmey_at_cmm.chem.upenn.edu>:
>
> > On 8/31/06, alexandra.marques_at_fc.up.pt <alexandra.marques_at_fc.up.pt> wrote:
> >>
> >> Hi
> >>
> >> I still have a problem with my script to write the list of water
> >> molecules for
> >> each frame. The problem is with the format of the output file, that is, for
> >> each frame, the waters are listed in the horizontal and vertical lines. For
> >> example:
> >>
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >> Frame 0, 4561 4593 4614 4639 5215 5216 5220 7142 7186 7201 8694
> >>
> >> As my trajectory has 5000 frames the output file is very big and hard to
> >> analyse. Does anyone knows how should I modify the script in order to the
> >> waters be listed for each frame only in one line?
> >> The script is:
> >>
> >> set outfile [open watersTDT.txt w]
> >> set sel [atomselect top "name O and resname WAT and within 3.5 of
> >> residue 255"]
> >> set n [molinfo top get numframes]
> >> for { set i 0 } { $i < $n } { incr i } {
> >> $sel frame $i
> >> $sel update
> >> set residuelist [$sel get residue]
> >> puts "list of waters within 3.5 of residue 255:"
> >> puts "$residuelist"
> >
> >> foreach res $residuelist {
> >> puts $outfile "Frame $i, $residuelist"
> >
> > why do you output $residuelist here and not $res ??
> > this loop is redundant. either write $residuelist once
> > or write $res in each line.
> >
> > i'm also wondering why you need to output this and cannot
> > do the subsequent analysis right from withing vmd anyways.
> >
> > almost anything you could do with an external program you
> > can do with tcl as well, and even though tcl is signifiantly slower
> > than compiled code, the time needed to write out and then
> > read in and parse formatted text is significant, too.
> >
> > cheers,
> > axel.
> >
> >
> >> }
> >> }
> >> close $outfile
> >>
> >>
> >> Thank you very much for your help
> >> Alex
> >>
> >>
> >> -------------------------------------------------------------
> >> A FCUP utiliza o sistema de webmail Horde/IMP (www.horde.org)
> >>
> >> Visite: http://www.fc.up.pt/
> >>
> >>
> >
> >
> > --
> > =======================================================================
> > 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.
> >
>
>
>
> -------------------------------------------------------------
> A FCUP utiliza o sistema de webmail Horde/IMP (www.horde.org)
>
> Visite: http://www.fc.up.pt/
>
>
> 

-- 
=======================================================================
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.