[gmx-developers] How to get the number of frames contained by an .xtc trajectory file??

[Mirco Wahab]

Hi Tsjerk,

> If you have a small C program for calculating the number of
> frames, please do post it. It might be interesting for others.

I tried a crude, but self-contained tool some time ago that never
really worked, but after your nice posting (python script) on this
topic I was able to make it work reliably. Maybe it's of use for something.

Regards and thanks for your hints,

M.

==> xtcplay.cxx
-------------- next part --------------
/*
 [2012-06-06]
  == Version 0.01
  = XTCPLAY - play around w/Gromacs xtc trajectory files
  = self-contained, no xtc lib - works on LITTLE ENDIAN machines (Intel, AMD)
  = xtc_brute_force_scan() -> prints file offsets of every frame in the xtc
  = xtc_dump_last_frame() -> displays various xtc info and writes last frame
  = written after reading a usenet posting from Tsjerk Wassenaar 
  = C++ with minimal error handling, just for demonstration
  = tested with MSVC11/Win64 and gcc 4.3.4/x64 ( g++ xtcplay.cxx -o xtcplay)
  = please send error info or improvements to mirco.wahab at chemie.tu-freiberg.de
*/

// if Intel, AMD etc., see 
// http://en.wikipedia.org/wiki/Endianness#Endianness_and_operating_systems_on_architectures
#define _LSB_FIRST_ 1  
// otherwise:
// #define _MSB_FIRST_ 1

#include <stdint.h>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <string>
using namespace std;

int xtc_brute_force_scan(const string& fname);
int xtc_dump_last_frame(const string& fname);
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  //

 int main(int argc, char* argv[])
{
 string fname(argv[1]);

 // int status = xtc_brute_force_scan(fname);
 int status = xtc_dump_last_frame(fname);

 return status;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  //

 uint32_t u4_from_buffer(char b[]) 
{
#if defined(_LSB_FIRST_)
 return uint8_t(b[3]) <<  0 | uint8_t(b[2]) <<  8
      | uint8_t(b[1]) << 16 | uint8_t(b[0]) << 24  ;
#else
 return * (uint32_t*) b;  // is this correct? I don't have a MSB first machine
#endif
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  //

 int64_t find_xtc_frame_tag(char *block, int64_t len1, char *tag)
{
 char *p=block, *q=block+len1, *t=tag;
 while(p < q) { // unrolled for speed
    if(   *(p+0) == *(t+0) && *(p+1) == *(t+1) && *(p+2) == *(t+2)
       && *(p+3) == *(t+3) && *(p+4) == *(t+4) && *(p+5) == *(t+5)
       && *(p+6) == *(t+6) && *(p+7) == *(t+7)  )  
       return p-block;
    ++p;
 }
 return -1;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  //

 int xtc_brute_force_scan(const string& fn)
{
 // read the whole xtc to memory - can your box handle this?
 ifstream fd;
 fd.open(fn.c_str(), ios::binary);
  if(!fd.is_open()) {
    perror(fn.c_str());
    exit(1);
 }

 fd.seekg(0, ios_base::end);       
 streamsize flen = fd.tellg();
 fd.seekg(0, ios_base::beg);       

 cout << "FILE SIZE = " << flen/(1024*1024) << " MB\nFRAME\tOFSETT" << endl;
 char *buf = new char[size_t(flen)];
 if(buf == NULL) {
    cerr << "out of memory" << endl;
    exit(1);
 }

 fd.read(buf, flen);
 streamsize rdlen = fd.gcount();
 if(rdlen) {
    streamsize cnt=0, pos=0, offs=0;
    while( (pos=find_xtc_frame_tag(buf+offs, rdlen-offs, buf)) != -1 ) {
       cnt++;
       cout << dec << cnt << "\t" << (pos+offs) << "/" << flen 
            << "\t(0x" << hex << (pos+offs) << ")" << endl; 
       offs += pos+92;
    } 
    cout << "TOTAL = " << dec << cnt << " FRAMES" << endl;
 }
 else {
    cerr << "short read" << endl;
    exit(1);
 }
 
 delete [] buf;
 fd.close();
 return 0;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  //

 int xtc_dump_last_frame(const string& fn)
{
 ifstream fd;
 fd.open(fn.c_str(), ios::binary);
 if(!fd.is_open()) {
    perror(fn.c_str());
    exit(1);
 }
 
 char header[92];
 fd.read(header, 92);
 uint32_t magic = u4_from_buffer(header);           // 1995
 uint32_t natom = u4_from_buffer(header+4);         // natoms
 uint32_t frmsz = 92 + u4_from_buffer(header+88);   // framesize
 cout << "FILE=" << fn.c_str() << ", MAGIC=" << magic 
      << ", NATOMS=" << natom << ", FRAMESIZE=" << frmsz << " BYTES, " ;

 streamsize offs = frmsz + (frmsz/4); // guess last frame position 
 fd.seekg(-offs, ios_base::end);      // move read pointer up there
 streamsize fpos = fd.tellg();
 streamsize nframes = fpos/frmsz+1;
 cout << "NFRAMES=" << nframes << endl;
 
 char *rdbuf = new char [size_t(offs+1)]; // prepare load of last frame
 fd.read(rdbuf, offs+1);                  // read last block from file
 streamsize rdlen = fd.gcount();
 fd.close();
                                          // find frame start
 int64_t pos = find_xtc_frame_tag(rdbuf, rdlen, header); 
 if(pos == -1) {
    cerr << "error: no frame header found for last frame" << endl;
    return 1;
 }

 cout << "write last frame of " << nframes << " (" << (rdlen-pos) 
      << " bytes)" << endl;
 string fname = "lastframe_" + fn;
 ofstream fh(fname.c_str(), ios::binary);
 fh.write(rdbuf+pos, rdlen-pos);
 fh.close();
 delete [] rdbuf;

 return 0;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  //