From bb13d32a02d7440909de2261bada82525b5ca9e1 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 30 Jun 2016 21:45:36 +0200 Subject: [PATCH 01/26] first dcd stumb cython file --- .../MDAnalysis/lib/formats/include/correl.h | 188 +++++ .../lib/formats/include/endianswap.h | 173 ++++ .../MDAnalysis/lib/formats/include/fastio.h | 454 +++++++++++ .../MDAnalysis/lib/formats/include/readdcd.h | 764 ++++++++++++++++++ package/MDAnalysis/lib/formats/libdcd.pyx | 59 ++ 5 files changed, 1638 insertions(+) create mode 100644 package/MDAnalysis/lib/formats/include/correl.h create mode 100644 package/MDAnalysis/lib/formats/include/endianswap.h create mode 100644 package/MDAnalysis/lib/formats/include/fastio.h create mode 100644 package/MDAnalysis/lib/formats/include/readdcd.h create mode 100644 package/MDAnalysis/lib/formats/libdcd.pyx diff --git a/package/MDAnalysis/lib/formats/include/correl.h b/package/MDAnalysis/lib/formats/include/correl.h new file mode 100644 index 00000000000..f3dacb042e7 --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/correl.h @@ -0,0 +1,188 @@ +/* -*- Mode: C; tab-width: 4; indent-tabs-mode:nil; -*- */ +/* vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 */ +/* + MDAnalysis --- http://mdanalysis.googlecode.com + Copyright (c) 2006-2014 Naveen Michaud-Agrawal, + Elizabeth J. Denning, Oliver Beckstein, + and contributors (see AUTHORS for the full list) + Released under the GNU Public Licence, v2 or any higher version + + Please cite your use of MDAnalysis in published work: + + N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and + O. Beckstein. MDAnalysis: A Toolkit for the Analysis of + Molecular Dynamics Simulations. J. Comput. Chem. 32 (2011), 2319--2327, + in press. +*/ + +#ifndef CORREL_H +#define CORREL_H + +#include +/* Python.h for 'typedef int Py_intptr_t;' (fixes Issue 19) */ +#include + +static void +copyseries(int frame, char *data, const Py_intptr_t *strides, + const float *tempX, const float *tempY, const float *tempZ, + const char* datacode, int numdata, const int* atomlist, const int* atomcounts, + int lowerb, double* aux) +{ + char code; + int index1 = 0, index2 = 0, index3 = 0, index4 = 0; + double x1, x2, y1, y2, z1, z2, x3, y3, z3, aux1, aux2; + int i = 0, j = 0, atomno = 0; + int dataset = 0; + int stride0, stride1; + + stride0 = strides[0]; + stride1 = strides[1]; + + /* If I eventually switch to using frame,property ordering for timeseries data + stride0 = strides[1]; + stride1 = strides[0]; + */ + + for (i=0;i 0) || (aux2 > 0 && aux1 < 0)) { + aux1 *= -1; + } + // Check if the dihedral has wrapped around 2 pi + aux2 = *(double*)(data + dataset*stride0 + (frame-1)*stride1); + if (fabs(aux1-aux2) > M_PI) { + if (aux1 > 0) { aux1 -= 2*M_PI; } + else { aux1 += 2*M_PI; } + } + *(double*)(data + dataset++*stride0 + frame*stride1) = aux1; + break; + case 'w': + /* dipole orientation of 3-site water: ^ d + index1 = oxygen, index2, index3 = hydrogen | + returns d ,O, + d = rO - (rH1 + rH2)/2 H | H + | + */ + index1 = atomlist[atomno++]-lowerb; // O + index2 = atomlist[atomno++]-lowerb; // H1 + index3 = atomlist[atomno++]-lowerb; // H2 + x1 = tempX[index1] - 0.5*(tempX[index2] + tempX[index3]); // dx + y1 = tempY[index1] - 0.5*(tempY[index2] + tempY[index3]); // dy + z1 = tempZ[index1] - 0.5*(tempZ[index2] + tempZ[index3]); // dz + *(double*)(data + dataset++*stride0 + frame*stride1) = x1; + *(double*)(data + dataset++*stride0 + frame*stride1) = y1; + *(double*)(data + dataset++*stride0 + frame*stride1) = z1; + break; + } + } +} + +// This accounts for periodic boundary conditions +// taken from MMTK +#define distance_vector_2(d, r1, r2, data) \ + { \ + double xh = 0.5*(data)[0]; \ + double yh = 0.5*(data)[1]; \ + double zh = 0.5*(data)[2]; \ + d[0] = r2[0]-r1[0]; \ + if (d[0] > xh) d[0] -= (data)[0]; \ + if (d[0] <= -xh) d[0] += (data)[0]; \ + d[1] = r2[1]-r1[1]; \ + if (d[1] > yh) d[1] -= (data)[1]; \ + if (d[1] <= -yh) d[1] += (data)[1]; \ + d[2] = r2[2]-r1[2]; \ + if (d[2] > zh) d[2] -= (data)[2]; \ + if (d[2] <= -zh) d[2] += (data)[2]; \ + } + +#endif diff --git a/package/MDAnalysis/lib/formats/include/endianswap.h b/package/MDAnalysis/lib/formats/include/endianswap.h new file mode 100644 index 00000000000..6f1ced37483 --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/endianswap.h @@ -0,0 +1,173 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2003 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: endianswap.h,v $ + * $Author: eamon $ $Locker: $ $State: Exp $ + * $Revision: 1.3 $ $Date: 2004/04/16 15:37:00 $ + * + *************************************************************************** + * DESCRIPTION: + * Byte swapping routines used in various plugins + * There are two versions of each routine, one that's safe to use in + * all cases (but is slow) and one that is only safe to use on memory + * addresses that are aligned to the word size that's being byte-swapped + * but are much much much faster. Use the aligned versions of these + * routines whenever possible. The 'ndata' length count parameters and + * internal loops should be safe to use on huge memory arrays on 64-bit + * machines. + * + ***************************************************************************/ + +#ifndef ENDIAN_SWAP_H +#define ENDIAN_SWAP_H + +/* works on unaligned 2-byte quantities */ +static void swap2_unaligned(void *v, long ndata) { + long i; + char * dataptr = (char *) v; + char tmp; + + for (i = 0; i < ndata-1; i += 2) { + tmp = dataptr[i]; + dataptr[i] = dataptr[i+1]; + dataptr[i+1] = tmp; + } +} + + +/* works on unaligned 4-byte quantities */ +static void swap4_unaligned(void *v, long ndata) { + long i; + char *dataptr; + char tmp; + + dataptr = (char *) v; + for (i=0; i>8)&0xff) | ((*N&0xff)<<8)); + } +} + + +/* Only works with aligned 4-byte quantities, will cause a bus error */ +/* on some platforms if used on unaligned data. */ +static void swap4_aligned(void *v, long ndata) { + int *data = (int *) v; + long i; + int *N; + for (i=0; i>24)&0xff) | ((*N&0xff)<<24) | + ((*N>>8)&0xff00) | ((*N&0xff00)<<8)); + } +} + + +/* Only works with aligned 8-byte quantities, will cause a bus error */ +/* on some platforms if used on unaligned data. */ +static void swap8_aligned(void *v, long ndata) { + /* Use int* internally to prevent bugs caused by some compilers */ + /* and hardware that would potentially load data into an FP reg */ + /* and hose everything, such as the old "jmemcpy()" bug in NAMD */ + int *data = (int *) v; + long i; + int *N; + int t0, t1; + + for (i=0; i>24)&0xff) | ((t0&0xff)<<24) | + ((t0>>8)&0xff00) | ((t0&0xff00)<<8)); + + t1 = N[1]; + t1=(((t1>>24)&0xff) | ((t1&0xff)<<24) | + ((t1>>8)&0xff00) | ((t1&0xff00)<<8)); + + N[0] = t1; + N[1] = t0; + } +} + +#if 0 +/* Other implementations that might be faster in some cases */ + +/* swaps the endianism of an eight byte word. */ +void mdio_swap8(double *i) { + char c; + char *n; + n = (char *) i; + c = n[0]; + n[0] = n[7]; + n[7] = c; + c = n[1]; + n[1] = n[6]; + n[6] = c; + c = n[2]; + n[2] = n[5]; + n[5] = c; + c = n[3]; + n[3] = n[4]; + n[4] = c; +} + +#endif + +#endif + diff --git a/package/MDAnalysis/lib/formats/include/fastio.h b/package/MDAnalysis/lib/formats/include/fastio.h new file mode 100644 index 00000000000..998b5ebd5fb --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/fastio.h @@ -0,0 +1,454 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2009 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: fastio.h,v $ + * $Author: johns $ $Locker: $ $State: Exp $ + * $Revision: 1.20 $ $Date: 2009/04/29 15:45:29 $ + * + *************************************************************************** + * DESCRIPTION: + * This is a simple abstraction layer for system-dependent I/O calls + * that allow plugins to do binary I/O using the fastest possible method. + * + * This code is intended for use by binary trajectory reader plugins that + * work with multi-gigabyte data sets, reading only binary data. + * + ***************************************************************************/ + +/* Compiling on windows */ +#if defined(_MSC_VER) + +#if 1 +/* use native Windows I/O calls */ +#define FASTIO_NATIVEWIN32 1 + +#include +#include + +typedef HANDLE fio_fd; +typedef LONGLONG fio_size_t; +typedef void * fio_caddr_t; + +typedef struct { + fio_caddr_t iov_base; + int iov_len; +} fio_iovec; + +#define FIO_READ 0x01 +#define FIO_WRITE 0x02 + +#define FIO_SEEK_CUR FILE_CURRENT +#define FIO_SEEK_SET FILE_BEGIN +#define FIO_SEEK_END FILE_END + +static int fio_win32convertfilename(const char *filename, char *newfilename, int maxlen) { + int i; + int len=strlen(filename); + + if ((len + 1) >= maxlen) + return -1; + + for (i=0; i + +typedef FILE * fio_fd; +typedef size_t fio_size_t; /* MSVC doesn't uinversally support ssize_t */ +typedef void * fio_caddr_t; /* MSVC doesn't universally support caddr_t */ + +typedef struct { + fio_caddr_t iov_base; + int iov_len; +} fio_iovec; + +#define FIO_READ 0x01 +#define FIO_WRITE 0x02 + +#define FIO_SEEK_CUR SEEK_CUR +#define FIO_SEEK_SET SEEK_SET +#define FIO_SEEK_END SEEK_END + +static int fio_open(const char *filename, int mode, fio_fd *fd) { + char * modestr; + FILE *fp; + + if (mode == FIO_READ) + modestr = "rb"; + + if (mode == FIO_WRITE) + modestr = "wb"; + + fp = fopen(filename, modestr); + if (fp == NULL) { + return -1; + } else { + *fd = fp; + return 0; + } +} + +static int fio_fclose(fio_fd fd) { + return fclose(fd); +} + +static fio_size_t fio_fread(void *ptr, fio_size_t size, + fio_size_t nitems, fio_fd fd) { + return fread(ptr, size, nitems, fd); +} + +static fio_size_t fio_readv(fio_fd fd, const fio_iovec * iov, int iovcnt) { + int i; + fio_size_t len = 0; + + for (i=0; i +#include +#include +#include + +typedef int fio_fd; +typedef off_t fio_size_t; /* off_t is 64-bits with LFS builds */ + +/* enable use of kernel readv() if available */ +#if defined(__sun) || defined(__APPLE_CC__) || defined(__linux) +#define USE_KERNEL_READV 1 +#endif + +typedef void * fio_caddr_t; + +#if defined(USE_KERNEL_READV) +#include +typedef struct iovec fio_iovec; +#else + +typedef struct { + fio_caddr_t iov_base; + int iov_len; +} fio_iovec; +#endif + + +#define FIO_READ 0x01 +#define FIO_WRITE 0x02 + +#define FIO_SEEK_CUR SEEK_CUR +#define FIO_SEEK_SET SEEK_SET +#define FIO_SEEK_END SEEK_END + +static int fio_open(const char *filename, int mode, fio_fd *fd) { + int nfd; + int oflag = 0; + + if (mode == FIO_READ) + oflag = O_RDONLY; + + if (mode == FIO_WRITE) + oflag = O_WRONLY | O_CREAT | O_TRUNC; + + nfd = open(filename, oflag, 0666); + if (nfd < 0) { + return -1; + } else { + *fd = nfd; + return 0; + } +} + +static int fio_fclose(fio_fd fd) { + return close(fd); +} + +static fio_size_t fio_fread(void *ptr, fio_size_t size, + fio_size_t nitems, fio_fd fd) { + int i; + fio_size_t len = 0; + int cnt = 0; + + for (i=0; i= 0) + return 0; /* success (emulate behavior of fseek) */ + else + return -1; /* failure (emulate behavior of fseek) */ +} + +static fio_size_t fio_ftell(fio_fd fd) { + return lseek(fd, 0, SEEK_CUR); +} + +#endif + + +/* higher level routines that are OS independent */ + +static int fio_write_int32(fio_fd fd, int i) { + return (fio_fwrite(&i, 4, 1, fd) != 1); +} + +static int fio_read_int32(fio_fd fd, int *i) { + return (fio_fread(i, 4, 1, fd) != 1); +} + +static int fio_write_str(fio_fd fd, const char *str) { + int len = strlen(str); + return (fio_fwrite((void *) str, len, 1, fd) != 1); +} + diff --git a/package/MDAnalysis/lib/formats/include/readdcd.h b/package/MDAnalysis/lib/formats/include/readdcd.h new file mode 100644 index 00000000000..c9d0f9999c9 --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/readdcd.h @@ -0,0 +1,764 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2003 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ + +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: readdcd.h,v $ + * $Author: johns $ $Locker: $ $State: Exp $ + * $Revision: 1.32 $ $Date: 2004/09/21 20:52:37 $ + * + ***************************************************************************/ + +#ifndef READ_DCD_H +#define READ_DCD_H + +#include +#include +#include +#include +#include "endianswap.h" +#include "fastio.h" + +/* DEFINE ERROR CODES THAT MAY BE RETURNED BY DCD ROUTINES */ +#define DCD_SUCCESS 0 /* No problems */ +#define DCD_EOF -1 /* Normal EOF */ +#define DCD_DNE -2 /* DCD file does not exist */ +#define DCD_OPENFAILED -3 /* Open of DCD file failed */ +#define DCD_BADREAD -4 /* read call on DCD file failed */ +#define DCD_BADEOF -5 /* premature EOF found in DCD file */ +#define DCD_BADFORMAT -6 /* format of DCD file is wrong */ +#define DCD_FILEEXISTS -7 /* output file already exists */ +#define DCD_BADMALLOC -8 /* malloc failed */ + +/* + * Read the header information from a dcd file. + * Input: fd - a file struct opened for binary reading. + * Output: 0 on success, negative error code on failure. + * Side effects: *natoms set to number of atoms per frame + * *nsets set to number of frames in dcd file + * *istart set to starting timestep of dcd file + * *nsavc set to timesteps between dcd saves + * *delta set to value of trajectory timestep + * *nfixed set to number of fixed atoms + * *freeind may be set to heap-allocated space + * *reverse set to one if reverse-endian, zero if not. + * *charmm set to internal code for handling charmm data. + */ +static int read_dcdheader(fio_fd fd, int *natoms, int *nsets, int *istart, int *nsavc, + double *delta, int *nfixed, int **freeind, + float **fixedcoords, int *reverse, int *charmm, + char **remarks, int *len_remarks); + +/* + * Read a dcd timestep from a dcd file + * Input: fd - a file struct opened for binary reading, from which the + * header information has already been read. + * natoms, nfixed, first, *freeind, reverse, charmm - the corresponding + * items as set by read_dcdheader + * first - true if this is the first frame we are reading. + * x, y, z: space for natoms each of floats. + * unitcell - space for six floats to hold the unit cell data. + * Not set if no unit cell data is present. + * Output: 0 on success, negative error code on failure. + * Side effects: x, y, z contain the coordinates for the timestep read. + * unitcell holds unit cell data if present. + */ +static int read_dcdstep(fio_fd fd, int natoms, float *x, float *y, float *z, + float *unitcell, int nfixed, int first, int *freeind, + float *fixedcoords, int reverse, int charmm); + +/* + * Read a subset of a timestep from a dcd file + * Input: fd - a file struct opened for binary reading, from which the + * header information has already been read + * natoms, nfixed, first, *freeind, reverse, charmm - the corresponding + * items as set by read_dcdheader + * first - true if this is the first frame we are reading. + * lowerb, upperb - the range of atoms to read data for + * x, y, z: space for upperb-lowerb+1 each of floats + * unitcell - space for six floats to hold the unit cell data. + * Not set if no unit cell data is present. + * Ouput: 0 on success, negative error code on failure. + * Side effects: x, y, z contain coordinates for the range of atoms + * unitcell holds unit cell data if present. + */ +static int read_dcdsubset(fio_fd fd, int natoms, int lowerb, int upperb, float *x, float *y, float *z, + float *unitcell, int nfixed, int first, int *freeind, + float *fixedcoords, int reverse, int charmm); + +/* + * Skip past a timestep. If there are fixed atoms, this cannot be used with + * the first timestep. + * Input: fd - a file struct from which the header has already been read + * natoms - number of atoms per timestep + * nfixed - number of fixed atoms + * charmm - charmm flags as returned by read_dcdheader + * Output: 0 on success, negative error code on failure. + * Side effects: One timestep will be skipped; fd will be positioned at the + * next timestep. + */ +static int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm, int numstep); + +/* + * clean up dcd data + * Input: nfixed, freeind - elements as returned by read_dcdheader + * Output: None + * Side effects: Space pointed to by freeind is freed if necessary. + */ +static void close_dcd_read(int *freeind, float *fixedcoords); + +/* + * Write a header for a new dcd file + * Input: fd - file struct opened for binary writing + * remarks - string to be put in the remarks section of the header. + * The string will be truncated to 70 characters. + * natoms, istart, nsavc, delta - see comments in read_dcdheader + * Output: 0 on success, negative error code on failure. + * Side effects: Header information is written to the dcd file. + */ +static int write_dcdheader(fio_fd fd, const char *remarks, int natoms, + int istart, int nsavc, double delta, int with_unitcell, + int charmm); + +/* + * Write a timestep to a dcd file + * Input: fd - a file struct for which a dcd header has already been written + * curframe: Count of frames written to this file, starting with 1. + * curstep: Count of timesteps elapsed = istart + curframe * nsavc. + * natoms - number of elements in x, y, z arrays + * x, y, z: pointers to atom coordinates + * Output: 0 on success, negative error code on failure. + * Side effects: coordinates are written to the dcd file. + */ +static int write_dcdstep(fio_fd fd, int curstep, int curframe, + int natoms, const float *x, const float *y, const float *z, + const double *unitcell, int charmm); + + + +#define DCD_IS_CHARMM 0x01 +#define DCD_HAS_4DIMS 0x02 +#define DCD_HAS_EXTRA_BLOCK 0x04 + +/* READ Macro to make porting easier */ +#define READ(fd, buf, size) \ + fio_fread(((void *) buf), (size), 1, (fd)) + + +/* WRITE Macro to make porting easier */ +#define WRITE(fd, buf, size) \ + fio_fwrite(((void *) buf), (size), 1, (fd)) + +/* XXX This is broken - fread never returns -1 */ +#define CHECK_FREAD(X, msg) if (X==-1) \ + { \ + return(DCD_BADREAD); \ + } + +#define CHECK_FEOF(X, msg) if (X==0) \ + { \ + return(DCD_BADEOF); \ + } + +static int read_dcdheader(fio_fd fd, int *N, int *NSET, int *ISTART, + int *NSAVC, double *DELTA, int *NAMNF, + int **FREEINDEXES, float **fixedcoords, int *reverseEndian, + int *charmm, char **remarks, int *len_remarks) +{ + int input_integer; /* buffer space */ + int ret_val; + char hdrbuf[84]; /* char buffer used to store header */ + int NTITLE; + + /* First thing in the file should be an 84 */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading first int from dcd file"); + CHECK_FEOF(ret_val, "reading first int from dcd file"); + + /* Check magic number in file header and determine byte order*/ + if (input_integer != 84) { + /* check to see if its merely reversed endianism */ + /* rather than a totally incorrect file magic number */ + swap4_aligned(&input_integer, 1); + + if (input_integer == 84) { + *reverseEndian=1; + } else { + /* not simply reversed endianism, but something rather more evil */ + return DCD_BADFORMAT; + } + } else { + *reverseEndian=0; + } + + /* Buffer the entire header for random access */ + ret_val = READ(fd, hdrbuf, 84); + CHECK_FREAD(ret_val, "buffering header"); + CHECK_FEOF(ret_val, "buffering header"); + + /* Check for the ID string "COORD" */ + if (hdrbuf[0] != 'C' || hdrbuf[1] != 'O' || + hdrbuf[2] != 'R' || hdrbuf[3] != 'D') { + return DCD_BADFORMAT; + } + + /* CHARMm-genereate DCD files set the last integer in the */ + /* header, which is unused by X-PLOR, to its version number. */ + /* Checking if this is nonzero tells us this is a CHARMm file */ + /* and to look for other CHARMm flags. */ + if (*((int *) (hdrbuf + 80)) != 0) { + (*charmm) = DCD_IS_CHARMM; + if (*((int *) (hdrbuf + 44)) != 0) + (*charmm) |= DCD_HAS_EXTRA_BLOCK; + + if (*((int *) (hdrbuf + 48)) == 1) + (*charmm) |= DCD_HAS_4DIMS; + } else { + (*charmm) = 0; + } + + /* Store the number of sets of coordinates (NSET) */ + (*NSET) = *((int *) (hdrbuf + 4)); + if (*reverseEndian) swap4_unaligned(NSET, 1); + + /* Store ISTART, the starting timestep */ + (*ISTART) = *((int *) (hdrbuf + 8)); + if (*reverseEndian) swap4_unaligned(ISTART, 1); + + /* Store NSAVC, the number of timesteps between dcd saves */ + (*NSAVC) = *((int *) (hdrbuf + 12)); + if (*reverseEndian) swap4_unaligned(NSAVC, 1); + + /* Store NAMNF, the number of fixed atoms */ + (*NAMNF) = *((int *) (hdrbuf + 36)); + if (*reverseEndian) swap4_unaligned(NAMNF, 1); + + /* Read in the timestep, DELTA */ + /* Note: DELTA is stored as a double with X-PLOR but as a float with CHARMm */ + if ((*charmm) & DCD_IS_CHARMM) { + float ftmp; + ftmp = *((float *)(hdrbuf+40)); /* is this safe on Alpha? */ + if (*reverseEndian) + swap4_aligned(&ftmp, 1); + + *DELTA = (double)ftmp; + } else { + (*DELTA) = *((double *)(hdrbuf + 40)); + if (*reverseEndian) swap8_unaligned(DELTA, 1); + } + + /* Get the end size of the first block */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading second 84 from dcd file"); + CHECK_FEOF(ret_val, "reading second 84 from dcd file"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != 84) { + return DCD_BADFORMAT; + } + + /* Read in the size of the next block */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading size of title block"); + CHECK_FEOF(ret_val, "reading size of title block"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (((input_integer-4) % 80) == 0) { + /* Read NTITLE, the number of 80 character title strings there are */ + ret_val = READ(fd, &NTITLE, sizeof(int)); + CHECK_FREAD(ret_val, "reading NTITLE"); + CHECK_FEOF(ret_val, "reading NTITLE"); + if (*reverseEndian) swap4_aligned(&NTITLE, 1); + *len_remarks = NTITLE*80; + *remarks = (char*)malloc(*len_remarks); + ret_val = fio_fread(*remarks, *len_remarks, 1, fd); + CHECK_FEOF(ret_val, "reading TITLE"); + + /* Get the ending size for this block */ + ret_val = READ(fd, &input_integer, sizeof(int)); + + CHECK_FREAD(ret_val, "reading size of title block"); + CHECK_FEOF(ret_val, "reading size of title block"); + } else { + return DCD_BADFORMAT; + } + + /* Read in an integer '4' */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading a '4'"); + CHECK_FEOF(ret_val, "reading a '4'"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != 4) { + return DCD_BADFORMAT; + } + + /* Read in the number of atoms */ + ret_val = READ(fd, N, sizeof(int)); + CHECK_FREAD(ret_val, "reading number of atoms"); + CHECK_FEOF(ret_val, "reading number of atoms"); + if (*reverseEndian) swap4_aligned(N, 1); + + /* Read in an integer '4' */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading a '4'"); + CHECK_FEOF(ret_val, "reading a '4'"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != 4) { + return DCD_BADFORMAT; + } + + *FREEINDEXES = NULL; + *fixedcoords = NULL; + if (*NAMNF != 0) { + (*FREEINDEXES) = (int *) calloc(((*N)-(*NAMNF)), sizeof(int)); + if (*FREEINDEXES == NULL) + return DCD_BADMALLOC; + + *fixedcoords = (float *) calloc((*N)*4 - (*NAMNF), sizeof(float)); + if (*fixedcoords == NULL) + return DCD_BADMALLOC; + + /* Read in index array size */ + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading size of index array"); + CHECK_FEOF(ret_val, "reading size of index array"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != ((*N)-(*NAMNF))*4) { + return DCD_BADFORMAT; + } + + ret_val = READ(fd, (*FREEINDEXES), ((*N)-(*NAMNF))*sizeof(int)); + CHECK_FREAD(ret_val, "reading size of index array"); + CHECK_FEOF(ret_val, "reading size of index array"); + + if (*reverseEndian) + swap4_aligned((*FREEINDEXES), ((*N)-(*NAMNF))); + + ret_val = READ(fd, &input_integer, sizeof(int)); + CHECK_FREAD(ret_val, "reading size of index array"); + CHECK_FEOF(ret_val, "reading size of index array"); + if (*reverseEndian) swap4_aligned(&input_integer, 1); + + if (input_integer != ((*N)-(*NAMNF))*4) { + return DCD_BADFORMAT; + } + } + + return DCD_SUCCESS; +} + +static int read_charmm_extrablock(fio_fd fd, int charmm, int reverseEndian, + float *unitcell) { + int i, input_integer; + + if ((charmm & DCD_IS_CHARMM) && (charmm & DCD_HAS_EXTRA_BLOCK)) { + /* Leading integer must be 48 */ + if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) + return DCD_BADREAD; + if (reverseEndian) swap4_aligned(&input_integer, 1); + if (input_integer == 48) { + double tmp[6]; + if (fio_fread(tmp, 48, 1, fd) != 1) return DCD_BADREAD; + if (reverseEndian) + swap8_aligned(tmp, 6); + for (i=0; i<6; i++) unitcell[i] = (float)tmp[i]; + } else { + /* unrecognized block, just skip it */ + if (fio_fseek(fd, input_integer, FIO_SEEK_CUR)) return DCD_BADREAD; + } + if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) return DCD_BADREAD; + } + + return DCD_SUCCESS; +} + +static int read_fixed_atoms(fio_fd fd, int N, int num_free, const int *indexes, + int reverseEndian, const float *fixedcoords, + float *freeatoms, float *pos) { + int i, input_integer; + + /* Read leading integer */ + if (fio_fread(&input_integer, sizeof(int), 1, fd) != 1) return DCD_BADREAD; + if (reverseEndian) swap4_aligned(&input_integer, 1); + if (input_integer != 4*num_free) return DCD_BADFORMAT; + + /* Read free atom coordinates */ + if (fio_fread(freeatoms, 4*num_free, 1, fd) != 1) return DCD_BADREAD; + if (reverseEndian) + swap4_aligned(freeatoms, num_free); + + /* Copy fixed and free atom coordinates into position buffer */ + memcpy(pos, fixedcoords, 4*N); + for (i=0; i 1) { + seekoffset *= numsteps; + } + + rc = fio_fseek(fd, seekoffset, FIO_SEEK_CUR); + if (rc == -1) return DCD_BADEOF; + + return DCD_SUCCESS; +} + +static int jump_to_dcdstep(fio_fd fd, int natoms, int nsets, int nfixed, int charmm, int header_size, int step) { + int rc; + if (step > nsets) { + return DCD_BADEOF; + } + // Calculate file offset + off_t extrablocksize, ndims, firstframesize, framesize; + off_t pos; + extrablocksize = charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0; + ndims = charmm & DCD_HAS_4DIMS ? 4 : 3; + firstframesize = (natoms+2) * ndims * sizeof(float) + extrablocksize; + framesize = (natoms-nfixed+2) * ndims * sizeof(float) + extrablocksize; + // Use zero indexing + if (step == 0) { + pos = header_size; + } + else { + pos = header_size + firstframesize + framesize * (step-1); + } + rc = fio_fseek(fd, pos, FIO_SEEK_SET); + if (rc == -1) return DCD_BADEOF; + return DCD_SUCCESS; +} + +#define NFILE_POS 8L +#define NSTEP_POS 20L + +static int write_dcdstep(fio_fd fd, int curframe, int curstep, int N, + const float *X, const float *Y, const float *Z, + const double *unitcell, int charmm) { + int out_integer; + + if (charmm) { + /* write out optional unit cell */ + if (unitcell != NULL) { + out_integer = 48; /* 48 bytes (6 doubles) */ + fio_write_int32(fd, out_integer); + WRITE(fd, unitcell, out_integer); + fio_write_int32(fd, out_integer); + } + } + + /* write out coordinates */ + out_integer = N*4; /* N*4 bytes per X/Y/Z array (N floats per array) */ + fio_write_int32(fd, out_integer); + WRITE(fd, X, out_integer); + fio_write_int32(fd, out_integer); + fio_write_int32(fd, out_integer); + WRITE(fd, Y, out_integer); + fio_write_int32(fd, out_integer); + fio_write_int32(fd, out_integer); + WRITE(fd, Z, out_integer); + fio_write_int32(fd, out_integer); + + /* update the DCD header information */ + fio_fseek(fd, NFILE_POS, FIO_SEEK_SET); + fio_write_int32(fd, curframe); + fio_fseek(fd, NSTEP_POS, FIO_SEEK_SET); + fio_write_int32(fd, curstep); + fio_fseek(fd, 0, FIO_SEEK_END); + + return DCD_SUCCESS; +} + +static int write_dcdheader(fio_fd fd, const char *remarks, int N, + int ISTART, int NSAVC, double DELTA, int with_unitcell, + int charmm) { + int out_integer; + float out_float; + char title_string[200]; + time_t cur_time; + struct tm *tmbuf; + char time_str[81]; + + out_integer = 84; + WRITE(fd, (char *) & out_integer, sizeof(int)); + strcpy(title_string, "CORD"); + WRITE(fd, title_string, 4); + fio_write_int32(fd, 0); /* Number of frames in file, none written yet */ + fio_write_int32(fd, ISTART); /* Starting timestep */ + fio_write_int32(fd, NSAVC); /* Timesteps between frames written to the file */ + fio_write_int32(fd, 0); /* Number of timesteps in simulation */ + fio_write_int32(fd, 0); /* NAMD writes NSTEP or ISTART - NSAVC here? */ + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + if (charmm) { + out_float = DELTA; + WRITE(fd, (char *) &out_float, sizeof(float)); + if (with_unitcell) { + fio_write_int32(fd, 1); + } else { + fio_write_int32(fd, 0); + } + } else { + WRITE(fd, (char *) &DELTA, sizeof(double)); + } + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + fio_write_int32(fd, 0); + if (charmm) { + fio_write_int32(fd, 24); /* Pretend to be Charmm version 24 */ + } else { + fio_write_int32(fd, 0); + } + fio_write_int32(fd, 84); + fio_write_int32(fd, 164); + fio_write_int32(fd, 2); + + strncpy(title_string, remarks, 80); + title_string[79] = '\0'; + WRITE(fd, title_string, 80); + + cur_time=time(NULL); + tmbuf=localtime(&cur_time); + strftime(time_str, 80, "REMARKS Created %d %B, %Y at %R", tmbuf); + WRITE(fd, time_str, 80); + + fio_write_int32(fd, 164); + fio_write_int32(fd, 4); + fio_write_int32(fd, N); + fio_write_int32(fd, 4); + + return DCD_SUCCESS; +} + +static void close_dcd_read(int *indexes, float *fixedcoords) { + free(indexes); + free(fixedcoords); +} + +#endif + diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx new file mode 100644 index 00000000000..ed403119a3e --- /dev/null +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -0,0 +1,59 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- http://www.MDAnalysis.org +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + +cdef extern from 'include/fastio.h': + ctypedef fio_fd: + static int fio_open(const char *filename, int mode, fio_fd *fd) + static int fio_fclose(fio_fd fd) + +#define FIO_READ 0x01 +#define FIO_WRITE 0x02 + +cdef enum: + FIO_READ = 0x01 + FIO_WRITE = 0x02 + +cdef class DCDFile: + cdef fio_fd fp + cdef readonly fname + cdef int is_open + + def __cinit__(self, fname, mode='r'): + self.fname = fname.encode('utf-8') + self.is_open = False + self.open(self.fname, mode) + + + def __dealloc__(self): + # call a close_dcd_read + self.close() + + def open(self, filename, mode): + # NOTE: to make handling easier lets disallow read/write mode + if mode == 'r': + fio_mode = FIO_READ + elif mode == 'w': + fio_mode = FIO_WRITE + else: + raise IOError("unkown mode '{}', use either r or w".format(mode)) + ok = fio_open(self.fname, fio_mode, self.fp) + if ok != 0: + raise IOError("couldn't open file: {}".format(filename)) + + def close(self): + ok = fio_close(self.fp) + if ok != 0: + raise IOError("couldn't close file: {}".format(self.fname)) From 3fda8a8cc69737809fc304c364669f7842c7b903 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sun, 7 Aug 2016 13:44:57 +0200 Subject: [PATCH 02/26] first working darft of dcd reimplmentation This is a good wirst step --- package/MDAnalysis/lib/formats/libdcd.pyx | 39 ++++++++++++++--------- package/setup.py | 19 ++++++----- 2 files changed, 33 insertions(+), 25 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index ed403119a3e..3b78825c04a 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -14,22 +14,30 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -cdef extern from 'include/fastio.h': - ctypedef fio_fd: - static int fio_open(const char *filename, int mode, fio_fd *fd) - static int fio_fclose(fio_fd fd) - -#define FIO_READ 0x01 -#define FIO_WRITE 0x02 +IF UNAME_SYSNAME == "Windows": + cdef extern from 'windows.h': + ctypedef HANDLE fio_fd + ctypedef LONGLONG fio_size_t; + ctypedef void * fio_caddr_t; +ELSE: + from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END, FILE + ctypedef FILE * fio_fd; + ctypedef size_t fio_size_t; + ctypedef void * fio_caddr_t; + _whence_vals = {"FIO_SEEK_SET": SEEK_SET, + "FIO_SEEK_CUR": SEEK_CUR, + "FIO_SEEK_END": SEEK_END} cdef enum: FIO_READ = 0x01 FIO_WRITE = 0x02 +cdef extern from 'include/fastio.h': + ctypedef struct fio_iovec: + pass + cdef class DCDFile: cdef fio_fd fp - cdef readonly fname - cdef int is_open def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -49,11 +57,12 @@ cdef class DCDFile: fio_mode = FIO_WRITE else: raise IOError("unkown mode '{}', use either r or w".format(mode)) - ok = fio_open(self.fname, fio_mode, self.fp) - if ok != 0: - raise IOError("couldn't open file: {}".format(filename)) + # ok = fio_open(self.fname, fio_mode, self.fp) + # if ok != 0: + # raise IOError("couldn't open file: {}".format(filename)) def close(self): - ok = fio_close(self.fp) - if ok != 0: - raise IOError("couldn't close file: {}".format(self.fname)) + pass + # ok = fio_close(self.fp) + # if ok != 0: + # raise IOError("couldn't close file: {}".format(self.fname)) diff --git a/package/setup.py b/package/setup.py index 834698e1966..17930f2424c 100755 --- a/package/setup.py +++ b/package/setup.py @@ -284,16 +284,16 @@ def extensions(config): # The callable is passed so that it is only evaluated at install time. include_dirs = [get_numpy_include] - dcd = MDAExtension('coordinates._dcdmodule', - ['MDAnalysis/coordinates/src/dcd.c'], - include_dirs=include_dirs + ['MDAnalysis/coordinates/include'], + dcd = MDAExtension('lib.formats.libdcd', + ['MDAnalysis/lib/formats/libdcd' + source_suffix], + include_dirs=include_dirs + ['MDAnalysis/lib/formats/include'], define_macros=define_macros, extra_compile_args=extra_compile_args) - dcd_time = MDAExtension('coordinates.dcdtimeseries', - ['MDAnalysis/coordinates/dcdtimeseries' + source_suffix], - include_dirs=include_dirs + ['MDAnalysis/coordinates/include'], - define_macros=define_macros, - extra_compile_args=extra_compile_args) + # dcd_time = MDAExtension('coordinates.dcdtimeseries', + # ['MDAnalysis/coordinates/dcdtimeseries' + source_suffix], + # include_dirs=include_dirs + ['MDAnalysis/coordinates/include'], + # define_macros=define_macros, + # extra_compile_args=extra_compile_args) distances = MDAExtension('lib.c_distances', ['MDAnalysis/lib/c_distances' + source_suffix], include_dirs=include_dirs + ['MDAnalysis/lib/include'], @@ -332,7 +332,7 @@ def extensions(config): sources=['MDAnalysis/lib/formats/cython_util' + source_suffix], include_dirs=include_dirs) - pre_exts = [dcd, dcd_time, distances, distances_omp, qcprot, + pre_exts = [dcd, distances, distances_omp, qcprot, transformation, libmdaxdr, util] cython_generated = [] if use_cython: @@ -511,4 +511,3 @@ def dynamic_author_list(): except OSError as err: print("Warning: failed to delete cythonized file {0}: {1}. " "Moving on.".format(cythonized, err.strerror)) - From 6ee3416e96403a42e2af93b7c764b93a5b7a2957 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Mon, 8 Aug 2016 07:31:07 +0200 Subject: [PATCH 03/26] open and close dcd files --- package/MDAnalysis/lib/formats/__init__.py | 3 ++- package/MDAnalysis/lib/formats/libdcd.pyx | 28 ++++++++++------------ 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/package/MDAnalysis/lib/formats/__init__.py b/package/MDAnalysis/lib/formats/__init__.py index 719832d8cd5..8d35b4f87ce 100644 --- a/package/MDAnalysis/lib/formats/__init__.py +++ b/package/MDAnalysis/lib/formats/__init__.py @@ -14,5 +14,6 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # from . import libmdaxdr +from . import libdcd -__all__ = ['libmdaxdr'] +__all__ = ['libmdaxdr', 'libdcd'] diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 3b78825c04a..1856aa587df 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -14,16 +14,14 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # +# This is not tested but should serve as an example how to enable windows +# support in the future IF UNAME_SYSNAME == "Windows": cdef extern from 'windows.h': ctypedef HANDLE fio_fd - ctypedef LONGLONG fio_size_t; - ctypedef void * fio_caddr_t; ELSE: from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END, FILE ctypedef FILE * fio_fd; - ctypedef size_t fio_size_t; - ctypedef void * fio_caddr_t; _whence_vals = {"FIO_SEEK_SET": SEEK_SET, "FIO_SEEK_CUR": SEEK_CUR, "FIO_SEEK_END": SEEK_END} @@ -33,11 +31,13 @@ cdef enum: FIO_WRITE = 0x02 cdef extern from 'include/fastio.h': - ctypedef struct fio_iovec: - pass + int fio_open(const char *filename, int mode, fio_fd *fd) + int fio_fclose(fio_fd fd) cdef class DCDFile: cdef fio_fd fp + cdef readonly fname + cdef int is_open def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -46,23 +46,21 @@ cdef class DCDFile: def __dealloc__(self): - # call a close_dcd_read self.close() - def open(self, filename, mode): - # NOTE: to make handling easier lets disallow read/write mode + def open(self, filename, mode='r'): if mode == 'r': fio_mode = FIO_READ elif mode == 'w': fio_mode = FIO_WRITE else: raise IOError("unkown mode '{}', use either r or w".format(mode)) - # ok = fio_open(self.fname, fio_mode, self.fp) - # if ok != 0: - # raise IOError("couldn't open file: {}".format(filename)) + ok = fio_open(self.fname, fio_mode, &self.fp) + if ok != 0: + raise IOError("couldn't open file: {}".format(filename)) def close(self): pass - # ok = fio_close(self.fp) - # if ok != 0: - # raise IOError("couldn't close file: {}".format(self.fname)) + ok = fio_fclose(self.fp) + if ok != 0: + raise IOError("couldn't close file: {}".format(self.fname)) From 9719b79a2e1ef7f69021a4114863e3e45a79042d Mon Sep 17 00:00:00 2001 From: Max Linke Date: Mon, 8 Aug 2016 09:02:50 +0200 Subject: [PATCH 04/26] add context manager support [skip ci] --- package/MDAnalysis/lib/formats/libdcd.pyx | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 1856aa587df..c4bb4cfb0e5 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -44,10 +44,19 @@ cdef class DCDFile: self.is_open = False self.open(self.fname, mode) - def __dealloc__(self): self.close() + def __enter__(self): + """Support context manager""" + return self + + def __exit__(self, exc_type, exc_val, exc_tb): + """Support context manager""" + self.close() + # always propagate exceptions forward + return False + def open(self, filename, mode='r'): if mode == 'r': fio_mode = FIO_READ From 4d1ec71c8350e401d5fed9bfa2564b60b7c77771 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Mon, 8 Aug 2016 09:04:33 +0200 Subject: [PATCH 05/26] remove pass statement --- package/MDAnalysis/lib/formats/libdcd.pyx | 1 - 1 file changed, 1 deletion(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index c4bb4cfb0e5..806c6bb1dc3 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -69,7 +69,6 @@ cdef class DCDFile: raise IOError("couldn't open file: {}".format(filename)) def close(self): - pass ok = fio_fclose(self.fp) if ok != 0: raise IOError("couldn't close file: {}".format(self.fname)) From 6159a2c024c573d50e393f0b24bc2a47c367cacc Mon Sep 17 00:00:00 2001 From: Max Linke Date: Mon, 8 Aug 2016 09:18:51 +0200 Subject: [PATCH 06/26] correctly handle is_open flag cases [skip ci] --- package/MDAnalysis/lib/formats/libdcd.pyx | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 806c6bb1dc3..e85e9c1468b 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -66,9 +66,14 @@ cdef class DCDFile: raise IOError("unkown mode '{}', use either r or w".format(mode)) ok = fio_open(self.fname, fio_mode, &self.fp) if ok != 0: - raise IOError("couldn't open file: {}".format(filename)) + raise IOError("couldn't open file: {}\n" + "ErrorCode: {}".format(self.fname, ok)) + self.is_open = True def close(self): - ok = fio_fclose(self.fp) - if ok != 0: - raise IOError("couldn't close file: {}".format(self.fname)) + if self.is_open: + ok = fio_fclose(self.fp) + self.is_open = False + if ok != 0: + raise IOError("couldn't close file: {}\n" + "ErrorCode: {}".format(self.fname, ok)) From b4642d2ac2a30aeff8c8903628593d9634fa4b16 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Tue, 9 Aug 2016 07:56:39 +0200 Subject: [PATCH 07/26] start dcd read header stub --- package/MDAnalysis/lib/formats/libdcd.pyx | 64 +++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index e85e9c1468b..d9da15afccb 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -30,15 +30,67 @@ cdef enum: FIO_READ = 0x01 FIO_WRITE = 0x02 +cdef enum: +DCD_ERRORS = { + 0: 'No Problem', + -1: 'Normal EOF', + -2: 'DCD file does not exist', + -3: 'Open of DCD file failed', + -4: 'read call on DCD file failed', + -5: 'premature EOF found in DCD file', + -6: 'format of DCD file is wrong', + -7: 'output file already exiss', + -8: 'malloc failed' +} + + cdef extern from 'include/fastio.h': int fio_open(const char *filename, int mode, fio_fd *fd) int fio_fclose(fio_fd fd) +# need to find a typedef for off_t +# fio_size_t fio_ftell(fio_fd fd) + +cdef extern from 'include/readdcd.h': + int read_dcdheader(fio_fd fd, int *natoms, int *nsets, int *istart, + int *nsavc, double *delta, int *nfixed, int **freeind, + float **fixedcoords, int *reverse, int *charmm, + char **remarks, int *len_remarks) + int read_dcdstep(fio_fd fd, int natoms, float *x, float *y, float *z, + float *unitcell, int nfixed, int first, int *freeind, + float *fixedcoords, int reverse, int charmm) + int read_dcdsubset(fio_fd fd, int natoms, int lowerb, int upperb, float *x, + float *y, float *z, float *unitcell, int nfixed, + int first, int *freeind, float *fixedcoords, int reverse, + int charmm) + int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm, int numstep) + void close_dcd_read(int *freeind, float *fixedcoords) + int write_dcdheader(fio_fd fd, const char *remarks, int natoms, + int istart, int nsavc, double delta, int with_unitcell, + int charmm) + int write_dcdstep(fio_fd fd, int curstep, int curframe, + int natoms, const float *x, const float *y, const float *z, + const double *unitcell, int charmm) + int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm, int numsteps) + cdef class DCDFile: cdef fio_fd fp cdef readonly fname cdef int is_open + cdef int natoms + cdef int nsets + cdef int istart + cdef int nsavc + cdef double delta + cdef int nfixed + cdef int **freeind + cdef float **fixedcords + cdef int reverse + cdef int charmm + cdef char **remarks + cdef int len_remarks + def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') self.is_open = False @@ -77,3 +129,15 @@ cdef class DCDFile: if ok != 0: raise IOError("couldn't close file: {}\n" "ErrorCode: {}".format(self.fname, ok)) + + def read_header(self): + if not self.is_open: + raise RuntimeError("No file open") + + ok = read_dcdheader(self.fd, &self.natoms, &self.nsets, + &self.istart, &self.nsavc, self.delta, + self.nfixed, self.freeind, self.fixedcoords, + self.reverse, *self.reverse, + self.remarks, &self.len_remarks) + if ok != 0: + raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) From d30d53ac012b6eeb225e027624c3223e6cd3cce0 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sun, 28 Aug 2016 18:48:54 +0200 Subject: [PATCH 08/26] add include guards for fastio.h --- .../MDAnalysis/lib/formats/include/fastio.h | 66 ++++++++++--------- 1 file changed, 35 insertions(+), 31 deletions(-) diff --git a/package/MDAnalysis/lib/formats/include/fastio.h b/package/MDAnalysis/lib/formats/include/fastio.h index 998b5ebd5fb..4eeda73e21f 100644 --- a/package/MDAnalysis/lib/formats/include/fastio.h +++ b/package/MDAnalysis/lib/formats/include/fastio.h @@ -22,10 +22,13 @@ * ***************************************************************************/ +#ifndef FASTIO_H +#define FASTIO_H + /* Compiling on windows */ #if defined(_MSC_VER) -#if 1 +#if 1 /* use native Windows I/O calls */ #define FASTIO_NATIVEWIN32 1 @@ -51,10 +54,10 @@ typedef struct { static int fio_win32convertfilename(const char *filename, char *newfilename, int maxlen) { int i; int len=strlen(filename); - + if ((len + 1) >= maxlen) return -1; - + for (i=0; i @@ -352,10 +355,10 @@ static int fio_open(const char *filename, int mode, fio_fd *fd) { int nfd; int oflag = 0; - if (mode == FIO_READ) + if (mode == FIO_READ) oflag = O_RDONLY; - if (mode == FIO_WRITE) + if (mode == FIO_WRITE) oflag = O_WRONLY | O_CREAT | O_TRUNC; nfd = open(filename, oflag, 0666); @@ -371,10 +374,10 @@ static int fio_fclose(fio_fd fd) { return close(fd); } -static fio_size_t fio_fread(void *ptr, fio_size_t size, +static fio_size_t fio_fread(void *ptr, fio_size_t size, fio_size_t nitems, fio_fd fd) { int i; - fio_size_t len = 0; + fio_size_t len = 0; int cnt = 0; for (i=0; i= 0) return 0; /* success (emulate behavior of fseek) */ - else + else return -1; /* failure (emulate behavior of fseek) */ } @@ -452,3 +455,4 @@ static int fio_write_str(fio_fd fd, const char *str) { return (fio_fwrite((void *) str, len, 1, fd) != 1); } +#endif // FASTIO_H From a36b04fea64103b7ec0baa6ab68c511ce7c655eb Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sun, 28 Aug 2016 19:08:42 +0200 Subject: [PATCH 09/26] read DCD Header and fixes other issue I set by accident --- package/MDAnalysis/lib/formats/libdcd.pyx | 51 +++++++---------------- 1 file changed, 14 insertions(+), 37 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index d9da15afccb..fec875a9ac6 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -14,23 +14,16 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -# This is not tested but should serve as an example how to enable windows -# support in the future -IF UNAME_SYSNAME == "Windows": - cdef extern from 'windows.h': - ctypedef HANDLE fio_fd -ELSE: - from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END, FILE - ctypedef FILE * fio_fd; - _whence_vals = {"FIO_SEEK_SET": SEEK_SET, - "FIO_SEEK_CUR": SEEK_CUR, - "FIO_SEEK_END": SEEK_END} +from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END +ctypedef int fio_fd; +_whence_vals = {"FIO_SEEK_SET": SEEK_SET, + "FIO_SEEK_CUR": SEEK_CUR, + "FIO_SEEK_END": SEEK_END} cdef enum: FIO_READ = 0x01 FIO_WRITE = 0x02 -cdef enum: DCD_ERRORS = { 0: 'No Problem', -1: 'Normal EOF', @@ -55,40 +48,23 @@ cdef extern from 'include/readdcd.h': int *nsavc, double *delta, int *nfixed, int **freeind, float **fixedcoords, int *reverse, int *charmm, char **remarks, int *len_remarks) - int read_dcdstep(fio_fd fd, int natoms, float *x, float *y, float *z, - float *unitcell, int nfixed, int first, int *freeind, - float *fixedcoords, int reverse, int charmm) - int read_dcdsubset(fio_fd fd, int natoms, int lowerb, int upperb, float *x, - float *y, float *z, float *unitcell, int nfixed, - int first, int *freeind, float *fixedcoords, int reverse, - int charmm) - int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm, int numstep) - void close_dcd_read(int *freeind, float *fixedcoords) - int write_dcdheader(fio_fd fd, const char *remarks, int natoms, - int istart, int nsavc, double delta, int with_unitcell, - int charmm) - int write_dcdstep(fio_fd fd, int curstep, int curframe, - int natoms, const float *x, const float *y, const float *z, - const double *unitcell, int charmm) - int skip_dcdstep(fio_fd fd, int natoms, int nfixed, int charmm, int numsteps) cdef class DCDFile: cdef fio_fd fp cdef readonly fname cdef int is_open - cdef int natoms cdef int nsets cdef int istart cdef int nsavc cdef double delta cdef int nfixed - cdef int **freeind - cdef float **fixedcords + cdef int *freeind + cdef float *fixedcoords cdef int reverse cdef int charmm - cdef char **remarks + cdef char *remarks cdef int len_remarks def __cinit__(self, fname, mode='r'): @@ -134,10 +110,11 @@ cdef class DCDFile: if not self.is_open: raise RuntimeError("No file open") - ok = read_dcdheader(self.fd, &self.natoms, &self.nsets, - &self.istart, &self.nsavc, self.delta, - self.nfixed, self.freeind, self.fixedcoords, - self.reverse, *self.reverse, - self.remarks, &self.len_remarks) + ok = read_dcdheader(self.fp, &self.natoms, &self.nsets, &self.istart, + &self.nsavc, &self.delta, &self.nfixed, &self.freeind, + &self.fixedcoords, &self.reverse, + &self.charmm, &self.remarks, &self.len_remarks) if ok != 0: raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) + + print(self.natoms) From 2c14c726fc0de6d6124ac7834ad99132142b01fd Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 1 Sep 2016 18:37:20 +0200 Subject: [PATCH 10/26] free memory when closing --- package/MDAnalysis/lib/formats/libdcd.pyx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index fec875a9ac6..2e0a508a41d 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -48,6 +48,7 @@ cdef extern from 'include/readdcd.h': int *nsavc, double *delta, int *nfixed, int **freeind, float **fixedcoords, int *reverse, int *charmm, char **remarks, int *len_remarks) + void close_dcd_read(int *freeind, float *fixedcoords); cdef class DCDFile: @@ -100,6 +101,11 @@ cdef class DCDFile: def close(self): if self.is_open: + print("closing", self.freeind == NULL, self.fixedcoords == NULL) + if self.freeind != NULL: + print("freeing memory") + close_dcd_read(self.freeind, self.fixedcoords); + ok = fio_fclose(self.fp) self.is_open = False if ok != 0: From 3fb31e502c39a73b940b5bd9d0987de8f6a1f725 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 1 Sep 2016 19:13:02 +0200 Subject: [PATCH 11/26] remarks handling --- package/MDAnalysis/lib/formats/libdcd.pyx | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 2e0a508a41d..547e09c2f86 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -20,6 +20,9 @@ _whence_vals = {"FIO_SEEK_SET": SEEK_SET, "FIO_SEEK_CUR": SEEK_CUR, "FIO_SEEK_END": SEEK_END} +from libc.stdint cimport uintptr_t +from libc.stdlib cimport free + cdef enum: FIO_READ = 0x01 FIO_WRITE = 0x02 @@ -65,8 +68,6 @@ cdef class DCDFile: cdef float *fixedcoords cdef int reverse cdef int charmm - cdef char *remarks - cdef int len_remarks def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -101,9 +102,7 @@ cdef class DCDFile: def close(self): if self.is_open: - print("closing", self.freeind == NULL, self.fixedcoords == NULL) if self.freeind != NULL: - print("freeing memory") close_dcd_read(self.freeind, self.fixedcoords); ok = fio_fclose(self.fp) @@ -116,11 +115,18 @@ cdef class DCDFile: if not self.is_open: raise RuntimeError("No file open") + cdef char *remarks + cdef int len_remarks + + print("reading", self.freeind, self.fixedcoords == NULL) ok = read_dcdheader(self.fp, &self.natoms, &self.nsets, &self.istart, &self.nsavc, &self.delta, &self.nfixed, &self.freeind, &self.fixedcoords, &self.reverse, - &self.charmm, &self.remarks, &self.len_remarks) + &self.charmm, &remarks, &len_remarks) + print("reading", self.freeind, self.fixedcoords == NULL) if ok != 0: raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) + if remarks != NULL: + free(remarks) print(self.natoms) From 3f717911001e4dcc28122df93509fbc45e4a5c17 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sat, 3 Sep 2016 18:17:19 +0200 Subject: [PATCH 12/26] comment all the things --- package/MDAnalysis/lib/formats/libdcd.pyx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 547e09c2f86..d6e64a1794c 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -102,6 +102,8 @@ cdef class DCDFile: def close(self): if self.is_open: + # In case there are fixed atoms we should free the memory again. + # Both pointers are guaranted to be non NULL if either one is. if self.freeind != NULL: close_dcd_read(self.freeind, self.fixedcoords); From cc1e380951a26bc73847710cefd93929d11ff572 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sat, 1 Oct 2016 21:22:06 +0200 Subject: [PATCH 13/26] convert remarks to python strings --- package/MDAnalysis/lib/formats/libdcd.pyx | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index d6e64a1794c..c93e2403973 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -117,18 +117,20 @@ cdef class DCDFile: if not self.is_open: raise RuntimeError("No file open") - cdef char *remarks - cdef int len_remarks + cdef char* c_remarks + cdef int len_remarks = 0 - print("reading", self.freeind, self.fixedcoords == NULL) ok = read_dcdheader(self.fp, &self.natoms, &self.nsets, &self.istart, &self.nsavc, &self.delta, &self.nfixed, &self.freeind, &self.fixedcoords, &self.reverse, - &self.charmm, &remarks, &len_remarks) - print("reading", self.freeind, self.fixedcoords == NULL) + &self.charmm, &c_remarks, &len_remarks) if ok != 0: raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) - if remarks != NULL: - free(remarks) + if c_remarks != NULL: + py_remarks = c_remarks[:len_remarks] + free(c_remarks) + else: + py_remarks = "" print(self.natoms) + return py_remarks From 5560bcbc73da8d635f344bb202f3bb3cf5aa7295 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sat, 1 Oct 2016 21:27:41 +0200 Subject: [PATCH 14/26] make n_atoms public in dcd --- package/MDAnalysis/lib/formats/libdcd.pyx | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index c93e2403973..43882a5a36e 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -47,7 +47,7 @@ cdef extern from 'include/fastio.h': # fio_size_t fio_ftell(fio_fd fd) cdef extern from 'include/readdcd.h': - int read_dcdheader(fio_fd fd, int *natoms, int *nsets, int *istart, + int read_dcdheader(fio_fd fd, int *n_atoms, int *nsets, int *istart, int *nsavc, double *delta, int *nfixed, int **freeind, float **fixedcoords, int *reverse, int *charmm, char **remarks, int *len_remarks) @@ -58,7 +58,7 @@ cdef class DCDFile: cdef fio_fd fp cdef readonly fname cdef int is_open - cdef int natoms + cdef readonly int n_atoms cdef int nsets cdef int istart cdef int nsavc @@ -71,6 +71,7 @@ cdef class DCDFile: def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') + self.n_atoms = 0 self.is_open = False self.open(self.fname, mode) @@ -120,17 +121,17 @@ cdef class DCDFile: cdef char* c_remarks cdef int len_remarks = 0 - ok = read_dcdheader(self.fp, &self.natoms, &self.nsets, &self.istart, + ok = read_dcdheader(self.fp, &self.n_atoms, &self.nsets, &self.istart, &self.nsavc, &self.delta, &self.nfixed, &self.freeind, &self.fixedcoords, &self.reverse, &self.charmm, &c_remarks, &len_remarks) if ok != 0: raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) + if c_remarks != NULL: py_remarks = c_remarks[:len_remarks] free(c_remarks) else: py_remarks = "" - print(self.natoms) return py_remarks From 1358c67a9ac21ad83ab81238222e62dd625e16e1 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sat, 1 Oct 2016 21:46:23 +0200 Subject: [PATCH 15/26] make delta public --- package/MDAnalysis/lib/formats/libdcd.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 43882a5a36e..5ab1bcf1b7a 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -62,7 +62,7 @@ cdef class DCDFile: cdef int nsets cdef int istart cdef int nsavc - cdef double delta + cdef readonly double delta cdef int nfixed cdef int *freeind cdef float *fixedcoords From 0f5b3b71f4662e507289210362ec514da6d61bb5 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sat, 1 Oct 2016 22:06:24 +0200 Subject: [PATCH 16/26] make more stuff public --- package/MDAnalysis/lib/formats/libdcd.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 5ab1bcf1b7a..142b2578570 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -60,10 +60,10 @@ cdef class DCDFile: cdef int is_open cdef readonly int n_atoms cdef int nsets - cdef int istart - cdef int nsavc + cdef readonly int istart + cdef readonly int nsavc cdef readonly double delta - cdef int nfixed + cdef readonly int nfixed cdef int *freeind cdef float *fixedcoords cdef int reverse From d16560dbda21fd9b341307f5a7b71a8aa6ac2bdf Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sat, 1 Oct 2016 22:11:14 +0200 Subject: [PATCH 17/26] periodic property --- package/MDAnalysis/lib/formats/libdcd.pyx | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 142b2578570..e32e7ecb1cc 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -27,6 +27,11 @@ cdef enum: FIO_READ = 0x01 FIO_WRITE = 0x02 +cdef enum: + DCD_IS_CHARMM = 0x01 + DCD_HAS_4DIMS = 0x02 + DCD_HAS_EXTRA_BLOCK = 0x04 + DCD_ERRORS = { 0: 'No Problem', -1: 'Normal EOF', @@ -134,4 +139,10 @@ cdef class DCDFile: else: py_remarks = "" + return py_remarks + + @property + def periodic(self): + return bool((self.charmm & DCD_IS_CHARMM) and + (self.charmm & DCD_HAS_EXTRA_BLOCK)) From 92bba4d090289e1994630eded6954e291bf7d87a Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sat, 1 Oct 2016 22:40:30 +0200 Subject: [PATCH 18/26] estimate number of frames --- package/MDAnalysis/lib/formats/libdcd.pyx | 34 ++++++++++++++++++++--- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index e32e7ecb1cc..291471d0b22 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -14,12 +14,22 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # +from os import path + from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END -ctypedef int fio_fd; _whence_vals = {"FIO_SEEK_SET": SEEK_SET, "FIO_SEEK_CUR": SEEK_CUR, "FIO_SEEK_END": SEEK_END} +# Tell cython about the off_t type. It doesn't need to match exactly what is +# defined since we don't expose it to python but the cython compiler needs to +# know about it. +cdef extern from 'sys/types.h': + ctypedef int off_t + +ctypedef int fio_fd; +ctypedef off_t fio_size_t + from libc.stdint cimport uintptr_t from libc.stdlib cimport free @@ -44,12 +54,10 @@ DCD_ERRORS = { -8: 'malloc failed' } - cdef extern from 'include/fastio.h': int fio_open(const char *filename, int mode, fio_fd *fd) int fio_fclose(fio_fd fd) -# need to find a typedef for off_t -# fio_size_t fio_ftell(fio_fd fd) + fio_size_t fio_ftell(fio_fd fd) cdef extern from 'include/readdcd.h': int read_dcdheader(fio_fd fd, int *n_atoms, int *nsets, int *istart, @@ -73,6 +81,8 @@ cdef class DCDFile: cdef float *fixedcoords cdef int reverse cdef int charmm + cdef readonly int n_dims + cdef readonly int n_frames def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -139,9 +149,25 @@ cdef class DCDFile: else: py_remarks = "" + self.n_dims = 3 if not self.charmm & DCD_HAS_4DIMS else 4 + self.n_frames = self._estimate_n_frames() return py_remarks + def _estimate_n_frames(self): + extrablocksize = 48 + 8 if self.charmm & DCD_HAS_EXTRA_BLOCK else 0 + firstframesize = self.n_atoms + 2 * self.n_dims * sizeof(float) + extrablocksize + framesize = ((self.n_atoms - self.nfixed + 2) * self.n_dims * sizeof(float) + + extrablocksize) + filesize = path.getsize(self.fname) + # It's safe to use ftell, even though ftell returns a long, because the + # header size is < 4GB. + header_size = fio_ftell(self.fp) + # TODO: check that nframessize is larger then 0, the c-implementation + # used to do that. + nframessize = filesize - header_size - firstframesize + return nframessize / framesize + 1 + @property def periodic(self): return bool((self.charmm & DCD_IS_CHARMM) and From 5ba11d29729ed5c8a1ced8313d5cda8604769b61 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Sat, 1 Oct 2016 23:10:32 +0200 Subject: [PATCH 19/26] work on read_next_frame --- package/MDAnalysis/lib/formats/libdcd.pyx | 45 +++++++++++++++++++++-- 1 file changed, 41 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 291471d0b22..f35682d4479 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -16,6 +16,8 @@ from os import path +cimport numpy as np + from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END _whence_vals = {"FIO_SEEK_SET": SEEK_SET, "FIO_SEEK_CUR": SEEK_CUR, @@ -30,6 +32,9 @@ cdef extern from 'sys/types.h': ctypedef int fio_fd; ctypedef off_t fio_size_t +ctypedef np.float32_t DTYPE_t +DTYPE = np.float32 + from libc.stdint cimport uintptr_t from libc.stdlib cimport free @@ -62,9 +67,13 @@ cdef extern from 'include/fastio.h': cdef extern from 'include/readdcd.h': int read_dcdheader(fio_fd fd, int *n_atoms, int *nsets, int *istart, int *nsavc, double *delta, int *nfixed, int **freeind, - float **fixedcoords, int *reverse, int *charmm, + float **fixedcoords, int *reverse_endian, int *charmm, char **remarks, int *len_remarks) - void close_dcd_read(int *freeind, float *fixedcoords); + void close_dcd_read(int *freeind, float *fixedcoords) + int read_dcdstep(fio_fd fd, int n_atoms, float *X, float *Y, float *Z, + float *unitcell, int num_fixed, + int first, int *indexes, float *fixedcoords, + int reverse_endian, int charmm) cdef class DCDFile: @@ -79,10 +88,11 @@ cdef class DCDFile: cdef readonly int nfixed cdef int *freeind cdef float *fixedcoords - cdef int reverse + cdef int reverse_endian cdef int charmm cdef readonly int n_dims cdef readonly int n_frames + cdef bint b_read_header def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -115,6 +125,7 @@ cdef class DCDFile: raise IOError("couldn't open file: {}\n" "ErrorCode: {}".format(self.fname, ok)) self.is_open = True + self.read_header = False def close(self): if self.is_open: @@ -138,7 +149,7 @@ cdef class DCDFile: ok = read_dcdheader(self.fp, &self.n_atoms, &self.nsets, &self.istart, &self.nsavc, &self.delta, &self.nfixed, &self.freeind, - &self.fixedcoords, &self.reverse, + &self.fixedcoords, &self.reverse_endian, &self.charmm, &c_remarks, &len_remarks) if ok != 0: raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) @@ -152,6 +163,7 @@ cdef class DCDFile: self.n_dims = 3 if not self.charmm & DCD_HAS_4DIMS else 4 self.n_frames = self._estimate_n_frames() + self.read_header = True return py_remarks def _estimate_n_frames(self): @@ -172,3 +184,28 @@ cdef class DCDFile: def periodic(self): return bool((self.charmm & DCD_IS_CHARMM) and (self.charmm & DCD_HAS_EXTRA_BLOCK)) + + + def _read_next_frame(self): + if not self.is_open: + raise RuntimeError("No file open") + if not self.read_header: + raise IOError("didn't read DCD header before reading frame") + + cdef np.ndarray xyz = np.empty((self.n_atoms, 3), dtype=DTYPE, + order='F') + cdef np.ndarray dimensions = np.empty(6, dtype=DTYPE) + + cdef DTYPE_t[::1] x = xyz[:, 0] + cdef DTYPE_t[::1] y = xyz[:, 1] + cdef DTYPE_t[::1] z = xyz[:, 2] + + ok = read_dcdstep(self.fp, self.n_atoms, &x[0], + &y[0], &z[0], + dimensions.data, self.nfixed, self.first, + self.freeind, self.fixedcoords, + self.reverse_endian, self.charmm) + if ok != 0: + raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) + + return xyz From b681465e3a8833a4b8c4351b5eecd224b2076430 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Tue, 4 Oct 2016 07:59:07 +0200 Subject: [PATCH 20/26] improve reading API --- package/MDAnalysis/lib/formats/libdcd.pyx | 26 +++++++++++++++++------ 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index f35682d4479..7089feef540 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -15,9 +15,12 @@ # from os import path +import numpy as np +from collections import namedtuple cimport numpy as np + from libc.stdio cimport SEEK_SET, SEEK_CUR, SEEK_END _whence_vals = {"FIO_SEEK_SET": SEEK_SET, "FIO_SEEK_CUR": SEEK_CUR, @@ -75,6 +78,7 @@ cdef extern from 'include/readdcd.h': int first, int *indexes, float *fixedcoords, int reverse_endian, int charmm) +DCDFrame = namedtuple('DCDFrame', 'x unitcell') cdef class DCDFile: cdef fio_fd fp @@ -93,6 +97,8 @@ cdef class DCDFile: cdef readonly int n_dims cdef readonly int n_frames cdef bint b_read_header + cdef int current_frame + cdef readonly remarks def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -125,7 +131,8 @@ cdef class DCDFile: raise IOError("couldn't open file: {}\n" "ErrorCode: {}".format(self.fname, ok)) self.is_open = True - self.read_header = False + self.current_frame = 0 + self.remarks = self._read_header() def close(self): if self.is_open: @@ -140,7 +147,8 @@ cdef class DCDFile: raise IOError("couldn't close file: {}\n" "ErrorCode: {}".format(self.fname, ok)) - def read_header(self): + + def _read_header(self): if not self.is_open: raise RuntimeError("No file open") @@ -162,8 +170,8 @@ cdef class DCDFile: self.n_dims = 3 if not self.charmm & DCD_HAS_4DIMS else 4 self.n_frames = self._estimate_n_frames() + self.b_read_header = True - self.read_header = True return py_remarks def _estimate_n_frames(self): @@ -186,10 +194,10 @@ cdef class DCDFile: (self.charmm & DCD_HAS_EXTRA_BLOCK)) - def _read_next_frame(self): + def read(self): if not self.is_open: raise RuntimeError("No file open") - if not self.read_header: + if not self.b_read_header: raise IOError("didn't read DCD header before reading frame") cdef np.ndarray xyz = np.empty((self.n_atoms, 3), dtype=DTYPE, @@ -200,12 +208,16 @@ cdef class DCDFile: cdef DTYPE_t[::1] y = xyz[:, 1] cdef DTYPE_t[::1] z = xyz[:, 2] + first_frame = self.current_frame == 0 + ok = read_dcdstep(self.fp, self.n_atoms, &x[0], &y[0], &z[0], - dimensions.data, self.nfixed, self.first, + dimensions.data, self.nfixed, first_frame, self.freeind, self.fixedcoords, self.reverse_endian, self.charmm) if ok != 0: raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) - return xyz + self.current_frame += 1 + + return DCDFrame(xyz, dimensions) From 794c533878cc2371a4a7cd37f48726a45f5ff232 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Tue, 4 Oct 2016 08:07:59 +0200 Subject: [PATCH 21/26] finish reading and iteration protocol --- package/MDAnalysis/lib/formats/libdcd.pyx | 46 ++++++++++++++++++++--- 1 file changed, 40 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 7089feef540..39a75b3bc11 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -94,11 +94,13 @@ cdef class DCDFile: cdef float *fixedcoords cdef int reverse_endian cdef int charmm + cdef str mode cdef readonly int n_dims cdef readonly int n_frames cdef bint b_read_header cdef int current_frame cdef readonly remarks + cdef int reached_eof def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -119,13 +121,36 @@ cdef class DCDFile: # always propagate exceptions forward return False + def __iter__(self): + self.close() + self.open(self.fname, self.mode) + return self + + def __next__(self): + if self.reached_eof: + raise StopIteration + return self.read() + + def __len__(self): + if not self.is_open: + raise RuntimeError('No file currently opened') + return self.n_frames + + def tell(self): + return self.current_frame + def open(self, filename, mode='r'): + if self.is_open: + self.close() + if mode == 'r': fio_mode = FIO_READ elif mode == 'w': fio_mode = FIO_WRITE else: raise IOError("unkown mode '{}', use either r or w".format(mode)) + self.mode = mode + ok = fio_open(self.fname, fio_mode, &self.fp) if ok != 0: raise IOError("couldn't open file: {}\n" @@ -133,6 +158,7 @@ cdef class DCDFile: self.is_open = True self.current_frame = 0 self.remarks = self._read_header() + self.reached_eof = False def close(self): if self.is_open: @@ -195,14 +221,17 @@ cdef class DCDFile: def read(self): + if self.reached_eof: + raise RuntimeError('Reached last frame in TRR, seek to 0') if not self.is_open: raise RuntimeError("No file open") - if not self.b_read_header: - raise IOError("didn't read DCD header before reading frame") + if self.mode != 'r': + raise RuntimeError('File opened in mode: {}. Reading only allow ' + 'in mode "r"'.format('self.mode')) cdef np.ndarray xyz = np.empty((self.n_atoms, 3), dtype=DTYPE, order='F') - cdef np.ndarray dimensions = np.empty(6, dtype=DTYPE) + cdef np.ndarray unitcell = np.empty(6, dtype=DTYPE) cdef DTYPE_t[::1] x = xyz[:, 0] cdef DTYPE_t[::1] y = xyz[:, 1] @@ -212,12 +241,17 @@ cdef class DCDFile: ok = read_dcdstep(self.fp, self.n_atoms, &x[0], &y[0], &z[0], - dimensions.data, self.nfixed, first_frame, + unitcell.data, self.nfixed, first_frame, self.freeind, self.fixedcoords, self.reverse_endian, self.charmm) - if ok != 0: + if ok != 0 and ok != -4: raise IOError("Reading DCD header failed: {}".format(DCD_ERRORS[ok])) + # we couldn't read any more frames. + if ok == -4: + self.reached_eof = True + raise StopIteration + self.current_frame += 1 - return DCDFrame(xyz, dimensions) + return DCDFrame(xyz, unitcell) From 629b8541146ae04ae390c8f7ec8a754b3ad93f4c Mon Sep 17 00:00:00 2001 From: Max Linke Date: Tue, 4 Oct 2016 08:21:53 +0200 Subject: [PATCH 22/26] add seek support --- package/MDAnalysis/lib/formats/libdcd.pyx | 53 ++++++++++++++++++++--- 1 file changed, 47 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/lib/formats/libdcd.pyx b/package/MDAnalysis/lib/formats/libdcd.pyx index 39a75b3bc11..d07dfb1c80d 100644 --- a/package/MDAnalysis/lib/formats/libdcd.pyx +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -66,6 +66,7 @@ cdef extern from 'include/fastio.h': int fio_open(const char *filename, int mode, fio_fd *fd) int fio_fclose(fio_fd fd) fio_size_t fio_ftell(fio_fd fd) + fio_size_t fio_fseek(fio_fd fd, fio_size_t offset, int whence) cdef extern from 'include/readdcd.h': int read_dcdheader(fio_fd fd, int *n_atoms, int *nsets, int *istart, @@ -101,6 +102,9 @@ cdef class DCDFile: cdef int current_frame cdef readonly remarks cdef int reached_eof + cdef int firstframesize + cdef int framesize + cdef int header_size def __cinit__(self, fname, mode='r'): self.fname = fname.encode('utf-8') @@ -198,21 +202,25 @@ cdef class DCDFile: self.n_frames = self._estimate_n_frames() self.b_read_header = True + # make sure fixed atoms have been read + self.read() + self.seek(0) + return py_remarks def _estimate_n_frames(self): extrablocksize = 48 + 8 if self.charmm & DCD_HAS_EXTRA_BLOCK else 0 - firstframesize = self.n_atoms + 2 * self.n_dims * sizeof(float) + extrablocksize - framesize = ((self.n_atoms - self.nfixed + 2) * self.n_dims * sizeof(float) + - extrablocksize) + self.firstframesize = self.n_atoms + 2 * self.n_dims * sizeof(float) + extrablocksize + self.framesize = ((self.n_atoms - self.nfixed + 2) * self.n_dims * sizeof(float) + + extrablocksize) filesize = path.getsize(self.fname) # It's safe to use ftell, even though ftell returns a long, because the # header size is < 4GB. - header_size = fio_ftell(self.fp) + self.header_size = fio_ftell(self.fp) # TODO: check that nframessize is larger then 0, the c-implementation # used to do that. - nframessize = filesize - header_size - firstframesize - return nframessize / framesize + 1 + nframessize = filesize - self.header_size - self.firstframesize + return nframessize / self.framesize + 1 @property def periodic(self): @@ -255,3 +263,36 @@ cdef class DCDFile: self.current_frame += 1 return DCDFrame(xyz, unitcell) + + def seek(self, frame): + """Seek to Frame. + + Please note that this function will generate internal file offsets if + they haven't been set before. For large file this means the first seek + can be very slow. Later seeks will be very fast. + + Parameters + ---------- + frame : int + seek the file to given frame + + Raises + ------ + IOError + If you seek for more frames than are available or if the + seek fails (the low-level system error is reported). + """ + if frame >= self.n_frames: + raise IOError('Trying to seek over max number of frames') + self.reached_eof = False + + cdef fio_size_t offset + if frame == 0: + offset = self.header_size + else: + offset = self.header_size + self.firstframesize + self.framesize * (frame - 1) + + ok = fio_fseek(self.fp, offset, _whence_vals["FIO_SEEK_SET"]) + if ok != 0: + raise IOError("DCD seek failed with system errno={}".format(ok)) + self.current_frame = frame From 65926f6ad4b8f27de5668fed22e9a3b9675c0acb Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Fri, 21 Oct 2016 16:45:15 +0100 Subject: [PATCH 23/26] TST: Added first draft of testing module for new Cython-based DCD file handling interface. --- .../MDAnalysisTests/formats/test_libdcd.py | 53 +++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 testsuite/MDAnalysisTests/formats/test_libdcd.py diff --git a/testsuite/MDAnalysisTests/formats/test_libdcd.py b/testsuite/MDAnalysisTests/formats/test_libdcd.py new file mode 100644 index 00000000000..0b7357de306 --- /dev/null +++ b/testsuite/MDAnalysisTests/formats/test_libdcd.py @@ -0,0 +1,53 @@ +from __future__ import print_function + +from nose.tools import raises +from numpy.testing import assert_equal, assert_array_equal +from numpy.testing import assert_array_almost_equal +from numpy.testing import assert_almost_equal + +from MDAnalysis.lib.formats.libdcd import DCDFile +from MDAnalysisTests.datafiles import DCD + +from unittest import TestCase +from MDAnalysisTests.tempdir import run_in_tempdir +import numpy as np + +class DCDReadFrameTest(TestCase): + + def setUp(self): + self.dcdfile = DCDFile(DCD) + + def tearDown(self): + del self.dcdfile + + def test_read_coords(self): + # confirm shape of coordinate data against result from previous + # MDAnalysis implementation of DCD file handling + dcd_frame = self.dcdfile.read() + xyz = dcd_frame[0] + assert_equal(xyz.shape, (3341, 3)) + + def test_read_unit_cell(self): + # confirm unit cell read against result from previous + # MDAnalysis implementation of DCD file handling + dcd_frame = self.dcdfile.read() + unitcell = dcd_frame[1] + expected = np.array([ 0., 0., 0., 90., 90., 90.], + dtype=np.float32) + assert_equal(unitcell, expected) + + def test_seek_over_max(self): + # should raise IOError if beyond 98th frame + with self.assertRaises(IOError): + self.dcdfile.seek(102) + + def test_seek_normal(self): + # frame seek within range is tested + new_frame = 91 + self.dcdfile.seek(new_frame) + assert_equal(self.dcdfile.current_frame, new_frame) + + def test_seek_negative(self): + # frame seek with negative number + with self.assertRaises(IOError): + self.dcdfile.seek(-78) From 5d669c581c382cd4907a77a7cad3130440c7a59a Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Sun, 23 Oct 2016 18:13:14 +0100 Subject: [PATCH 24/26] TST: All libdcd unit tests now pass, except one related to unit cell dimensions. --- testsuite/MDAnalysisTests/formats/test_libdcd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/formats/test_libdcd.py b/testsuite/MDAnalysisTests/formats/test_libdcd.py index 0b7357de306..4e407b888a6 100644 --- a/testsuite/MDAnalysisTests/formats/test_libdcd.py +++ b/testsuite/MDAnalysisTests/formats/test_libdcd.py @@ -45,7 +45,7 @@ def test_seek_normal(self): # frame seek within range is tested new_frame = 91 self.dcdfile.seek(new_frame) - assert_equal(self.dcdfile.current_frame, new_frame) + assert_equal(self.dcdfile.tell(), new_frame) def test_seek_negative(self): # frame seek with negative number From 111c44a1d65dcaa6081c677a91eceffd1b8ad6f6 Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Sun, 23 Oct 2016 18:53:39 +0100 Subject: [PATCH 25/26] TST: Added several more unit tests to test_libdcd.py. --- .../MDAnalysisTests/formats/test_libdcd.py | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/testsuite/MDAnalysisTests/formats/test_libdcd.py b/testsuite/MDAnalysisTests/formats/test_libdcd.py index 4e407b888a6..19550237677 100644 --- a/testsuite/MDAnalysisTests/formats/test_libdcd.py +++ b/testsuite/MDAnalysisTests/formats/test_libdcd.py @@ -51,3 +51,51 @@ def test_seek_negative(self): # frame seek with negative number with self.assertRaises(IOError): self.dcdfile.seek(-78) + + def test_iteration(self): + self.dcdfile.__next__() + self.dcdfile.__next__() + self.dcdfile.__next__() + expected_frame = 3 + assert_equal(self.dcdfile.tell(), expected_frame) + + def test_zero_based_frames(self): + expected_frame = 0 + assert_equal(self.dcdfile.tell(), expected_frame) + + def test_length_traj(self): + expected = 98 + assert_equal(len(self.dcdfile), expected) + + def test_context_manager(self): + frame = 22 + with self.dcdfile as f: + f.seek(frame) + assert_equal(f.tell(), frame) + + @raises(ValueError) + def test_open_wrong_mode(self): + DCDFile('foo', 'e') + + @raises(IOError) + def test_raise_not_existing(self): + DCDFile('foo') + + def test_n_atoms(self): + assert_equal(self.dcdfile.n_atoms, 3341) + + @raises(RuntimeError) + @run_in_tempdir() + def test_read_write_mode_file(self): + with DCDFile('foo', 'w') as f: + f.read() + + @raises(RuntimeError) + def test_read_closed(self): + self.dcdfile.close() + self.dcdfile.read() + + def test_iteration_2(self): + with self.dcdfile as f: + for frame in f: + pass From 78a4257fdedda08a7b73a78ebd39b4c1261388fa Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Mon, 24 Oct 2016 16:26:44 +0100 Subject: [PATCH 26/26] TST: updated unit tests to reflect appropriate Error raising for test_read_write_mode_file and test_open_wrong_mode. --- testsuite/MDAnalysisTests/formats/test_libdcd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testsuite/MDAnalysisTests/formats/test_libdcd.py b/testsuite/MDAnalysisTests/formats/test_libdcd.py index 19550237677..1e7cbf11c11 100644 --- a/testsuite/MDAnalysisTests/formats/test_libdcd.py +++ b/testsuite/MDAnalysisTests/formats/test_libdcd.py @@ -73,7 +73,7 @@ def test_context_manager(self): f.seek(frame) assert_equal(f.tell(), frame) - @raises(ValueError) + @raises(IOError) def test_open_wrong_mode(self): DCDFile('foo', 'e') @@ -84,7 +84,7 @@ def test_raise_not_existing(self): def test_n_atoms(self): assert_equal(self.dcdfile.n_atoms, 3341) - @raises(RuntimeError) + @raises(IOError) @run_in_tempdir() def test_read_write_mode_file(self): with DCDFile('foo', 'w') as f: