VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Fri Mar 12 2010 - 09:20:07 CST

On Thu, 2010-03-11 at 15:54 -0500, Florentina Tofoleanu wrote:
> Hello,

hi florentina,

> You were right Axel, the first two frames in the dcd were causing the
> error. I eventually got the new "wrapped" trajectory and calculated
> the HBs with water pbwithin 10A of a segment of protein. I wanted to
> compare the no. of HBs between outer and inner segments of protein and
> water around them. A curious result was that I got more HBs for the
> less exposed segments, and increasing the region "pbwithin" from 10 to
> 12, 14, 16 and 18 A around the segments I got similar results for
> different segments, still less HBs for the outer ones. The fact is
> these outer segments were sticking out of the water during the
> simulation. Could the lower number of HBs for a more exposed to water
> segment be an error of the script or of VMD? Why does the region I
> take water into consideration, from 10 to 18A, matter when I have a
> cutoff distance for the HBs of 3.1A?

please remember, that pbwithin selects atoms at their
*original* position. using a pbwithin selection in combination
with "measure hbonds" is currently not a good idea, since
measure hbonds is not PBC aware. it also looks as if your
system is not very large, so with a large cutoff for pbwithin,
you could be including atoms that are much closer.

in these cases it is often very helpful to visually check
the result of a selection. e.g. create a new representation
set it to a vdw rep with an unusual fixed color (orange)
and transparent material, then try out your selection, and
check whether you get the atoms selected that you expect.

nothing will be "sticking" out in a PBC simulation unless
your trajectory is severely messed up. it may just be an
indication that your water shell is fairly small and that
you cannot get away with using pbc wrap globally, but that
you have to wrap the waters separately for each region
of your protein. however, you should also be careful with
your results in that case, since too small a water area
around your protein, can negatively affect the applicability
of your results.

cheers,
   axel.

