VMD-L Mailing List

From: Irene Newhouse (einew_at_hotmail.com)
Date: Mon Mar 24 2008 - 19:26:00 CDT

Thank you very much! This does work, except that VMD doesn't like "mol new"
because, the error message says, there's no mol loaded. So I changed that to "mol load",
obsolescent though it may be.
 
Irene> Date: Mon, 24 Mar 2008 18:43:28 -0400> From: akohlmey_at_cmm.chem.upenn.edu> To: einew_at_hotmail.com> CC: vmd-l_at_ks.uiuc.edu> Subject: Re: vmd-l: script to read dcd frame-at-a-time & write a selection> > 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> _________________________________________________________________
_________________________________________________________________
Windows Live Hotmail is giving away Zunes.
http://www.windowslive-hotmail.com/ZuneADay/?locale=en-US&ocid=TXT_TAGLM_Mobile_Zune_V3