WIP: Dcd cython

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']

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 <math.h>
+/* Python.h for 'typedef int Py_intptr_t;'  (fixes Issue 19) */
+#include <Python.h>
+
+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<numdata;i++) {
+    code = datacode[i];
+    switch (code) {
+    case 'm':
+      x1 = y1 = z1 = 0.0;
+      aux2 = 0;
+      for (j=0;j<atomcounts[i];j++) {
+        index1 = atomlist[atomno]-lowerb;
+        aux1 = aux[atomno++];
+        aux2 += aux1;
+        x1 += tempX[index1]*aux1;
+        y1 += tempY[index1]*aux1;
+        z1 += tempZ[index1]*aux1;
+      }
+      *(double*)(data + dataset++*stride0 + frame*stride1) = x1/aux2;
+      *(double*)(data + dataset++*stride0 + frame*stride1) = y1/aux2;
+      *(double*)(data + dataset++*stride0 + frame*stride1) = z1/aux2;
+      break;
+    case 'x':
+      index1 = atomlist[atomno++]-lowerb;
+      *(double*)(data+dataset++*stride0+frame*stride1) = tempX[index1];
+      break;
+    case 'y':
+      index1 = atomlist[atomno++]-lowerb;
+      *(double*)(data+dataset++*stride0+frame*stride1) = tempY[index1];
+      break;
+    case 'z':
+      index1 = atomlist[atomno++]-lowerb;
+      *(double*)(data+dataset++*stride0+frame*stride1) = tempZ[index1];
+      break;
+    case 'v':
+      index1 = atomlist[atomno++]-lowerb;
+      *(double*)(data + dataset++*stride0 + frame*stride1) = tempX[index1];
+      *(double*)(data + dataset++*stride0 + frame*stride1) = tempY[index1];
+      *(double*)(data + dataset++*stride0 + frame*stride1) = tempZ[index1];
+      break;
+    case 'a':
+      index1 = atomlist[atomno++]-lowerb;
+      index2 = atomlist[atomno++]-lowerb;
+      index3 = atomlist[atomno++]-lowerb;
+      x1 = tempX[index1]-tempX[index2];
+      y1 = tempY[index1]-tempY[index2];
+      z1 = tempZ[index1]-tempZ[index2];
+      x2 = tempX[index3]-tempX[index2];
+      y2 = tempY[index3]-tempY[index2];
+      z2 = tempZ[index3]-tempZ[index2];
+      aux1 = sqrt(x1*x1+y1*y1+z1*z1); 
+      aux2 = sqrt(x2*x2+y2*y2+z2*z2);
+      *(double*)(data + dataset++*stride0 + frame*stride1) = acos((x1*x2+y1*y2+z1*z2)/(aux1*aux2));
+      break;
+    case 'd':
+      index1 = atomlist[atomno++]-lowerb;
+      index2 = atomlist[atomno++]-lowerb;
+      *(double*)(data + dataset++*stride0+frame*stride1) = tempX[index2]-tempX[index1];
+      *(double*)(data + dataset++*stride0+frame*stride1) = tempY[index2]-tempY[index1];
+      *(double*)(data + dataset++*stride0+frame*stride1) = tempZ[index2]-tempZ[index1];
+      break;
+    case 'r':
+      index1 = atomlist[atomno++]-lowerb;
+      index2 = atomlist[atomno++]-lowerb;
+      x1 = tempX[index2]-tempX[index1];
+      y1 = tempY[index2]-tempY[index1];
+      z1 = tempZ[index2]-tempZ[index1];
+      *(double*)(data + dataset++*stride0+frame*stride1) = sqrt(x1*x1+y1*y1+z1*z1);
+      break;
+    case 'h':
+      index1 = atomlist[atomno++]-lowerb;
+      index2 = atomlist[atomno++]-lowerb;
+      index3 = atomlist[atomno++]-lowerb;
+      index4 = atomlist[atomno++]-lowerb;
+      x1 = tempX[index2]-tempX[index1]; 
+      y1 = tempY[index2]-tempY[index1]; 
+      z1 = tempZ[index2]-tempZ[index1];
+      x2 = tempX[index3]-tempX[index2]; 
+      y2 = tempY[index3]-tempY[index2]; 
+      z2 = tempZ[index3]-tempZ[index2];
+      x3 = tempX[index3]-tempX[index4]; 
+      y3 = tempY[index3]-tempY[index4]; 
+      z3 = tempZ[index3]-tempZ[index4];
+      double nx1, ny1, nz1, nx2, ny2, nz2;
+      // v1 x v2
+      nx1 = (y1*z2) - (y2*z1); ny1 = (z1*x2)-(z2*x1); nz1 = (x1*y2)-(x2*y1);
+      // v3 x v2
+      nx2 = (y3*z2) - (y2*z3); ny2 = (z3*x2)-(z2*x3); nz2 = (x3*y2)-(x2*y3);
+      double a, b;
+      a = sqrt(nx1*nx1+ny1*ny1+nz1*nz1); b = sqrt(nx2*nx2+ny2*ny2+nz2*nz2);
+      // normalized the cross products
+      nx1 /= a; ny1 /= a; nz1 /= a; nx2 /= b; ny2 /= b; nz2 /= b;
+      // find the angle
+      aux1 = acos(nx1*nx2+ny1*ny2+nz1*nz2);
+      // and the sign of the angle
+      aux2 = (nx2*x1+ny2*y1+nz2*z1);
+      if ((aux2 < 0 && aux1 > 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

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<ndata; i++) {
+    tmp = dataptr[0];
+    dataptr[0] = dataptr[3];
+    dataptr[3] = tmp;
+    tmp = dataptr[1];
+    dataptr[1] = dataptr[2];
+    dataptr[2] = tmp;
+    dataptr += 4;
+  }
+}
+
+
+/* works on unaligned 8-byte quantities */
+static void swap8_unaligned(void *v, long ndata) {
+  char *data = (char *) v;
+  long i;
+  char byteArray[8];
+  char *bytePointer;
+
+  for (i=0; i<ndata; i++) {
+    bytePointer = data + (i<<3);
+    byteArray[0]  =  *bytePointer;
+    byteArray[1]  =  *(bytePointer+1);
+    byteArray[2]  =  *(bytePointer+2);
+    byteArray[3]  =  *(bytePointer+3);
+    byteArray[4]  =  *(bytePointer+4);
+    byteArray[5]  =  *(bytePointer+5);
+    byteArray[6]  =  *(bytePointer+6);
+    byteArray[7]  =  *(bytePointer+7);
+
+    *bytePointer     = byteArray[7];
+    *(bytePointer+1) = byteArray[6];
+    *(bytePointer+2) = byteArray[5];
+    *(bytePointer+3) = byteArray[4];
+    *(bytePointer+4) = byteArray[3];
+    *(bytePointer+5) = byteArray[2];
+    *(bytePointer+6) = byteArray[1];
+    *(bytePointer+7) = byteArray[0];
+  }
+}
+
+
+/* Only works with aligned 2-byte quantities, will cause a bus error */
+/* on some platforms if used on unaligned data.                      */
+static void swap2_aligned(void *v, long ndata) {
+  short *data = (short *) v;
+  long i;
+  short *N; 
+
+  for (i=0; i<ndata; i++) {
+    N = data + i;
+    *N=(((*N>>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<ndata; i++) {
+    N = data + i;
+    *N=(((*N>>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<ndata; i++) {
+    N = data + (i<<1);
+    t0 = N[0];
+    t0=(((t0>>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
+