Skip to content

Commit

Permalink
DCD.py changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@jdetle
jdetle committed Jul 27, 2016
1 parent 2a5fc32 commit b7b58ad
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 36 deletions.
11 changes: 8 additions & 3 deletions package/MDAnalysis/coordinates/DCD.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -498,7 +498,8 @@ def _read_frame(self, frame):
ts.frame = frame
return ts

def timeseries(self, asel, start=None, stop=None, step=None, format='afc'):
def timeseries(self, asel, start=None, stop=-1, step=None, skip=None,
format='afc'):
"""Return a subset of coordinate data for an AtomGroup
:Arguments:
Expand All @@ -514,7 +515,11 @@ def timeseries(self, asel, start=None, stop=None, step=None, format='afc'):
where the shape is (frame, number of atoms,
coordinates)
"""
start, stop, skip = self.check_slice_indices(start, stop, step)
if skip is not None:
step = skip
start, stop, step = self.check_slice_indices(start, stop, step)
# for old cython _read_timeseries code
skip = step
if len(asel) == 0:
raise NoDataError("Timeseries requires at least one atom to analyze")
if len(format) != 3 and format not in ['afc', 'acf', 'caf', 'cfa', 'fac', 'fca']:
Expand All @@ -523,7 +528,7 @@ def timeseries(self, asel, start=None, stop=None, step=None, format='afc'):
# Check if the atom numbers can be grouped for efficiency, then we can read partial buffers
# from trajectory file instead of an entire timestep
# XXX needs to be implemented
return self._read_timeseries(atom_numbers, start, stop, step, format)
return self._read_timeseries(atom_numbers, start, stop, skip, format)

def correl(self, timeseries, start=None, stop=None, step=None):
"""Populate a TimeseriesCollection object with timeseries computed from the trajectory
Expand Down
65 changes: 32 additions & 33 deletions package/MDAnalysis/coordinates/src/dcd.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -535,7 +535,7 @@ __read_timeseries(PyObject *self, PyObject *args)
for (i=0;i<n_frames;i++) {
if (skip > 1) {
// Figure out how many steps to skip
numskip = skip-1;
numskip = skip;
if(remaining_frames < numskip){
numskip = 1;
}
Expand All @@ -554,50 +554,49 @@ __read_timeseries(PyObject *self, PyObject *args)

dcd->first = 0;
dcd->setsread += numskip;
remaining_frames = stop - dcd->setsread++;
remaining_frames = stop - dcd->setsread;
if (rc < 0) {
// return an exception
PyErr_SetString(PyExc_IOError, "Error reading frame from DCD file");
goto error;
}
}

// Copy into Numeric array only those atoms we are interested in
for (j=0;j<n_atoms;j++) {
index = atomlist[j]-lowerb;
}
/*
* coord[a][b][c] = *(float*)(coord->data + a*coord->strides[0] + b*coord->strides[1] + c*coord->strides[2])
*/
if (strncasecmp(format, "afc", 3) == 0) {
*(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 0*coord->strides[2]) = tempX[index];
*(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 1*coord->strides[2]) = tempY[index];
*(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 2*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "acf", 3) == 0) {
*(double*)(coord->data + j*coord->strides[0] + 0*coord->strides[1] + i*coord->strides[2]) = tempX[index];
*(double*)(coord->data + j*coord->strides[0] + 1*coord->strides[1] + i*coord->strides[2]) = tempY[index];
*(double*)(coord->data + j*coord->strides[0] + 2*coord->strides[1] + i*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "fac", 3) == 0) {
*(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 0*coord->strides[2]) = tempX[index];
*(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 1*coord->strides[2]) = tempY[index];
*(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 2*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "fca", 3) == 0) {
*(double*)(coord->data + i*coord->strides[0] + 0*coord->strides[1] + j*coord->strides[2]) = tempX[index];
*(double*)(coord->data + i*coord->strides[0] + 1*coord->strides[1] + j*coord->strides[2]) = tempY[index];
*(double*)(coord->data + i*coord->strides[0] + 2*coord->strides[1] + j*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "caf", 3) == 0) {
*(double*)(coord->data + 0*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempX[index];
*(double*)(coord->data + 1*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempY[index];
*(double*)(coord->data + 2*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "cfa", 3) == 0) {
*(double*)(coord->data + 0*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempX[index];
*(double*)(coord->data + 1*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempY[index];
*(double*)(coord->data + 2*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempZ[index];
}

/*
* coord[a][b][c] = *(float*)(coord->data + a*coord->strides[0] + b*coord->strides[1] + c*coord->strides[2])
*/
if (strncasecmp(format, "afc", 3) == 0) {
*(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 0*coord->strides[2]) = tempX[index];
*(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 1*coord->strides[2]) = tempY[index];
*(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 2*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "acf", 3) == 0) {
*(double*)(coord->data + j*coord->strides[0] + 0*coord->strides[1] + i*coord->strides[2]) = tempX[index];
*(double*)(coord->data + j*coord->strides[0] + 1*coord->strides[1] + i*coord->strides[2]) = tempY[index];
*(double*)(coord->data + j*coord->strides[0] + 2*coord->strides[1] + i*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "fac", 3) == 0) {
*(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 0*coord->strides[2]) = tempX[index];
*(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 1*coord->strides[2]) = tempY[index];
*(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 2*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "fca", 3) == 0) {
*(double*)(coord->data + i*coord->strides[0] + 0*coord->strides[1] + j*coord->strides[2]) = tempX[index];
*(double*)(coord->data + i*coord->strides[0] + 1*coord->strides[1] + j*coord->strides[2]) = tempY[index];
*(double*)(coord->data + i*coord->strides[0] + 2*coord->strides[1] + j*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "caf", 3) == 0) {
*(double*)(coord->data + 0*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempX[index];
*(double*)(coord->data + 1*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempY[index];
*(double*)(coord->data + 2*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempZ[index];
} else if (strncasecmp(format, "cfa", 3) == 0) {
*(double*)(coord->data + 0*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempX[index];
*(double*)(coord->data + 1*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempY[index];
*(double*)(coord->data + 2*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempZ[index];
}
}
// Check if we've been interupted by the user
if (PyErr_CheckSignals() == 1) goto error;

}
// Reset trajectory
rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET);
dcd->setsread = 0;
Expand Down

0 comments on commit b7b58ad

Please sign in to comment.