6
$\begingroup$

I am trying to read a Nanoscale Molecular Dynamics (NAMD), Dihedral Coordinates (DCD) file using Python. Following the code given at here to read DCD files using Fortran, I have written the following code snippet to read the header section of the DCD file.

import struct
import typing


class Trajectory:
    def __init__(self, psfname: str, dcdfile: str) -> None:
        self.psfName = psfname
        self.dcdFile = dcdfile
        self.nAtoms = 0
        self.nFrames = 0
        self.readDCD()

    def readDCD(self) -> None:
        with open(self.dcdFile, "rb") as dcd:
            header = dcd.read(92)
            test = struct.unpack("4s", header[0:4])[0]
            print(header)


trajectory = Trajectory(
    "/home/ssm_grp/project2/jpcb.paper/stripped/additive/adenine/stripped.ade.psf",
    "/home/ssm_grp/project2/jpcb.paper/stripped/additive/adenine/stripped.ade.full.add.dcd",
)

The Fotran code gave the value 'CORD' which is expected, while the python code gave b'T\x00\x00\x00'

Edit: I figured that this bitstring b'T\x00\x00\x00' corresponds to the value 84, indicating that the DCD file is in little-endian format. Now I am confused as to why is that four bits not read in Fortran. Any reference on why this is happening, or how to read DCD files in python would be appreciated.

$\endgroup$
9
  • 1
    $\begingroup$ Why do you not want to use third party packages? MDAnalysis is quite good, and cross-platform and it can work with DCD files. $\endgroup$
    – S R Maiti
    Commented Mar 31, 2023 at 13:00
  • $\begingroup$ @SRMaiti, I have been working with MDAnalysis codes for my analysis scripts. This is just a side project for me. $\endgroup$ Commented Mar 31, 2023 at 14:59
  • $\begingroup$ @HemanthHaridas could you please add the input DCD file? And what is the correct output you expect ? Also please add the full forms of the words NAMD and DCD $\endgroup$ Commented Mar 31, 2023 at 20:38
  • 1
    $\begingroup$ You can check the source code for chemfiles: github.com/chemfiles/chemfiles/blob/master/src/formats/DCD.cpp, they implement DCD reader in C++, you maybe able to translate that into python. (Fortran is much less understandable than C++ imo) $\endgroup$
    – S R Maiti
    Commented Mar 31, 2023 at 20:39
  • $\begingroup$ @VandanRevanur, I would not be able to share the DCD file because it is from an unpublished work. NAMD stands for Nanoscale Molecular Dynamics and DCD stands for Dihedral Coordinates $\endgroup$ Commented Apr 1, 2023 at 3:07

1 Answer 1

3
$\begingroup$

With some digging around, I found this webpage which gave some hints to the specific padding required to read the various data from the DCD header. Even though the data given in the above mentioned webpage was not completely correct, I could just print out the specific bits to figure out the rest. Given below is a working code that can read the header section of the DCD file.

import struct
import typing
import numpy 
import os

class Trajectory:
    def __init__(self, psf: str, dcd: str) -> None:
        self.psf = psf
        self.dcd = dcd
        self.__readDCD__()

    def __readDCD__(self):
        with open(self.dcd, "rb") as DCDReader:
            header = DCDReader.read(272)
            endian = struct.unpack("i", header[0:4])[0]                 # Must be equal to 84
            cordst = struct.unpack("4s", header[4:8])[0]                # CORD
            
            # Trajectory information
            nframe = struct.unpack("i", header[8:12])[0]                # Number of frames
            iframe = struct.unpack("i", header[12:16])[0]               # First frame in DCD file
            stride = struct.unpack("i", header[16:20])[0]               # Stride between two consecutive frames
            lframe = struct.unpack("i", header[20:24])[0]               # Last frame in DCD file
            nframe = int((lframe - iframe) / stride)

            # List of five zeros => Does not do anything
            fzero1 = struct.unpack("i", header[24:28])[0]
            fzero2 = struct.unpack("i", header[28:32])[0]
            fzero3 = struct.unpack("i", header[32:36])[0]
            fzero4 = struct.unpack("i", header[36:40])[0]
            fzero5 = struct.unpack("i", header[40:44])[0]

            timest = struct.unpack("f", header[44:48])[0]               # Timestep
            testst = struct.unpack("i", header[48:52])[0]               # 

            nzero1 = struct.unpack("i", header[52:56])[0]
            nzero2 = struct.unpack("i", header[56:60])[0]
            nzero3 = struct.unpack("i", header[60:64])[0]
            nzero4 = struct.unpack("i", header[64:68])[0]
            nzero5 = struct.unpack("i", header[68:72])[0]
            nzero6 = struct.unpack("i", header[72:76])[0]
            nzero7 = struct.unpack("i", header[76:80])[0]
            nzero8 = struct.unpack("i", header[80:84])[0]
            nzero9 = struct.unpack("i", header[84:88])[0]
           
            ntitle = struct.unpack("i", header[96:100])[0]              # number of title strings
            tsizes = 100 + ntitle*80                                    # size of title strings            
            tbsize = struct.unpack("i", header[tsizes:tsizes+4])[0]     # size of title block
            constv = struct.unpack("i", header[tsizes+4:tsizes+8])[0]   # value 4
            natoms = struct.unpack("i", header[tsizes+8:tsizes+12])[0]  # Number of atoms
            # Read the first 272 bytes of data and parse the header information

            # Just a hack to find the number of remaining bytes
            DCDReader.seek(0, os.SEEK_END)
            remainingSize   = DCDReader.tell() - 272
            DCDReader.seek(272, 0)
            
            firstframe  = DCDReader.read(60)
            testvl  = struct.unpack("i", firstframe[4:8])[0]    # size of pbc length
            pbcvecx = struct.unpack("d", firstframe[8:16])[0]   # PBC length in x
            pbcvecy = struct.unpack("d", firstframe[24:32])[0] # PBC length in y
            pbcvecz = struct.unpack("d", firstframe[48:56])[0] # PBC length in z
            gamma   = struct.unpack("d", firstframe[16:24])[0]  # gamma
            beta    = struct.unpack("d", firstframe[40:48])[0]  # beta
            alpha   = struct.unpack("d", firstframe[32:40])[0]  # alpha
            print(testvl, pbcvecx, pbcvecy, pbcvecz, gamma, beta, alpha)
#            print(consta, testvl, pbcvecx, pbcvecy, pbcvecz, numpy.arccos(alpha)*180/numpy.pi, numpy.arccos(beta)*180/numpy.pi, numpy.arccos(gamma)*180/numpy.pi)
```
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .