Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

How can one import the coordinates stored in a binary CHARMM/NAMD DCD file?

This is the structure of the file:

HDR     NSET    ISTRT   NSAVC   5-ZEROS NATOM-NFREAT    DELTA   9-ZEROS
`CORD'  #files  step 1  step    zeroes  (zero)          timestep  (zeroes)
                        interval
C*4     INT     INT     INT     5INT    INT             DOUBLE  9INT
==========================================================================
NTITLE          TITLE
INT (=2)        C*MAXTITL
                (=32)
==========================================================================
NATOM
#atoms
INT
==========================================================================
X(I), I=1,NATOM         (DOUBLE)
Y(I), I=1,NATOM         
Z(I), I=1,NATOM         
==========================================================================

A sample pdb and dcd file is here.

A C plugin for reading the DCD can be found here. As with any file format without a rigid formal specification there is a lot of variation and different edge cases...


Surprisingly googling did not turn up any ready made recipes. I'm hoping somebody has something on the shelf:)

Of course this is trivial to implement using BinaryReadList, so if nobody else answers, I will answer this question myself, there by still making it useful to the community.

share|improve this question
    
Is the dcd format related to the pdb format? –  rcollyer Apr 10 '12 at 15:06
    
Well, normally one needs a combination of a protein structure file - psf (lrz.de/~heller/ego/manual/node88.html) or pdb file, which contains the atom types and connectivity and a dcd file which contains the coordinates at different times (the trajectory). But the internal structure of the two formats is completely different. –  Ajasja Apr 10 '12 at 17:50
add comment

1 Answer 1

up vote 4 down vote accepted

I ported most of the matlab package to Mathematica. Here is the result

(*Some utility functions*)
SetAttributes[MapShowIt, {HoldAll, Listable}];
MapShowIt[code__] := MapShowIt[{code}];
MapShowIt[code_] := With[{y = code}, Defer[code = y]]

ClearAll[Puts]
SetAttributes[Puts, HoldAll];
Options[Puts] = {DisplayFunction -> Shallow};
$VerbosePrint = False;
Puts[msg_, data_, opt : OptionsPattern[]] := 
  If[$VerbosePrint, Print[msg, OptionValue[DisplayFunction][data]]];
Puts[msg_List] := If[$VerbosePrint, Print[MapShowIt[msg]]];
Puts[msg_] := If[$VerbosePrint, Print[msg]];
(*warning this is still unsafe ... it should be done as \
answered  herehttp://mathematica.stackexchange.com/a/4189/745*)

Unprotect[Dot];
SetAttributes[Dot, HoldRest];
Dot[h_, key_Symbol | key_String] := 
  System`Utilities`HashTableGet[h, ToString[Unevaluated[key]]];
Dot /: Set[Dot[h_, key_Symbol | key_String], value_] := (
   Quiet[
      System`Utilities`HashTableRemove[h, ToString[Unevaluated[key]]]
    , System`Utilities`HashTableRemove::norem];
   System`Utilities`HashTableAdd[h, ToString[Unevaluated[key]], 
    value]);
Protect[Dot];
On[Assert];

ClearAll[readDCDHeader];
Options[readDCDHeader] = {"Verbose" -> False};
readDCDHeader::nonintframes = 
  "Number of frames calculated from files size (`1`) is not an \
integer!";
readDCDHeader::diffframes = 
  "Header claims `1` frames, but there are `2` frames!";
readDCDHeader[fileName_, opts : OptionsPattern[]] :=
  Module[{str, magicnum, h, i, newsize, newsize1, numlines, fint}, 
   Block[{$VerbosePrint = OptionValue["Verbose"]},
      h = System`Utilities`HashTable[];

      str = OpenRead[fileName, BinaryFormat -> True];


      fint = BinaryRead[str, "Integer32"];
      (*if we don't read an 84 then try to reverse the endianes
      If[fint=!=84,(*then*)
        Close[str];
       str= OpenRead[fileName,BinaryFormat -> True,
    ByteOrdering->-$ByteOrdering];]*)
      Assert[fint == 84, "First integer must be 84"];
      h.stream = str;
      magicnum = BinaryReadList[str, "Character8", 4];

      Assert[magicnum == {"C", "O", "R", "D"}, "CORD not present"];

      h.nset = BinaryRead[str, "Integer32"];
      h.istart = BinaryRead[str, "Integer32"];
      h.nsavc = BinaryRead[str, "Integer32"];

      (*read free indexes*)
      SetStreamPosition[str, 40];
      h.numfree = BinaryRead[str, "Integer32"];

      Puts[{h.nset, h.istart, h.nsavc, h.numfree}];

      (*find out if is charmm DCD*)
      SetStreamPosition[str, 84];
      i = BinaryRead[str, "Integer32"];
      If[i == 0, (*then*)
         h.charmm = False;
       ,(*else*)
       h.charmm = True;

       (* check for extra block*)
       SetStreamPosition[str, 48];
       i = BinaryRead[str, "Integer32"];
       h.charmm$extrablock = (i == 1);

                     SetStreamPosition[str, 52];
                      i = BinaryRead[str, "Integer32"];
       h.charmm$4dims = (i == 1);
       ];

      Puts[{h.charmm, h.charmm$extrablock, h.charmm$4dims}];

      (*read the timestep*)  
      SetStreamPosition[str, 44];
      If[h.charmm, (*then*)
       h.DELTA = BinaryRead[str, "Real32"];
       ,(*else*)
       h.DELTA = BinaryRead[str, "Real64"];];
      h.step = h.DELTA;

      (*get the title*)
      SetStreamPosition[str, 92];
      newsize = BinaryRead[str, "Integer32"];
      numlines = BinaryRead[str, "Integer32"];
      (*TODO check for curoupted Ntitle values*)
      h.title = 
     StringJoin@BinaryReadList[str, "Character8", numlines*80];
      newsize1 = BinaryRead[str, "Integer32"];
      Assert[newsize == newsize1];
      Puts[h.title];
      i = BinaryRead[str, "Integer32"]; 
    Assert[i == 4, "4 must be read before num of atoms"];
      h.numatoms = BinaryRead[str, "Integer32"];
      h.N = h.numatoms;
      i = BinaryRead[str, "Integer32"]; 
    Assert[i == 4, "4 must be read after num of atoms"];
      Puts[{h.DELTA, h.N}];

      (*love this comment from the original matdcd package*)
      (*stuff with freeindexes.  Just smile and nod.*)
      If[ h.numfree =!= 0, (*then*)
          i = BinaryRead[str, "Integer32"];  (* should be N-
     NAMNF*4*)

     h.freeindexes = BinaryReadList[str, "Integer32", h.N - h.numfree];
          i = BinaryRead[str, "Integer32"];  (* should be N-
     NAMNF*4*)
       ];

    h.headerend = StreamPosition[str];  
    (*calculate one frame size in bytes*)
    h.framesize = 3*(h.numatoms*4 + 2*4(*for the blocksize*))
        + If[h.charmm$extrablock, 12*4 + 2*4, 0]
        + If[h.charmm$4dims, +h.numatoms*4 + 4*2, 0];

    h.numframes = (FileByteCount[fileName] - h.headerend)/h.framesize;
    (* Warn if noninteger frame number and if the actuaal frames \
differ from h.nset*)
    If[Head[h.numframes] =!= Integer,(*then*)
       Message[readDCDHeader::nonintframes, h.numframes]];
    If[Head[h.numframes] != h.nset,(*then*)
       Message[readDCDHeader::diffframes, h.nset, h.numframes]];

    Return[h];
    ]];


ClearAll[readDCDStep];
Options[readDCDStep] = {"Verbose" -> False,

   "Atoms" -> All(*or a one based list of atoms to take*)};
readDCDStep[h_System`Utilities`HashTable, opts : OptionsPattern[]] :=


  Module[{x, y, z, str, blocksize, ind}, 
   Block[{$VerbosePrint = OptionValue["Verbose"]},
      ind = OptionValue["Atoms"];
      Assert[h.numfree == 0, 
     "Fixed atoms anad free indices are not supported"];
      str = h.stream;

      If[h.charmm && h.charmm$extrablock, (*then*)
       (*unit cell info*) 
         blocksize = BinaryRead[str, "Integer32"];
         Puts[
      "Skipping unit info cords. (blocksize: " <> 
       ToString[blocksize] <> ")"];
         Skip[str, "Byte", blocksize];

         Assert[blocksize == BinaryRead[str, "Integer32"], 
      "Wrong blocksize in extra block "];
       ];
     (* Get x coordinates *)
      blocksize = BinaryRead[str, "Integer32"];
      Puts[
     "Getting x cords. (blocksize: " <> ToString[blocksize] <> ")"];
      x = BinaryReadList[str, "Real32", blocksize/4]; 
      If[Head[ind] == List, x = Part[x, ind]];
      Puts["x:\n", x];
      Assert[blocksize == BinaryRead[str, "Integer32"], 
     "Wrong blocksize in x coords"];

      (* Get y coordinates *)
      blocksize = BinaryRead[str, "Integer32"];
      Puts[
     "Getting y cords. (blocksize: " <> ToString[blocksize] <> ")"];
      y = BinaryReadList[str, "Real32", blocksize/4];
      If[Head[ind] == List, y = Part[y, ind]];
      Puts["y:\n", y];
      Assert[blocksize == BinaryRead[str, "Integer32"], 
     "Wrong blocksize in y coords"];

      (* Get z coordinates *)
      Puts[
     "Getting z cords. (blocksize: " <> ToString[blocksize] <> ")"];
      blocksize = BinaryRead[str, "Integer32"];
      z = BinaryReadList[str, "Real32", blocksize/4];
      If[Head[ind] == List, z = Part[z, ind]];
      Puts["z:\n", z];
      Assert[blocksize == BinaryRead[str, "Integer32"], 
     "Wrong blocksize in z coords"];

      (*skip 4th dimension if it exists*)
      If[h.charmm && h.charmm$4dims,(*then*)
         Puts["Skipping w cords."];   
         blocksize = BinaryRead[str, "Integer32"];
                       Skip[str, "Byte", blocksize];

     Assert[blocksize == BinaryRead[str, "Integer32"], 
      "Wrong blocksize in 4th dim"];
       ];
    Assert[Length[x] == Length[y] == Length[z], 
     "Wrong size of x or y or z"];  
    Return[Developer`ToPackedArray[Transpose@{x, y, z}]];
    ]];

