VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Mon Mar 24 2008 - 17:43:28 CDT

On Mon, 24 Mar 2008, Irene Newhouse wrote:

irene,

IN> I'm still battling with trying to read a dcd file a frame at a time,
IN> making a selection & writing that out as a pdb file [because there
IN> may be differing numbers of atoms in the selection along the
IN> trajectory].

IN> I want this to be a self-contained script so that I can submit it in batch mode to analyze my entire, looong trajectory. Each trajectory is a trimer, we want to know how each one of those is behaving, & we have 6 separate trajectories. I don't have access to a computer than can read the entire trajectory. In fact, the cluster on which I intend to batch process can't do that either.

IN> I'm having a very hard time getting the dcd file read. This is what
IN> I have so far:

IN> #initialize
IN> set numfr 350
IN> set pref h5a100
IN> mol load parm7 $pref.prmtop rst7 $pref.inpcrd

mol load is obsolescent. please use 'mol new' instead.
also mol new returns the molecule id, so to not have to
rely on 'top' you can do:

set mol [mol new parm7 $pref.prmtop]

you don't need to load the restart since you don't want to analyze
that but your .dcd instead, right?

IN> set brdg [atomselect top "(same residue as water and within 5 of resid 1447 to 1449) \

this then would become:

set brdg [atomselect $mol ...

this is a bit more reliable, as you don't have to keep
track of what is the current top molecule. most of the
time it works, but since you want something really self-contained,
i'd recommend to go this route.

IN> or (resid 1447 to 1449) or (resid 220 to 228) or (resid 130 to 135) or (resid 180 to 190)"]

IN> #loop over all frames
IN> for {set frame 0} {$frame < $numfr} {incr frame} {
IN> # read the current frame & update the selection
IN> puts "Frame $frame ..."

# delete all remaining frames so there is no memory "leak"
animate delete beg 0 $mol

IN> animate read dcd $pref.dcd waitfor 1 top

now here you want to use 'mol addfile' instead:

mol addfile $pref.dcd type dcd first $frame last $frame waitfor all molid $mol

IN> $brdg frame $frame

since you have only one frame loaded you'd run:

$brdg frame 0

# and to recompute the selection add:

$brdg update

IN>
IN>
IN> # write the current frame into a pdb file
IN>
IN> $brdg writepdb [format $pref-%04d.pdb $frame]
IN>
IN> #end loop over all frames
IN> }
IN>
IN> I am debugging with an extraction from the full trajectory, that's 350 frames long. The last 2 lines of the VMD console read:
IN>
IN> read_dcdheader: premature end of file
IN> ERROR: could not readh h5a100.dcd

this has nothing to do with the rest of your script.
most likely either your dcd file is corrupted.

cheers,
  axel.

IN>
IN> How do I code reading 1 frame of the dcd file?
IN>
IN> Thanks!
IN> Irene Newhouse
IN> _________________________________________________________________
IN> Test your Star IQ
IN> http://club.live.com/red_carpet_reveal.aspx?icid=redcarpet_HMTAGMAR 

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