From 303bdbf9d2771b87c814897c741fdf614c4650ea Mon Sep 17 00:00:00 2001 From: robrobo Date: Thu, 23 Apr 2026 13:28:42 +0200 Subject: [PATCH] parse jumps now correctly uses fractional coordinates and should work for triclinic boxes. this should lead to correct behavior in all cases except for maybe changing angles in triclinic boxes. old behavior can still be used --- src/mdevaluate/reader.py | 55 ++++++++++++++++++++++++++-------------- 1 file changed, 36 insertions(+), 19 deletions(-) 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)