> Or does anybody have any explanation for it? How does the wrap
> function modify the way water "sees" where the protein atoms are?
>
> Again, the script for getting the HBs:
>
> set outfile [open output.dat w];
>
> set numframes [molinfo $mol get numframes]
>
> set refA2 [atomselect $mol "protein and segname A"]
> set refB2 [atomselect $mol "water pbwithin 18 of protein and segname
> A"]
>
> for {set frame 0} {$frame < $numframes} {incr frame} {
>
> #$refA2 frame $frame
> $refB2 frame $frame
>
> #$refA2 update
> $refB2 update
>
>
> set nhb2 [llength [lindex [measure hbonds 3.1 60 $refA2 $refB2] 0]]
>
> puts $outfile "$frame $nhb2"
>
> }
>
> close $outfile
>
> On Thu, Mar 11, 2010 at 12:51 PM, Axel Kohlmeyer <akohlmey_at_gmail.com>
> wrote:
> On Thu, 2010-03-11 at 12:40 -0500, Florentina Tofoleanu wrote:
>
> > Thank you Axel for your answer. I added that line,
> >
> > pbc wrap -molid $mol -compound res -centersel $refA -first
> $frame
> > -last $frame
> >
> > but I get this error:
> >
> > -error Suspicious pbc side length (a=0.000000 b=0.000000
> c=0.000000).
> > Have you forgotten to set the pbc parameters?
> > domain error: argument not in valid range
>
>
> this means, that one of more frames of your .dcd file have
> no information about the size of your simulation cell.
> btw, this would also cause problems in the case of using
> a selection containing "pbwithin".
>
> check out the pbctools plugin documentation on how to
> print, read, and apply cell size information using pbctools.
>
> cheers,
> axel.
>
> >
>
>
> > Only "pbc wrap -molid $mol -compound res -centersel $refA"
> works fine,
> > but just for the last frame. I wrote this little script to
> get the
> > result only for the "wrap" function. Could you explain me
> why it
> > doesn't change all the frames even though it is in a "for"
> cycle? Is
> > it related to the inutility of "$refA frame $frame", being
> the same
> > selection in every frame you need to specify the frames you
> want to
> > change?
> >
> > set mol [mol new file.psf type psf waitfor all]
> > mol addfile file.dcd type dcd waitfor all molid $mol
> >
> > set numframes [molinfo $mol get numframes]
> > set all [atomselect top all]
> > set refA [atomselect $mol "protein"]
> >
> > for {set frame 0} {$frame < $numframes} {incr frame} {
> >
> > #$refA frame $frame
> >
> > pbc wrap -molid $mol -compound res -centersel $refA
> > #pbc wrap -molid $mol -compound res -centersel $refA -first
> $frame
> > -last $frame
> >
> > animate write dcd wrap.dcd sel $all
> > }
> >
> > Could you tell me what I need to change to make it work?
> >
> > Regards,
> >
> > Florentina
> >
> > On Wed, Mar 10, 2010 at 6:41 PM, Axel Kohlmeyer
> <akohlmey_at_gmail.com>
> > wrote:
> > On Wed, Mar 10, 2010 at 5:27 PM, Florentina
> Tofoleanu
> > <florentina.tofoleanu_at_gmail.com> wrote:
> > > Hello,
> >
> > hi florentina,
> >
> > [...]
> >
> > > $refA update
> >
> > this update is not really needed. this selection
> will not
> > change.
> >
> > > $refB update
> > >
> > > set nhb [llength [lindex [measure hbonds 3.1 60
> $refA $refB]
> > 0]]
> >
> >
> > here is the culprit. while $refB will select
> positions with
> > PBC in mind,
> > it actually does not include multiple copies of the
> > coordinates and thus
> > measure hbonds will only compute the distances to
> the original
> > coordinates of all particles. i.e. it is not PBC
> aware and
> > thus you
> > get in both cases the same result.
> >
> > if you cell is large enough, you might want to try
> running
> > pbc wrap with providing $refA (or rather its center
> of mass)
> > as the reference point for the center of the
> wrapping.
> > this way your "within" selection should be ok.
> >
> > >
> > > puts $outfile "$frame $nhb"
> > >
> > > }
> > >
> > > close $outfile
> > >
> > >
> > > Is there an error in the script, or just the
> pbwithin
> > command won't do?
> > >
> > > I've seen in this post
> > >
> >
> http://www.ks.uiuc.edu/Research/vmd/mailing_list/vmd-l/14664.html
> > > "(Pbwithin) does not, however, modify the
> coordinates - they
> > will still be
> > > the coordinates of the atom in the primary image,
> so if you
> > want to compute
> > > distances, you need to take into account the
> periodic
> > boundary conditions
> > > manually. " Could Olaf, or anybody else, explain
> what taking
> > into account
> > > the pbc manually mean and how can I change my
> script?
> >
> >
> > how about adding?
> >
> > pbc wrap -molid $mol -compound res -centersel $refA
> -first
> > $frame -last $frame
> >
> > (add -orthorhombic for more speed, if your cell is
> > orthorhombic)
> >
> > cheers,
> > axel.
> >
> >
> > >
> > > Thank you,
> > >
> > > Florentina
> > >
> > >
> > > On Sat, Sep 19, 2009 at 9:55 AM, Axel Kohlmeyer
> > <akohlmey_at_gmail.com> wrote:
> > >>
> > >> 2009/9/19 Thomas Evangelidis
> <te8624_at_mbg.duth.gr>:
> > >> > Hi Katherine,
> > >> >
> > >> > lets say you want to count all the waters(TIP3)
> within 5
> > A of resids
> > >> > 1-50
> > >> > taking into account waters in the nearest
> periodic cells,
> > then the
> > >> > command
> > >> > should be:
> > >> >
> > >> > llength [[atomselect top "(oxygen pbwithin 5 of
> (protein
> > and resid 1 to
> > >> > 50))
> > >> > and resname TIP3"] list]
> > >>
> > >> folks,
> > >>
> > >> when giving code examples, please make sure that
> you give
> > examples
> > >> of _good_ coding. this piece of code has a high
> risk of
> > creating a memory
> > >> leak
> > >> and is inefficient.
> > >>
> > >> the same can be done by:
> > >>
> > >> set sel [atomselect top "(oxygen pbwithin 5 of
> (protein and
> > resid 1 to
> > >> 50)) and resname TIP3"]
> > >> puts "found [$sel num] water molecules"
> > >> $sel delete
> > >>
> > >> and this code uses the "num" option to do the
> counting in c
> > rather
> > >> than creating a tcl
> > >> list and counting its length and show that the
> functions
> > created by
> > >> atomselect have to
> > >> be deleted after use (it is better to always
> manually
> > delete them
> > >> rather than counting
> > >> on tcl to clean up).
> > >>
> > >> > I have an alpha version of a plugin to truncate
> > trajectories and keep
> > >> > only a
> > >> > specified number of waters within R Angstroms
> of an
> > atomselection, it
> > >> > might
> > >> > be helpful to your analysis.
> > >>
> > >> that reminds me. having you been able to sort out
> your
> > issues with bigdcd?
> > >> did the improved re-rentrant version of bigdcd
> work for
> > you?
> > >>
> > >>
> > >> thanks,
> > >> axel.
> > >> >
> > >> > regards,
> > >> > Tom
> > >> >
> > >> >
> > >> >> Dear VMD users:
> > >> >> My sytem contains 18 proteins in a unit cell.
> I also
> > have several 20ns
> > >> >> trajectories (*.dcd) of this system.
> > >> >> I would like to calculate the distribution of
> the water
> > molecules
> > >> >> around
> > >> >> each one of the proteins and around the
> different
> > amoinoacids that
> > >> >> compose
> > >> >> them during the trajectories.
> > >> >> I've been trying to find another VMD mail list
> posted
> > message
> > >> >> containing
> > >> >> this question and somewhere in the VMD user's
> guide but
> > I haven't found
> > >> >> anything related.
> > >> >> I've worked with some scripts in VMD, using
> the
> > atomselect command and
> > >> >> doing
> > >> >> loops along trajectories.
> > >> >> What I need is some help with the right
> command to count
> > water
> > >> >> molecules
> > >> >> (TIP3) around a specific part of the protein.
> > >> >> Thank you in advance.
> > >> >> KP
> > >> >>
> > >> >
> > >> >
> > >> > ----- End message from kparra_at_mail.usf.edu
> -----
> > >> >
> > >> >
> > >> >
> > >> >
> > >>
> > >>
> > >>
> > >> --
> > >> Dr. Axel Kohlmeyer akohlmey_at_gmail.com
> > >> Institute for Computational Molecular Science
> > >> College of Science and Technology
> > >> Temple University, Philadelphia PA, USA.
> > >
> > >
> > >
> > > --
> > > Florentina Tofoleanu
> > > Postgraduate Research Fellow
> > > Theoretical and Computational Biophysics Group
> > > University College Dublin
> > > School of Physics, Rm. 110
> > > Belfield, Dublin 4, Ireland
> > >
> >
> >
> >
> >
> > --
> > Dr. Axel Kohlmeyer akohlmey_at_gmail.com
> >
> > http://sites.google.com/site/akohlmey/
> >
> > Institute for Computational Molecular Science
> >
> >
> > Temple University, Philadelphia PA, USA.
> >
> >
> >
> >
> > --
> > Florentina Tofoleanu
> > Postgraduate Research Fellow
> > Theoretical and Computational Biophysics Group
> > University College Dublin
> > School of Physics, Rm. 110
> > Belfield, Dublin 4, Ireland
>
>
> --
>
> Dr. Axel Kohlmeyer akohlmey_at_gmail.com
> http://sites.google.com/site/akohlmey/
>
> Institute for Computational Molecular Science
> Temple University, Philadelphia PA, USA.
>
>
>
>
>
>
> --
> Florentina Tofoleanu
> Postgraduate Research Fellow
> Theoretical and Computational Biophysics Group
> University College Dublin
> School of Physics, Rm. 110
> Belfield, Dublin 4, Ireland 

-- 
Dr. Axel Kohlmeyer  akohlmey_at_gmail.com
http://sites.google.com/site/akohlmey/
Institute for Computational Molecular Science
Temple University, Philadelphia PA, USA.