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/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..4eeda73e21f --- /dev/null +++ b/package/MDAnalysis/lib/formats/include/fastio.h @@ -0,0 +1,458 @@ +/*************************************************************************** + *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. + * + ***************************************************************************/ + +#ifndef FASTIO_H +#define FASTIO_H + +/* 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); +} + +#endif // FASTIO_H 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..d07dfb1c80d --- /dev/null +++ b/package/MDAnalysis/lib/formats/libdcd.pyx @@ -0,0 +1,298 @@ +# -*- 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 +# + +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, + "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 + +ctypedef np.float32_t DTYPE_t +DTYPE = np.float32 + +from libc.stdint cimport uintptr_t +from libc.stdlib cimport free + +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', + -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) + 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, + int *nsavc, double *delta, int *nfixed, int **freeind, + float **fixedcoords, int *reverse_endian, int *charmm, + char **remarks, int *len_remarks) + 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) + +DCDFrame = namedtuple('DCDFrame', 'x unitcell') + +cdef class DCDFile: + cdef fio_fd fp + cdef readonly fname + cdef int is_open + cdef readonly int n_atoms + cdef int nsets + cdef readonly int istart + cdef readonly int nsavc + cdef readonly double delta + cdef readonly int nfixed + cdef int *freeind + 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 + cdef int firstframesize + cdef int framesize + cdef int header_size + + 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) + + 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 __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" + "ErrorCode: {}".format(self.fname, ok)) + self.is_open = True + self.current_frame = 0 + self.remarks = self._read_header() + self.reached_eof = False + + 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); + + ok = fio_fclose(self.fp) + self.is_open = False + 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") + + cdef char* c_remarks + cdef int len_remarks = 0 + + ok = read_dcdheader(self.fp, &self.n_atoms, &self.nsets, &self.istart, + &self.nsavc, &self.delta, &self.nfixed, &self.freeind, + &self.fixedcoords, &self.reverse_endian, + &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 = "" + + 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 + + # 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 + 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. + self.header_size = fio_ftell(self.fp) + # TODO: check that nframessize is larger then 0, the c-implementation + # used to do that. + nframessize = filesize - self.header_size - self.firstframesize + return nframessize / self.framesize + 1 + + @property + def periodic(self): + return bool((self.charmm & DCD_IS_CHARMM) and + (self.charmm & DCD_HAS_EXTRA_BLOCK)) + + + 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 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 unitcell = 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] + + first_frame = self.current_frame == 0 + + ok = read_dcdstep(self.fp, self.n_atoms, &x[0], + &y[0], &z[0], + unitcell.data, self.nfixed, first_frame, + self.freeind, self.fixedcoords, + self.reverse_endian, self.charmm) + 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, 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 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)) - diff --git a/testsuite/MDAnalysisTests/formats/test_libdcd.py b/testsuite/MDAnalysisTests/formats/test_libdcd.py new file mode 100644 index 00000000000..1e7cbf11c11 --- /dev/null +++ b/testsuite/MDAnalysisTests/formats/test_libdcd.py @@ -0,0 +1,101 @@ +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.tell(), new_frame) + + 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(IOError) + 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(IOError) + @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