Re: Get number of timesteps from DCD via script

From: Joshua D. Moore (joshuadmoore_at_gmail.com)
Date: Tue Feb 08 2011 - 16:56:02 CST

Here's a small fortran code which should output the header. I checked it on
a dcd file written by LAMMPS and it seems to work. I didn't check on one
written by NAMD. I'm not sure if NAMD is outputting all this information or
not. I would guess that they are.

Let me know if it solves your problem. This is basically what Andrew said
except the numbers in ICNTRL actually mean something. I think your problem
with delta was that it is actually a single precision real and not an
integer as on that page I was showing you. That page has several mistakes.
I've pointed this out a few times. delta is real(4) and the coordinates are
single precision and not double. The PBC are double and the page is right
about that.

Best.

Josh

program read_dcd_header
implicit none
character*80 dcdfilename
character*4 HDR
integer ICNTRL(20)
integer NTITL
integer j
integer natoms
character*80 TITLE(100)
real(4) :: delta

write(*,*) "What is your dcd file name?"
read(*,*) dcdfilename

open(unit=243,file=dcdfilename,form='unformatted',status="unknown")
read(243) HDR,(ICNTRL(j),j=1,9),delta,(ICNTRL(j),j=11,20)
read(243) NTITL,(TITLE(j),j=1,NTITL)
read(243) natoms

write(*,*) "# of configurations in dcd file:", ICNTRL(1)
write(*,*) "time step of first configuration:", ICNTRL(2)
write(*,*) "frequency of configurations in dcd file, dcd_freq:", ICNTRL(3)
write(*,*) "last time step in dcd file:", ICNTRL(4)
write(*,*) "timestep of simulation", delta
write(*,*) "natoms in dcd file", natoms

end program

On Tue, Feb 8, 2011 at 6:24 AM, Bjoern Olausson <namdlist_at_googlemail.com>wrote:

> On Monday 07 February 2011 19:53:56 Joshua D. Moore wrote:
> > catdcd can tell you the number of frames in the trajectory
> >
> > http://www.ks.uiuc.edu/Development/MDTools/catdcd/
> >
> > The other information is written in the header of the dcd file. You can
> > print that out if you wanted to using a simple C or fortran code.
> >
> > The format is here with one correct in that the the PBCs are written in
> > double precision while the coordinates are in single precision (not
> double
> > as is written on the page).
> >
> > http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/dcdplugin.html
> >
> > The format can also be ascertained fairly easily from the plugin in the
> VMD
> > source code.
> >
> > Josh
> >
> Interesting, thanks!
>
> So when I understand the header correctly I can get NSET (Number of overall
> frames?), integration time (DELTA), write out frequency (NSAVC) and the
> first
> timestep (not sure about the meaning of ISTRT - step 1).
>
> With these information every thing else I want is at ones hand :-)
>
> Since I am a noob regarding C or Fortran I'll played around with python and
> came up with the appended script. But after the "NSAVC" header, something
> pops
> up which is not documented on the above website. Does anyone know what it
> exactly refers to?
> Furthermore the "DELTA" looks messed up. It should be 2.0 fs.
>
> Any clue what went wrong?
>
> http://paste.pocoo.org/show/FeirO5JuypU7XzrUmQtK/
> ----------------------------------------------------
> #!/bin/env python
> # -*- coding: utf-8 -*-
> #
> #
>
> import sys, struct
>
> #HDR NSET ISTRT NSAVC 5-ZEROS
> NATOM-NFREAT
> DELTA 9-ZEROS
> #`CORD' #files step 1 step interv. zeroes
> (zero)
> timestep (zeroes)
> #C*4 INT INT INT
> 5INT INT
> DOUBLE 9INT
> #==========================================================================
> #NTITLE TITLE
> #INT (=2) C*MAXTITL (=32)
> #==========================================================================
> #NATOM
> ##atoms
> #INT
> #==========================================================================
> #X(I), I=1,NATOM (SINGLE)
> #Y(I), I=1,NATOM
> #Z(I), I=1,NATOM
> #==========================================================================
>
>
> f = open(sys.argv[1], 'rb')
> HDR = struct.unpack('i4c',f.read(struct.calcsize('i4c')))
> NSET = struct.unpack('i', f.read(struct.calcsize('i')))[0]
> ISTART = struct.unpack('i', f.read(struct.calcsize('i')))[0]
> NSAVC = struct.unpack('i', f.read(struct.calcsize('i')))[0]
> ##NOTE: The following seems not to be documented on
> ##NOTE: http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/dcdplugin.html
> UNDOC = struct.unpack('i', f.read(struct.calcsize('i')))[0]
> ZEROS5 = struct.unpack('4i', f.read(struct.calcsize('4i')))
> NFREAT = struct.unpack('i', f.read(struct.calcsize('i')))
> DELTA = struct.unpack('d',f.read(struct.calcsize('d')))
> ZEROS9 = struct.unpack('8i', f.read(struct.calcsize('8i')))
>
> ISTOP = ISTART + (NSET * NSAVC)
> NSTP = ISTOP - ISTART
>
> print "%(ID)-13s %(DATA)s" %{"ID": "HDR:", "DATA": HDR}
> print "%(ID)-13s %(DATA)i" %{"ID": "Coor Sets:", "DATA": NSET}
> print "%(ID)-13s %(DATA)i" %{"ID": "First Step:", "DATA": ISTART}
> print "%(ID)-13s %(DATA)i" %{"ID": "Save Int.:", "DATA": NSAVC}
> print "%(ID)-13s %(DATA)i" %{"ID": "UNDOC:", "DATA": UNDOC}
> print "%(ID)-13s %(DATA)s" %{"ID": "5-ZEROS:", "DATA": ZEROS5}
> print "%(ID)-13s %(DATA)s" %{"ID": "NATOM-NFREAT:", "DATA": NFREAT}
> print "%(ID)-13s %(DATA)s" %{"ID": "Int. Time:", "DATA": DELTA}
> print "%(ID)-13s %(DATA)s" %{"ID": "9-ZEROS:", "DATA": ZEROS9}
> print "\n%(ID)-13s %(DATA)i" %{"ID": "Last Step:", "DATA": ISTOP}
> print "%(ID)-13s %(DATA)i" %{"ID": "Num Steps:", "DATA": NSTP}
>
> ----------------------------------------------------
>
> Cheers and kind regards,
> Bjoern
>
> --
> Bjoern Olausson
> Martin-Luther-Universität Halle-Wittenberg
> Fachbereich Biochemie/Biotechnologie
> Kurt-Mothes-Str. 3
> 06120 Halle/Saale
>
> Phone: +49-345-55-24942
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:56:36 CST