diff --git a/src/mdevaluate/reader.py b/src/mdevaluate/reader.py index 2f1bb72..e8f5ce0 100755 --- a/src/mdevaluate/reader.py +++ b/src/mdevaluate/reader.py @@ -187,35 +187,52 @@ def nojump_save_filename(reader: BaseReader): return full_path_fallback -def parse_jumps(trajectory: Coordinates): - prev = trajectory[0].whole +def parse_jumps(trajectory: Coordinates, whole: bool=True, fractional_inverted: bool=True): + if whole: + prev = trajectory[0].whole + else: + prev = trajectory[0] box = prev.box + if fractional_inverted: + s_prev = prev @ np.linalg.inv(box) + SparseData = namedtuple("SparseData", ["data", "row", "col"]) jump_data = ( SparseData(data=array("b"), row=array("l"), col=array("l")), SparseData(data=array("b"), row=array("l"), col=array("l")), SparseData(data=array("b"), row=array("l"), col=array("l")), ) + for i, curr in enumerate(trajectory): if i % 500 == 0: logger.debug("Parse jumps Step: %d", i) - r3 = np.subtract(curr, prev) - delta_z = np.array(np.rint(np.divide(r3[:, 2], box[2][2])), dtype=np.int8) - r2 = np.subtract( - r3, - (np.rint(np.divide(r3[:, 2], box[2][2])))[:, np.newaxis] - * box[2][np.newaxis, :], - ) - delta_y = np.array(np.rint(np.divide(r2[:, 1], box[1][1])), dtype=np.int8) - r1 = np.subtract( - r2, - (np.rint(np.divide(r2[:, 1], box[1][1])))[:, np.newaxis] - * box[1][np.newaxis, :], - ) - delta_x = np.array(np.rint(np.divide(r1[:, 0], box[0][0])), dtype=np.int8) - delta = np.array([delta_x, delta_y, delta_z]).T - prev = curr - box = prev.box + if not fractional_inverted: + r3 = np.subtract(curr, prev) + delta_z = np.array(np.rint(np.divide(r3[:, 2], box[2][2])), dtype=np.int8) + r2 = np.subtract( + r3, + (np.rint(np.divide(r3[:, 2], box[2][2])))[:, np.newaxis] + * box[2][np.newaxis, :], + ) + delta_y = np.array(np.rint(np.divide(r2[:, 1], box[1][1])), dtype=np.int8) + r1 = np.subtract( + r2, + (np.rint(np.divide(r2[:, 1], box[1][1])))[:, np.newaxis] + * box[1][np.newaxis, :], + ) + delta_x = np.array(np.rint(np.divide(r1[:, 0], box[0][0])), dtype=np.int8) + delta = np.array([delta_x, delta_y, delta_z]).T + prev = curr + box = prev.box + else: + s_curr = curr @ np.linalg.inv(curr.box) + + ds = s_curr - s_prev + delta = np.array(np.rint(ds), dtype=np.int8) + + s_prev = s_curr + + for d in range(3): (col,) = np.where(delta[:, d] != 0) jump_data[d].col.extend(col)