Re: Get number of timesteps from DCD via script

From: Bjoern Olausson (namdlist_at_googlemail.com)
Date: Wed Feb 09 2011 - 06:37:20 CST

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 : Wed Feb 29 2012 - 15:56:37 CST