Re: Get number of timesteps from DCD via script

From: Joshua D. Moore (joshuadmoore_at_gmail.com)
Date: Wed Feb 09 2011 - 06:41:06 CST

Just paste what I sent into a file named "read_dcd_header.f90"

Then compile with

f90 -o read_dcd_header.exe read_dcd_header.f90

Then execute with

./read_dcd_header.exe

When it asks you for the name of your dcd file, type in your dcd filename.

On Wed, Feb 9, 2011 at 7:37 AM, Bjoern Olausson <namdlist_at_googlemail.com>wrote:

> On Tuesday 08 February 2011 23:56:02 you wrote:
> > 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
> >
>
> Mhhh, as I sayed, I have noe clue abotu Fortran and it is a hard time
> reading
> the code and editing it. I understand the above code the following way:
>
> read(243) HDR,(ICNTRL(j),j=1,9),delta,(ICNTRL(j),j=11,20)
>
> Read 4 Chars --> COOR
> Read 1 int(4) 1--> # of configurations
> Read 1 int(4) 2--> time step of first configuration
> Read 1 int(4) 3--> frequency of configurations
> Read 1 int(4) 4--> last time step in dcd file
> Read 1 int(4) 5--> 0
> Read 1 int(4) 6--> 0
> Read 1 int(4) 7--> 0
> Read 1 int(4) 8--> 0
> Read 1 int(4) 9--> 0
> Read 1 real(4)1--> timestep of simulation
> Read 1 int(4) 11--> 0
> Read 1 int(4) 12--> 0
> Read 1 int(4) 13--> 0
> Read 1 int(4) 14--> 0
> Read 1 int(4) 15--> 0
> Read 1 int(4) 16--> 0
> Read 1 int(4) 17--> 0
> Read 1 int(4) 18--> 0
> Read 1 int(4) 19--> 0
> Read 1 int(4) 20--> 0
> READ 1 int(4) 1--> Length of title
> READ 1 char(80) --> title
> READ 1 int(4) --> Number of Atoms
>
>
> Where is the PBC here?
>
> This is, at least in my opinion, exactly what I am doing in my python
> script.
> Except that I am not sure how to treat the "real" type i python. They types
> available are listed here:
> http://docs.python.org/library/struct.html#format-characters
>
> Sorry to bug you with this, but it is frustrating without a proper
> documentation and not being able to easily read Fortran.
>
> Could anyone complement/correct the header structure:
> HDR NSET ISTART NSAVC ILAST ZEROS5 DELTA ZEROS9
> i4c i i i i 5i real 9i
>
> After the DELTA my script breaks and I am not sure if it id because of a
> wrong
> type or some other missing value.
>
> http://paste.pocoo.org/show/334948/
>
> #!/bin/env python
> # -*- coding: utf-8 -*-
> #
> #
>
> import sys, struct
>
> #Info: SIMULATION PARAMETERS:
> #Info: TIMESTEP 2
> #Info: NUMBER OF STEPS 1500000
> #Info: STEPS PER CYCLE 10
> #Info: PERIODIC CELL BASIS 1 93.1932 0 0
> #Info: PERIODIC CELL BASIS 2 0 93.1932 0
> #Info: PERIODIC CELL BASIS 3 0 0 143.704
> #Info: PERIODIC CELL CENTER 0 0 0
>
> 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]
> ILAST = struct.unpack('i', f.read(struct.calcsize('i')))
> ZEROS5 = struct.unpack('5i', f.read(struct.calcsize('5i')))
> DELTA = struct.unpack('f',f.read(struct.calcsize('f')))
> ZEROS9 = struct.unpack('9i', f.read(struct.calcsize('9i')))
> NTITL = struct.unpack('i',f.read(struct.calcsize('i')))
> NTITL = struct.unpack('d',f.read(struct.calcsize('d')))
> TITLE = struct.unpack('i200c',f.read(struct.calcsize('i200c')))
> NATOM = struct.unpack('i',f.read(struct.calcsize('i')))
>
> ISTOP = ISTART + (NSET * NSAVC)
> NSTP = ISTOP - ISTART
>
> print "%(ID)-13s %(DATA)s" %{"ID": "HDR:", "DATA": HDR}
> print "%(ID)-13s %(DATA)s" %{"ID": "Coor Sets:", "DATA": NSET}
> print "%(ID)-13s %(DATA)s" %{"ID": "First Step:", "DATA": ISTART}
> print "%(ID)-13s %(DATA)s" %{"ID": "Save Int.:", "DATA": NSAVC}
> print "%(ID)-13s %(DATA)s" %{"ID": "Last Step:", "DATA": ILAST}
> print "%(ID)-13s %(DATA)s" %{"ID": "5-ZEROS:", "DATA": ZEROS5}
> print "%(ID)-13s %(DATA)s" %{"ID": "Int. Time:", "DATA": DELTA}
> print "%(ID)-13s %(DATA)s" %{"ID": "9-ZEROS:", "DATA": ZEROS9}
> print "%(ID)-13s %(DATA)s" %{"ID": "Title chars:", "DATA": NTITL}
> print "%(ID)-13s %(DATA)s" %{"ID": "Title:", "DATA": TITLE}
> print "%(ID)-13s %(DATA)s" %{"ID": "Atoms:", "DATA": NATOM}
>
> print "\n%(ID)-13s %(DATA)i" %{"ID": "Last Step:", "DATA": ISTOP}
> print "%(ID)-13s %(DATA)i" %{"ID": "Num Steps:", "DATA": NSTP}
>
>
> And the output:
> blub_at_repulsion $ python dcd-header.py
> ../bdvm/namd/membrane30_400-bdvm-rna-
> smd/run002/memb30_run400-bdvm-rna_npgt.dcd
>
> HDR: (84, 'C', 'O', 'R', 'D')
> Coor Sets: 1000
> First Step: 1500500
> Save Int.: 500
> Last Step: (2000000,)
> 5-ZEROS: (0, 0, 0, 0, 0)
> Int. Time: (0.040909659117460251,)
> 9-ZEROS: (1, 0, 0, 0, 0, 0, 0, 0, 0)
> Title chars: (3.4800730975980617e-312,)
> Title: (2, 'R', 'E', 'M', 'A', 'R', 'K', 'S', ' ', 'F', 'I', 'L',
> 'E',
> 'N', 'A', 'M', 'E', '=', 'm', 'e', 'm', 'b', '3', '0', '_', 'r', 'u', 'n',
> '4', '0', '0', '-', 'b', 'd', 'v', 'm', '-', 'r', 'n', 'a', '_', 'n', 'p',
> 'g', 't', '.', 'd', 'c', 'd', ' ', 'C', 'R', 'E', 'A', 'T', 'E', 'D', ' ',
> 'B', 'Y', ' ', 'N', 'A', 'M', 'D', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',
> '
> ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', 'R', 'E', 'M', 'A', 'R', 'K', 'S', '
> ',
> 'D', 'A', 'T', 'E', ':', ' ', '0', '7', '/', '0', '5', '/', '1', '0', ' ',
> 'C', 'R', 'E', 'A', 'T', 'E', 'D', ' ', 'B', 'Y', ' ', 'U', 'S', 'E', 'R',
> ':', ' ', 'b', 'l', 'u', 'b', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',
> '
> ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', '
> ',
> ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', '\xa4', '\x00',
> '\x00',
> '\x00', '\x04', '\x00', '\x00', '\x00', '\xa9', '\xfb', '\x01', '\x00',
> '\x04', '\x00', '\x00', '\x00', '0', '\x00', '\x00', '\x00', '\xfe', '{',
> '<',
> '\xa1', 'E', 'N', 'W', '@', '\x00', '\x00', '\x00', '\x00', '\x00', '\x00',
> '\x00', '\x00', '\xfe', '{', '<', '\xa1')
> Atoms: (1079463493,)
>
> Last Step: 2000500
> Num Steps: 500000
>
>
>
> --
> 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 : Mon Dec 31 2012 - 23:19:47 CST