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

[Tsjerk Wassenaar]

Hey Mirco,

That's pretty cool. I had a python script for finding all frame
positions to write a reversed trajectory. Very handy with trjconv -pbc
nojump to make a trajectory smoothly converge to the end point. I
avoided reading in the whole trajectory. Maybe it's of some use to you
(it even has comments).

Cheers,

Tsjerk

###

#!/usr/bin/env python

# Reverse an XTC trajectory
#
# (c)2011 Tsjerk A. Wassenaar
#   University of Groningen
#

from struct import unpack
import sys
import os

def i(x): return sum([ord(x[j])<<(24-j*8) for j in range(4)])

def strseek(stream,s,bufsize=10000):
    v = len(s)
    x = stream.read(bufsize)
    n = 0
    while len(x) >= v:
        m = x.find(s)
        if m > -1:
            # If we found the tag, the position is the total length
            # read plus the index at which the tag was found
            n += m
            yield n
            # Now we strip the part up to and including the tag
            x = x[m+v:]
            n += v
        elif len(x) > v:
            # If the tag was not found, strip the string, keeping only
            # the last v-1 characters, as these could combine with the
            # next part to be read in.
            # n is incremented with the number of characters taken from
            # the current string x (len(x)-v+1)
            n += len(x)-v+1
            x = x[1-v:]
        if len(x) <= v:
            x += stream.read(bufsize)

# Get the tag to identify the start of a frame
f   = open(sys.argv[1])
tag = f.read(8)                     # Tag: magic number and number of atoms
n   = 92 + i(f.read(84)[-4:])       # Size of frame in bytes
f.close()

# Find all positions at which the tag is found and the frame lengths
frames  = [ i for i in strseek(open(sys.argv[1]),tag) ]
nf      = len(frames)
lengths = [ j-i for i,j in zip(frames,frames[1:]+[nf]) ]

# Reverse the lists
frames.reverse()
lengths.reverse()

# Write the frames in reversed order
f = open(sys.argv[1])
o = open(sys.argv[2],"w")
for pos,bytes in zip(frames,lengths):
    f.seek(pos)
    o.write(f.read(bytes))

# DONE


On Wed, Jun 6, 2012 at 12:49 AM, Mirco Wahab
<mirco.wahab at chemie.tu-freiberg.de> wrote:
> 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
>
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.



-- 
Tsjerk A. Wassenaar, Ph.D.

post-doctoral researcher
Molecular Dynamics Group
* Groningen Institute for Biomolecular Research and Biotechnology
* Zernike Institute for Advanced Materials
University of Groningen
The Netherlands