VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Wed May 20 2009 - 14:05:27 CDT

On Wed, 2009-05-20 at 19:21 +0100, Florentina Tofoleanu wrote:
> Dear Axel,

dear florentina,

i'm keeping the cc to vmd-l, so that others can
see the resolution as well.

> Thank you so much for the quick reply. You made clear many things I
> used but didn't fully understand. :-)

you are welcome.

> Ok, so my new script would be:

> set mol [mol new X.psf type psf waitfor all]
> mol addfile X.dcd type dcd waitfor all molid $mol
>
> set outfile [open com.dat w];
>
> set numframes [molinfo $mol get numframes]
>
> set refA [atomselect $mol "chain A and not hydrogen"]
> set refB [atomselect $mol "chain B and not hydrogen"]
>
> for {set frame 0} {$frame < $numframes} {incr frame} {
>
> $refA frame $frame
> $refB frame $frame
>
> set com1 [measure center $refA weight mass]
> set com2 [measure center $refB weight mass]
>
> puts $outfile "$frame [veclength [vecsub $com1 $com2]]"
> }
> close $outfile
>
>
> which gave me this error:
>
> "measure center: bad weight sum, would cause divide by zero"
>
> Could you tell me what's still wrong with the script?

i don't think there is anything wrong with the script by
itself, but rather it looks as if your .psf file would
not contain any masses. is this true?

you didn't see it before, because the reader for .pdb files
tries to guess masses from the atom names. note that it may
guess wrong on occasion (typical example: is HG1 a gamma
hydrogen or a Mercury atom?). if you trust the guessing to
work correctly, you can change the first two lines into:

set mol1 [mol new X.pdb type pdb waitfor all]
set sel1 [atomselect $mol1 all]
set mol [mol new X.psf type psf waitfor all]
mol addfile X.dcd type dcd waitfor all molid $mol
set sel2 [atomselect $mol all]
$sel2 set mass [$sel1 get mass]

# cleanup
$sel1 delete
$sel2 delete
mol delete $mol1

or rather stick with a geometric center (i.e. use no weighting).
since you exclude the hydrogens the result should be pretty
much the same (carbons are "still" 75% of the weight of oxygens).

cheers,
   axel.

>
> Regards,
>
> Florentina
>
> On Wed, May 20, 2009 at 6:09 PM, Axel Kohlmeyer
> <akohlmey_at_cmm.chem.upenn.edu> wrote:
> On Wed, 2009-05-20 at 17:34 +0100, Florentina Tofoleanu wrote:
>
> florentina,
>
> > I'm using this script to calculate it, but I'm getting the
> same value
> > for each of the frames:
> >
> >
> > mol load pdb X.pdb
> > mol load psf X.psf
> > mol load dcd X.dcd
>
>
> this is deprecated syntax and you may not get
> the desired result, since each load will create
> a new molecules and thus the selections referring
> to molecule 0 will only refer to the .pdb file which
> has only one frame, whereas the referring to the top
> molecule will give you the number of frames of the
> dcd file.
>
> what you want to do instead is:
>
> set mol [mol new X.psf type psf waitfor all]
> mol addfile X.dcd type dcd waitfor all molid $mol
>
> >
> > set outfile [open com.dat w];
> >
> > set numframes [molinfo top get numframes];
> >
> > set refA [atomselect 0 "chain A and not hydrogen"]
> > set refB [atomselect 0 "chain B and not hydrogen"]
>
>
> and this should become:
>
> set numframes [molinfo $mol get numframes]
> set refA [atomselect $mol "chain A and not hydrogen"]
> set refB [atomselect $mol "chain B and not hydrogen"]
>
>
> > for {set frame 0} {$frame < $numframes} {incr frame} {
> > $refA frame $frame
> > $refB frame $frame
> > $refA update
> > $refB update
>
>
> the update is not needed (only if you use "within"
> or selections based on position in space).
>
> > set com1 [measure center $refA weight mass]
> > set com2 [measure center $refB weight mass]
>
> > lappend simdata [veclength [vecsub $com1 $com2]]
>
> this is probably not what you want, too, as it will
> put all COM values in one line
>
> more likely you want:
>
> puts $outfile "$frame [veclength [vecsub $com1 $com2]]"
>
> > }
>
> ...and get rid of this line.
>
> > puts $outfile “$simdata”;
>
>
> > close $outfile
>
>
> > I've read the messages posted before related to this issue,
> but can't
> > really find the problem with this script. Could someone
> please show me
> > what's wrong with it and help me correcting it?
>
>
> which is hopefully done.
>
> cheers,
> axel.
>
>
> > Regards,
> > Florentina
> >
> >
> > --
> > Florentina Tofoleanu
> > Postgraduate Research Fellow
> > Theoretical and Computational Biophysics Group
> > University College Dublin
> > School of Physics, Rm. 110
> > Belfield, Dublin 4, Ireland
>
>
> --
> =======================================================================
> 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.
>
>
>
>
> --
> Florentina Tofoleanu
> Postgraduate Research Fellow
> Theoretical and Computational Biophysics Group
> University College Dublin
> School of Physics, Rm. 110
> Belfield, Dublin 4, Ireland 

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