ClearAll[CloseDCD];
CloseDCD[h_System`Utilities`HashTable] := Close[h.stream];

ClearAll[ImportDCD];
Options[ImportDCD] = Evaluate[Options[readDCDStep]];
ImportDCD[fileName_, options : OptionsPattern[]] :=
  Module[{h, data, ropts, opts}, 
   Block[{$VerbosePrint = OptionValue["Verbose"]},
      opts = 
     DeleteDuplicates[{options}~Join~Options[ImportDCD], 
      ToString[First[#1]] == ToString[First[#2]] &];
      Puts[opts];  
      h = readDCDHeader[fileName, Verbose -> OptionValue["Verbose"]];

      ropts = Evaluate[FilterRules[opts, Options[readDCDStep]]];
      data = (readDCDStep[h, ropts])
                  & /@ Range[h.numframes];
      CloseDCD[h];

      Return[data]
    ]];

Some example usage:

  • Read the header

    h = readDCDHeader["wat.dcd", Verbose -> True]; CloseDCD[h];

  • Read one step

    h = readDCDHeader["wat.dcd", Verbose -> True]; data = readDCDStep[h, Verbose -> True]; Short[data] CloseDCD[h];

  • Import the whole trajectory

    data = ImportDCD["wat.dcd", Verbose -> False]; Dimensions[data] (*and for fun plot it:*)

    Manipulate[ pdata = Sphere[#] & /@ data[[i]]; Graphics3D[pdata] , {i, 1, Length[data], 1}]

Mathematica graphics

Its also possible to specify a list of one based indexes of which atoms to load using the "Atoms" option.

I really should put this on github and convert it into a package...

share|improve this answer
1  
Have you thought about adapting it to the ImportConverter paradigm? With the functionality written, it shouldn't be terribly difficult to do. –  rcollyer Apr 26 '12 at 15:30
    
A good idea and it should be fairly trivial to implement. –  Ajasja Apr 27 '12 at 8:54
    
You should host this at a more suitable place, such as GitHub. Very useful package, works as advertised! –  Szabolcs Dec 2 '13 at 23:40
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.