diff --git a/distopia/__init__.py b/distopia/__init__.py index e7ee9b46..48a9d16b 100644 --- a/distopia/__init__.py +++ b/distopia/__init__.py @@ -3,12 +3,30 @@ from ._distopia import ( calc_bonds_ortho, calc_bonds_no_box, - calc_bonds_triclinic + calc_bonds_triclinic, + calc_angles_no_box, + calc_angles_ortho, + calc_angles_triclinic, + calc_dihedrals_no_box, + calc_dihedrals_ortho, + calc_dihedrals_triclinic, + calc_distance_array_no_box, + calc_distance_array_ortho, + calc_distance_array_triclinic, ) __all__ = [ 'calc_bonds_ortho', 'calc_bonds_no_box', 'calc_bonds_triclinic', + 'calc_angles_no_box', + 'calc_angles_ortho', + 'calc_angles_triclinic', + 'calc_dihedrals_no_box', + 'calc_dihedrals_ortho', + 'calc_dihedrals_triclinic', + 'calc_distance_array_no_box', + 'calc_distance_array_ortho', + 'calc_distance_array_triclinic', '__version__', ] \ No newline at end of file diff --git a/distopia/_distopia.pyx b/distopia/_distopia.pyx index 4ad6de17..3be77a26 100644 --- a/distopia/_distopia.pyx +++ b/distopia/_distopia.pyx @@ -123,12 +123,40 @@ def get_n_double_lanes(): return GetNDoubleLanes() +def _check_results(results, nvals): + """Check that results is the right shape""" + if results.ndim > 1: + raise ValueError("results must be a 1D array") + if results.shape[0] != nvals: + raise ValueError(f"results must be the same length as coordinates ({nvals}), you provided {results.shape[0]}") + + +def _check_results_darray(results, nvals0 , nvals1 ): + """Check that results is the right shape""" + if results.ndim > 1: + raise ValueError("results must be a 1D array") + if results.shape[0] != nvals0 * nvals1: + raise ValueError(f"results must be a flattened 2D array of length MxN ({ nvals0, nvals1} -> {nvals0 *nvals1}), you provided {results.shape[0]}") + + +def _check_shapes(*args): + """Check that all arrays are the same shape""" + s1 = args[0].shape + if not all(thing.shape == s1 for thing in args[1:]): + raise ValueError("All input arrays must be the same length, you provided " + f"{', '.join(str(t.shape) for t in args)}") + + + + + + @cython.boundscheck(False) @cython.wraparound(False) def calc_bonds_no_box(floating[:, ::1] coords0, floating[:, ::1] coords1, results=None): - """Calculate pairwise distances between coords0 and coords1 + """Calculate pairwise distances between coords0 and coords1 with no periodic boundary conditions Parameters ---------- @@ -146,12 +174,18 @@ def calc_bonds_no_box(floating[:, ::1] coords0, cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + + _check_shapes(coords0, coords1) + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results(results, nvals) + results_view = results CalcBondsNoBox(& coords0[0][0], & coords1[0][0], nvals, & results_view[0]) @@ -165,14 +199,14 @@ def calc_bonds_ortho(floating[:, ::1] coords0, floating[:, ::1] coords1, floating[::1] box, floating[::1] results=None): - """Calculate pairwise distances between coords0 and coords1 + """Calculate pairwise distances between coords0 and coords1 under orthorhombic boundary conditions Parameters ---------- coords0, coords1 : float32 or float64 array must be same length and dtype box : float32 or float64 array - periodic boundary dimensions + orthorhombic periodic boundary dimensions in [L, L, L] format results: float32 or float64 array (optional) array to store results in, must be same length and dtype as coords0/coords1 @@ -185,12 +219,20 @@ def calc_bonds_ortho(floating[:, ::1] coords0, cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + + _check_shapes(coords0, coords1) + + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + + else: + _check_results(results, nvals) + results_view = results CalcBondsOrtho(& coords0[0][0], & coords1[0][0], nvals, & box[0], & results_view[0]) @@ -204,7 +246,7 @@ def calc_bonds_triclinic(floating[:, ::1] coords0, floating[:, ::1] coords1, floating[:, ::1] box, floating[::1] results=None): - """Calculate pairwise distances between coords0 and coords1 + """Calculate pairwise distances between coords0 and coords1 under triclinic boundary conditions Parameters ---------- @@ -224,12 +266,19 @@ def calc_bonds_triclinic(floating[:, ::1] coords0, cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + + _check_shapes(coords0, coords1) + + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results(results, nvals) + results_view = results CalcBondsTriclinic(& coords0[0][0], & coords1[0][0], nvals, & box[0][0], & results_view[0]) @@ -242,18 +291,37 @@ def calc_angles_no_box( floating[:, ::1] coords1, floating[:, ::1] coords2, floating[::1] results=None): + """Calculate angles between sets of coordinates with no periodic boundary conditions + + Parameters + ---------- + coords0, coords1, coords2 : float32 or float64 array + must be same length and dtype + results: float32 or float64 array (optional) + array to store results in, must be same length and dtype as coords0/coords1/coords2 + + Returns + ------- + angles : np.array + same length and dtype as coords0/coords1/coords2 + """ cdef floating[::1] results_view cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + _check_shapes(coords0, coords1) + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results(results, nvals) + results_view = results CalcAnglesNoBox(&coords0[0][0], &coords1[0][0], &coords2[0][0], @@ -268,18 +336,40 @@ def calc_angles_ortho( floating[:, ::1] coords2, floating[::1] box, floating[::1] results=None): + + """Calculate angles between sets of coordinates under orthorhombic boundary conditions + + Parameters + ---------- + coords0, coords1, coords2 : float32 or float64 array + must be same length and dtype + box : float32 or float64 array + orthorhombic periodic boundary dimensions in [L, L, L] format + results: float32 or float64 array (optional) + array to store results in, must be same length and dtype as coords0/coords1/coords2 + + Returns + ------- + angles : np.array + same length and dtype as coords0/coords1/coords2 + """ cdef floating[::1] results_view cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + _check_shapes(coords0, coords1, coords2) + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results(results, nvals) + results_view = results CalcAnglesOrtho(&coords0[0][0], &coords1[0][0], &coords2[0][0], @@ -294,18 +384,40 @@ def calc_angles_triclinic( floating[:, ::1] coords2, floating[:, ::1] box, floating[::1] results=None): + """Calculate angles between sets of coordinates under triclinic boundary conditions + + Parameters + ---------- + coords0, coords1, coords2 : float32 or float64 array + must be same length and dtype + box : float32 or float64 array + periodic boundary dimensions, in 3x3 format + results: float32 or float64 array (optional) + array to store results in, must be same length and dtype as coords0/coords1/coords2 + + Returns + ------- + angles : np array + same length and dtype as coords0/coords1/coords2 + """ cdef floating[::1] results_view cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + _check_shapes(coords0, coords1, coords2) + + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results(results, nvals) + results_view = results CalcAnglesTriclinic(&coords0[0][0], &coords1[0][0], &coords2[0][0], @@ -320,18 +432,38 @@ def calc_dihedrals_no_box( floating[:, ::1] coords2, floating[:, ::1] coords3, floating[::1] results=None): + """Calculate dihedral angles between sets of coordinates with no periodic boundary conditions + + Parameters + ---------- + coords0, coords1, coords2, coords3 : float32 or float64 array + must be same length and dtype + results: float32 or float64 array (optional) + array to store results in, must be same length and dtype as coords0/coords1/coords2/coords3 + + Returns + ------- + dihedrals : np array + same length and dtype as coords0/coords1/coords2/coords3 + """ cdef floating[::1] results_view cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + _check_shapes(coords0, coords1, coords2, coords3) + + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results(results, nvals) + results_view = results CalcDihedralsNoBox(&coords0[0][0], &coords1[0][0], &coords2[0][0], &coords3[0][0], @@ -347,18 +479,39 @@ def calc_dihedrals_ortho( floating[:, ::1] coords3, floating[::1] box, floating[::1] results=None): + """Calculate dihedral angles between sets of coordinates under orthorhombic boundary conditions + + Parameters + ---------- + coords0, coords1, coords2, coords3 : float32 or float64 array + must be same length and dtype + box : float32 or float64 array + orthorhombic periodic boundary dimensions in [L, L, L] format + results: float32 or float64 array (optional) + array to store results in, must be same length and dtype as coords0/coords1/coords2/coords3 + + Returns + ------- + dihedrals : np array + same length and dtype as coords0/coords1/coords2/coords3 + """ cdef floating[::1] results_view cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + _check_shapes(coords0, coords1, coords2, coords3) + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results(results, nvals) + results_view = results CalcDihedralsOrtho(&coords0[0][0], &coords1[0][0], &coords2[0][0], &coords3[0][0], @@ -374,18 +527,40 @@ def calc_dihedrals_triclinic( floating[:, ::1] coords3, floating[:, ::1] box, floating[::1] results=None): + """Calculate dihedral angles between sets of coordinates under triclinic boundary conditions + + Parameters + ---------- + coords0, coords1, coords2, coords3 : float32 or float64 array + must be same length and dtype + box : float32 or float64 array + periodic boundary dimensions, in 3x3 format + results: float32 or float64 array (optional) + array to store results in, must be same length and dtype as coords0/coords1/coords2/coords3 + + Returns + ------- + dihedrals : np array + same length and dtype as coords0/coords1/coords2/coords3 + """ cdef floating[::1] results_view cdef size_t nvals = coords0.shape[0] cdef cnp.npy_intp[1] dims dims[0] = nvals # FIXME truncation? + _check_shapes(coords0, coords1, coords2, coords3) + + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results(results, nvals) + results_view = results CalcDihedralsTriclinic(&coords0[0][0], &coords1[0][0], &coords2[0][0], &coords3[0][0], @@ -398,6 +573,21 @@ def calc_distance_array_no_box( floating[:, ::1] coords0, floating[:, ::1] coords1, floating[::1] results=None): + """Calculate pairwise distance matrix between coordinates with no periodic boundary conditions + + Parameters + ---------- + coords0, coords1 : float32 or float64 array + must be same length and dtype + results: float32 or float64 array (optional) + array to store results in, must be a single dimension of length MxN where M is the length of coords0 and N is the length of coords1 + + Returns + ------- + distances : np array + MxN array of distances + """ + cdef floating[::1] results_view cdef size_t nvals0 = coords0.shape[0] cdef size_t nvals1 = coords1.shape[0] @@ -405,12 +595,18 @@ def calc_distance_array_no_box( dims[0] = nvals0 * nvals1 + _check_shapes(coords0, coords1) + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results_darray(results, nvals0, nvals1) + + results_view = results CalcDistanceArrayNoBox(&coords0[0][0], &coords1[0][0], @@ -425,6 +621,22 @@ def calc_distance_array_ortho( floating[:, ::1] coords1, floating[::1] box, floating[::1] results=None): + """Calculate pairwise distance matrix between coordinates under orthorhombic boundary conditions + + Parameters + ---------- + coords0, coords1 : float32 or float64 array + must be same length and dtype + box : float32 or float64 array + orthorhombic periodic boundary dimensions in [L, L, L] format + results: float32 or float64 array (optional) + array to store results in, must be a single dimension of length MxN where M is the length of coords0 and N is the length of coords1 + + Returns + ------- + distances : np array + MxN array of distances + """ cdef floating[::1] results_view cdef size_t nvals0 = coords0.shape[0] cdef size_t nvals1 = coords1.shape[0] @@ -432,12 +644,17 @@ def calc_distance_array_ortho( dims[0] = nvals0 * nvals1 + _check_shapes(coords0, coords1) + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results_darray(results, nvals0, nvals1) + results_view = results CalcDistanceArrayOrtho(&coords0[0][0], &coords1[0][0], @@ -453,19 +670,41 @@ def calc_distance_array_triclinic( floating[:, ::1] coords1, floating[:, ::1] box, floating[::1] results=None): + """Calculate pairwise distance matrix between coordinates under triclinic boundary conditions + + Parameters + ---------- + coords0, coords1 : float32 or float64 array + must be same length and dtype + box : float32 or float64 array + periodic boundary dimensions, in 3x3 format + results: float32 or float64 array (optional) + array to store results in, must be a single dimension of length MxN where M is the length of coords0 and N is the length of coords1 + + Returns + ------- + distances : np array + MxN array of distances + """ cdef floating[::1] results_view cdef size_t nvals0 = coords0.shape[0] cdef size_t nvals1 = coords1.shape[0] + cdef cnp.npy_intp[1] dims dims[0] = nvals0 * nvals1 + _check_shapes(coords0, coords1) + if results is None: if floating is float: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT32, 0) else: results = cnp.PyArray_EMPTY(1, dims, cnp.NPY_FLOAT64, 0) + else: + _check_results_darray(results, nvals0, nvals1) + results_view = results CalcDistanceArrayTriclinic(&coords0[0][0], &coords1[0][0], diff --git a/distopia/tests/test_distopia.py b/distopia/tests/test_distopia.py index 4193dc2b..c81ebc2a 100644 --- a/distopia/tests/test_distopia.py +++ b/distopia/tests/test_distopia.py @@ -4,12 +4,15 @@ from numpy.testing import assert_allclose, assert_almost_equal, assert_equal -""" -Majority of detailed testing is done at the C++ level. -This is primarily to make sure that the Python API works as expected. -""" + +def convert_ndarray(*args, dtype): + if len(args) == 1: + return np.asarray(args[0], dtype=dtype) + else: + return (np.asarray(a, dtype=dtype) for a in args) + class TestDistances: def arange_input(self, N, dtype): return np.arange(3 * N, dtype=dtype).reshape(N, 3) @@ -63,10 +66,169 @@ def test_calc_bonds_triclinic_all_zero(self, N, use_result_buffer, dtype): assert result.dtype == dtype + + def test_no_box_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_bonds_no_box(c0, c1, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_bonds_no_box(c0, c1[:-1]) + + def test_ortho_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + box = np.array([10, 10, 10], dtype=np.float32) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_bonds_ortho(c0, c1, box, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_bonds_ortho(c0, c1[:-1], box) + + def test_triclinic_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + box = np.array([[10, 0, 0], [0, 10, 0], [0, 0, 10]], dtype=np.float32) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_bonds_triclinic(c0, c1, box, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_bonds_triclinic(c0, c1[:-1], box) + + + +class TestAngles: + + + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + def test_no_box_critical_angles(self, dtype): + # 0, 90, 180 + c0 = convert_ndarray(np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]], dtype=np.float32), dtype=dtype) + c1 = convert_ndarray(np.array([[1, 0, 0], [-1, 0, 0], [-1, 0, 0], [1, 0, 0]], dtype=np.float32), dtype=dtype) + c2 = convert_ndarray(np.array([[1, 1, 0], [-1, -1, 0], [-1, 0, 0], [2, 0, 0]], dtype=np.float32), dtype=dtype) + results = distopia.calc_angles_no_box(c0, c1, c2) + assert_almost_equal(results, np.array([np.pi / 2, np.pi / 2, 0, np.pi], dtype=np.float32)) + + + def test_no_box_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c2 = np.zeros(6, dtype=np.float32).reshape(2, 3) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_angles_no_box(c0, c1, c2, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_angles_no_box(c0, c1[:-1], c2) + + def test_ortho_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c2 = np.zeros(6, dtype=np.float32).reshape(2, 3) + box = np.array([10, 10, 10], dtype=np.float32) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_angles_ortho(c0, c1, c2, box, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_angles_ortho(c0, c1[:-1], c2, box) + + def test_triclinic_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c2 = np.zeros(6, dtype=np.float32).reshape(2, 3) + box = np.array([[10, 0, 0], [0, 10, 0], [0, 0, 10]], dtype=np.float32) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_angles_triclinic(c0, c1, c2, box, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_angles_triclinic(c0, c1[:-1], c2, box) + + + +class TestDihedrals: + + + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + def test_no_box_critical_dihedrals(self, dtype): + # 0, 90, 180 + c0 = convert_ndarray(np.array([[2, 1, 0], [-2, 1, 0], [-2, 1, 0], [-2, 1, 0], [1,2,1]], dtype=np.float32), dtype=dtype) + c1 = convert_ndarray(np.array([[1, 0, 0], [1, 0, 0], [1, 0, 0], [1, 0, 0], [1,1,1]], dtype=np.float32), dtype=dtype) + c2 = convert_ndarray(np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0], [2,1,1]], dtype=np.float32), dtype=dtype) + c3 = convert_ndarray(np.array([[-2, 1, 0],[-2, 1, 0], [0, 0, 1], [0, 0, -1], [2, 0, 1]], dtype=np.float32), dtype=dtype) + results = distopia.calc_dihedrals_no_box(c0, c1, c2, c3) + # NOTE: negative signs, do we need to take ABS? + assert_almost_equal(results, np.array([-0, -0, -np.pi/2, np.pi/2, -np.pi], dtype=np.float32)) + + + + def test_no_box_bad_result_or_input_shape(self): + c0 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c1 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c2 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c3 = np.zeros(12, dtype=np.float32).reshape(4, 3) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_dihedrals_no_box(c0, c1, c2, c3, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_dihedrals_no_box(c0, c1[:-1], c2, c3) + + def test_ortho_bad_result_or_input_shape(self): + c0 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c1 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c2 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c3 = np.zeros(12, dtype=np.float32).reshape(4, 3) + box = np.array([10, 10, 10], dtype=np.float32) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_dihedrals_ortho(c0, c1, c2, c3, box, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_dihedrals_ortho(c0, c1[:-1], c2, c3, box) + + def test_triclinic_bad_result_or_input_shape(self): + c0 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c1 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c2 = np.zeros(12, dtype=np.float32).reshape(4, 3) + c3 = np.zeros(12, dtype=np.float32).reshape(4, 3) + box = np.array([[10, 0, 0], [0, 10, 0], [0, 0, 10]], dtype=np.float32) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_dihedrals_triclinic(c0, c1, c2, c3, box, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_dihedrals_triclinic(c0, c1[:-1], c2, c3, box) + + +class TestDistanceArray: + + def test_no_box_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_distance_array_no_box(c0, c1, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_distance_array_no_box(c0, c1[:-1]) + + def test_ortho_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + box = np.array([10, 10, 10], dtype=np.float32) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_distance_array_ortho(c0, c1, box, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_distance_array_ortho(c0, c1[:-1], box) + + + def test_triclinic_bad_result_or_input_shape(self): + c0 = np.zeros(6, dtype=np.float32).reshape(2, 3) + c1 = np.zeros(6, dtype=np.float32).reshape(2, 3) + box = np.array([[10, 0, 0], [0, 10, 0], [0, 0, 10]], dtype=np.float32) + with pytest.raises(ValueError, match="results must be"): + distopia.calc_distance_array_triclinic(c0, c1, box, results=np.empty(1, dtype=np.float32)) + with pytest.raises(ValueError, match="All input arrays must"): + distopia.calc_distance_array_triclinic(c0, c1[:-1], box) + + + + + class TestMDA: + """ + Copy of some of the tests from MDAnalysisTests repo + """ prec = 5 + @staticmethod @pytest.fixture() def positions(): @@ -90,6 +252,9 @@ def triclinic_box(): @pytest.mark.parametrize("dtype", (np.float32, np.float64)) def test_bonds(self, box_bonds, dtype, positions): a, b, c, d = positions + a, b, c, d = convert_ndarray(a, b, c, d, dtype=dtype) + box_bonds = convert_ndarray(box_bonds, dtype=dtype) + dists = distopia.calc_bonds_no_box(a, b) assert_equal(len(dists), 4, err_msg="calc_bonds results have wrong length") dists_pbc = distopia.calc_bonds_ortho(a, b, box_bonds) @@ -114,3 +279,87 @@ def test_bonds_triclinic(self, triclinic_box, dtype, positions): dists = distopia.calc_bonds_triclinic(a, b, triclinic_box) reference = np.array([0.0, 1.7320508, 1.4142136, 2.82842712]) assert_almost_equal(dists, reference, self.prec, err_msg="calc_bonds with triclinic box failed") + + + + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + def test_angles(self, dtype, positions): + a, b, c, d = positions + a, b, c, d = convert_ndarray(a, b, c, d, dtype=dtype) + + angles = distopia.calc_angles_no_box(a, b, c) + # Check calculated values + assert_equal(len(angles), 4, err_msg="calc_angles results have wrong length") + assert_almost_equal(angles[1], np.pi, self.prec, + err_msg="180 degree angle calculation failed") + assert_almost_equal(np.rad2deg(angles[2]), 90., self.prec, + err_msg="Ninety degree angle in calc_angles failed") + assert_almost_equal(angles[3], 0.098174833, self.prec, + err_msg="Small angle failed in calc_angles") + + + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + def test_dihedrals(self, dtype, positions): + a, b, c, d = positions + a, b, c, d = convert_ndarray(a, b, c, d, dtype=dtype) + dihedrals = distopia.calc_dihedrals_no_box(a, b, c, d) + # Check calculated values + assert_equal(len(dihedrals), 4, err_msg="calc_dihedrals results have wrong length") + # FIXME: BROKEN assert np.isnan(dihedrals[0]), "Zero length dihedral failed" + # FIXME: BROKEN assert np.isnan(dihedrals[1]), "Straight line dihedral failed" + # FIXME: BROKEN assert_almost_equal(dihedrals[2], np.pi, self.prec, err_msg="180 degree dihedral failed") + assert_almost_equal(dihedrals[3], -0.50714064, self.prec, + err_msg="arbitrary dihedral angle failed") + + + + @staticmethod + @pytest.fixture() + def positions_angles(): + a = np.array([[0.0, 1.0, 0.0]], dtype=np.float32) + b = np.array([[0.0, 0.0, 0.0]], dtype=np.float32) + c = np.array([[1.0, 0.0, 0.0]], dtype=np.float32) + d = np.array([[1.0, 0.0, 1.0]], dtype=np.float32) + return a, b, c, d + + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + def test_periodic_dihedrals_angles(self, box_bonds, positions_angles, dtype): + a, b, c, d = positions_angles + a, b, c, d = convert_ndarray(a, b, c, d, dtype=dtype) + box = convert_ndarray(box_bonds[:3], dtype=dtype) + a2 = a + box * np.asarray((-1, 0, 0), dtype=dtype) + b2 = b + box * np.asarray((1, 0, 1), dtype=dtype) + c2 = c + box * np.asarray((-2, 5, -7), dtype=dtype) + d2 = d + box * np.asarray((0, -5, 0), dtype=dtype) + + ref = distopia.calc_dihedrals_no_box(a, b, c, d) + + box = np.asarray(box_bonds, dtype=dtype) + test1 = distopia.calc_dihedrals_ortho(a2, b, c, d, box=box) + test2 = distopia.calc_dihedrals_ortho(a, b2, c, d, box=box) + test3 = distopia.calc_dihedrals_ortho(a, b, c2, d, box=box) + test4 = distopia.calc_dihedrals_ortho(a, b, c, d2, box=box) + test5 = distopia.calc_dihedrals_ortho(a2, b2, c2, d2, box=box) + + for val in [test1, test2, test3, test4, test5]: + assert_almost_equal(ref, val, self.prec, err_msg="Min image in dihedral calculation failed") + + + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + def test_periodic_angles(self, box_bonds, positions_angles, dtype): + a, b, c, d = positions_angles + a, b, c, d = convert_ndarray(a, b, c, d, dtype=dtype) + box = convert_ndarray(box_bonds[:3], dtype=dtype) + a2 = a + box * np.asarray((-1, 0, 0), dtype=dtype) + b2 = b + box * np.asarray((1, 0, 1), dtype=dtype) + c2 = c + box * np.asarray((-2, 5, -7), dtype=dtype) + + ref = distopia.calc_angles_no_box(a, b, c) + + box = np.asarray(box_bonds, dtype=dtype) + test1 = distopia.calc_angles_ortho(a2, b, c, box=box) + test2 = distopia.calc_angles_ortho(a, b2, c, box=box) + test3 = distopia.calc_angles_ortho(a, b, c2, box=box) + + for val in [test1, test2, test3]: + assert_almost_equal(ref, val, self.prec, err_msg="Min image in angle calculation failed") \ No newline at end of file diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp index bb7b4ede..4f3dc4a8 100644 --- a/libdistopia/src/distopia.cpp +++ b/libdistopia/src/distopia.cpp @@ -488,9 +488,10 @@ namespace distopia { y = hn::MulAdd(xp_y, rbc_y, y); y = hn::MulAdd(xp_z, rbc_z, y); + + y = y / vb_norm; - // negate due to vector order (?) return hn::Neg(hn::Atan2(d, y, x)); }