From 51eb3d508038352924b4577cf327bb3552243bd5 Mon Sep 17 00:00:00 2001 From: Michelle Richer Date: Tue, 30 Apr 2024 15:11:34 -0400 Subject: [PATCH 1/4] Add transform= keyword to cbasis integrals --- gbasis/integrals/libcint.py | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/gbasis/integrals/libcint.py b/gbasis/integrals/libcint.py index 3fb368d0..d3ae1146 100644 --- a/gbasis/integrals/libcint.py +++ b/gbasis/integrals/libcint.py @@ -670,7 +670,7 @@ def make_int1e( ) # Make instance-bound integral method - def int1e(notation="physicist", origin=None, inv_origin=None): + def int1e(notation="physicist", transform=None, origin=None, inv_origin=None): # Handle ``notation`` argument if notation not in ("physicist", "chemist"): raise ValueError("``notation`` must be one of 'physicist' or 'chemist'") @@ -755,11 +755,17 @@ def int1e(notation="physicist", origin=None, inv_origin=None): # Apply permutation out = out[self._permutations, :][:, self._permutations] - # Return normalized integrals + # Normalize integrals if self.coord_type == "cartesian": - return np.einsum(norm_einsum, self._ovlp_minhalf, self._ovlp_minhalf, out) - else: - return out + out = np.einsum(norm_einsum, self._ovlp_minhalf, self._ovlp_minhalf, out) + + # Apply transformation + if transform is not None: + out = np.tensordot(transform, out, (1, 0)) + out = np.tensordot(transform, out, (1, 1)) + out = np.swapaxes(out, 0, 1) + + return out # Return instance-bound integral method return int1e @@ -823,7 +829,7 @@ def make_int2e( ) # Make instance-bound integral method - def int2e(notation="physicist", origin=None, inv_origin=None): + def int2e(notation="physicist", transform=None, origin=None, inv_origin=None): # Handle ``notation`` argument if notation == "physicist": physicist = True @@ -962,9 +968,9 @@ def int2e(notation="physicist", origin=None, inv_origin=None): out = out[:, :, self._permutations] out = out[:, :, :, self._permutations] - # Return normalized integrals + # Normalize integrals if self.coord_type == "cartesian": - return np.einsum( + out = np.einsum( norm_einsum, self._ovlp_minhalf, self._ovlp_minhalf, @@ -972,8 +978,16 @@ def int2e(notation="physicist", origin=None, inv_origin=None): self._ovlp_minhalf, out, ) - else: - return out + + # Apply transformation + if transform is not None: + out = np.tensordot(transform, out, (1, 0)) + out = np.tensordot(transform, out, (1, 1)) + out = np.tensordot(transform, out, (1, 2)) + out = np.tensordot(transform, out, (1, 3)) + out = np.swapaxes(np.swapaxes(out, 0, 3), 1, 2) + + return out # Return instance-bound integral method return int2e From 7d935404cbf7cb5db4ffa20ca12792cba28574c9 Mon Sep 17 00:00:00 2001 From: Michelle Richer Date: Wed, 1 May 2024 11:34:25 -0400 Subject: [PATCH 2/4] Expose transform= keyword --- gbasis/integrals/libcint.py | 100 +++++++++++++++++++++++++++--------- 1 file changed, 77 insertions(+), 23 deletions(-) diff --git a/gbasis/integrals/libcint.py b/gbasis/integrals/libcint.py index d3ae1146..e642a751 100644 --- a/gbasis/integrals/libcint.py +++ b/gbasis/integrals/libcint.py @@ -958,10 +958,6 @@ def int2e(notation="physicist", transform=None, origin=None, inv_origin=None): if constant is not None: out *= constant - # Transpose integrals in `out` array to proper notation - if physicist: - out = out.transpose(0, 2, 1, 3) - # Apply permutation out = out[self._permutations] out = out[:, self._permutations] @@ -979,6 +975,10 @@ def int2e(notation="physicist", transform=None, origin=None, inv_origin=None): out, ) + # Transpose integrals in `out` array to proper notation + if physicist: + out = out.transpose(0, 2, 1, 3) + # Apply transformation if transform is not None: out = np.tensordot(transform, out, (1, 0)) @@ -992,7 +992,7 @@ def int2e(notation="physicist", transform=None, origin=None, inv_origin=None): # Return instance-bound integral method return int2e - def overlap_integral(self, notation="physicist"): + def overlap_integral(self, notation="physicist", transform=None): r""" Compute the overlap integrals. @@ -1000,6 +1000,12 @@ def overlap_integral(self, notation="physicist"): ---------- notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1007,9 +1013,9 @@ def overlap_integral(self, notation="physicist"): Integral array. """ - return self._ovlp(notation=notation) + return self._ovlp(notation=notation, transform=transform) - def kinetic_energy_integral(self, notation="physicist"): + def kinetic_energy_integral(self, notation="physicist", transform=None): r""" Compute the kinetic energy integrals. @@ -1017,6 +1023,12 @@ def kinetic_energy_integral(self, notation="physicist"): ---------- notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1024,9 +1036,9 @@ def kinetic_energy_integral(self, notation="physicist"): Integral array. """ - return self._kin(notation=notation) + return self._kin(notation=notation, transform=transform) - def nuclear_attraction_integral(self, notation="physicist"): + def nuclear_attraction_integral(self, notation="physicist", transform=None): r""" Compute the nuclear attraction integrals. @@ -1034,6 +1046,12 @@ def nuclear_attraction_integral(self, notation="physicist"): ---------- notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1041,9 +1059,9 @@ def nuclear_attraction_integral(self, notation="physicist"): Integral array. """ - return self._nuc(notation=notation) + return self._nuc(notation=notation, transform=transform) - def electron_repulsion_integral(self, notation="physicist"): + def electron_repulsion_integral(self, notation="physicist", transform=None): r""" Compute the electron repulsion integrals. @@ -1051,6 +1069,12 @@ def electron_repulsion_integral(self, notation="physicist"): ---------- notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1058,9 +1082,9 @@ def electron_repulsion_integral(self, notation="physicist"): Integral array. """ - return self._eri(notation=notation) + return self._eri(notation=notation, transform=transform) - def r_inv_integral(self, origin=None, notation="physicist"): + def r_inv_integral(self, origin=None, notation="physicist", transform=None): r""" Compute the :math:`1/\left|\mathbf{r} - \mathbf{R}_\text{inv}\right|` integrals. @@ -1070,6 +1094,12 @@ def r_inv_integral(self, origin=None, notation="physicist"): Origin about which to evaluate integrals. notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1077,9 +1107,9 @@ def r_inv_integral(self, origin=None, notation="physicist"): Integral array. """ - return self._rinv(inv_origin=origin, notation=notation) + return self._rinv(inv_origin=origin, notation=notation, transform=transform) - def momentum_integral(self, origin=None, notation="physicist"): + def momentum_integral(self, origin=None, notation="physicist", transform=None): r""" Compute the momentum integrals. @@ -1089,6 +1119,12 @@ def momentum_integral(self, origin=None, notation="physicist"): Origin about which to evaluate integrals. notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1096,9 +1132,9 @@ def momentum_integral(self, origin=None, notation="physicist"): Integral array. """ - return self._mom(origin=origin, notation=notation) + return self._mom(origin=origin, notation=notation, transform=transform) - def angular_momentum_integral(self, origin=None, notation="physicist"): + def angular_momentum_integral(self, origin=None, notation="physicist", transform=None): r""" Compute the angular momentum integrals. @@ -1108,6 +1144,12 @@ def angular_momentum_integral(self, origin=None, notation="physicist"): Origin about which to evaluate integrals. notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1116,9 +1158,9 @@ def angular_momentum_integral(self, origin=None, notation="physicist"): """ raise NotImplementedError("Angular momentum integral doesn't work; see Issue #149") - # return self._amom(origin=origin, notation=notation) + # return self._amom(origin=origin, notation=notation, transform=transform) - def point_charge_integral(self, point_coords, point_charges, notation="physicist"): + def point_charge_integral(self, point_coords, point_charges, notation="physicist", transform=None): r""" Compute the point charge integrals. @@ -1130,6 +1172,12 @@ def point_charge_integral(self, point_coords, point_charges, notation="physicist Charges of point charges. notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1141,13 +1189,13 @@ def point_charge_integral(self, point_coords, point_charges, notation="physicist out = np.zeros((self.nbfn, self.nbfn, len(point_charges)), dtype=c_double, order="F") # Compute 1/|r - r_{inv}| for each charge for icharge, (coord, charge) in enumerate(zip(point_coords, point_charges)): - val = self._rinv(inv_origin=coord, notation=notation) + val = self._rinv(inv_origin=coord, notation=notation, transform=transform) val *= -charge out[:, :, icharge] = val # Return integrals in `out` array return out - def moment_integral(self, orders, origin=None, notation="physicist"): + def moment_integral(self, orders, origin=None, notation="physicist", transform=None): r""" Compute the moment integrals. @@ -1159,6 +1207,12 @@ def moment_integral(self, orders, origin=None, notation="physicist"): Origin about which to evaluate integrals. notation : ("physicist" | "chemist"), default="physicist" Axis order convention. + transform : np.ndarray(K, K_cont) + Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear + combinations of contractions (e.g. MO). + Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` + and index 0 of the array for contractions. + Default is no transformation. Returns ------- @@ -1178,9 +1232,9 @@ def moment_integral(self, orders, origin=None, notation="physicist"): try: for i, order in enumerate(orders): if sum(order) == 0: - out[:, :, i] = self._ovlp(notation=notation) + out[:, :, i] = self._ovlp(notation=notation, transform=transform) else: - out[:, :, i] = self._moments[tuple(order)](origin=origin, notation=notation) + out[:, :, i] = self._moments[tuple(order)](origin=origin, notation=notation, transform=transform) except KeyError: raise ValueError( "Invalid order; can use up to order 4 for any XYZ component," From d51fec68a0cde721fbab03945dcb3c16f1e08342 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 1 May 2024 15:35:06 +0000 Subject: [PATCH 3/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- gbasis/integrals/libcint.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/gbasis/integrals/libcint.py b/gbasis/integrals/libcint.py index e642a751..e0c74565 100644 --- a/gbasis/integrals/libcint.py +++ b/gbasis/integrals/libcint.py @@ -1160,7 +1160,9 @@ def angular_momentum_integral(self, origin=None, notation="physicist", transform raise NotImplementedError("Angular momentum integral doesn't work; see Issue #149") # return self._amom(origin=origin, notation=notation, transform=transform) - def point_charge_integral(self, point_coords, point_charges, notation="physicist", transform=None): + def point_charge_integral( + self, point_coords, point_charges, notation="physicist", transform=None + ): r""" Compute the point charge integrals. @@ -1234,7 +1236,9 @@ def moment_integral(self, orders, origin=None, notation="physicist", transform=N if sum(order) == 0: out[:, :, i] = self._ovlp(notation=notation, transform=transform) else: - out[:, :, i] = self._moments[tuple(order)](origin=origin, notation=notation, transform=transform) + out[:, :, i] = self._moments[tuple(order)]( + origin=origin, notation=notation, transform=transform + ) except KeyError: raise ValueError( "Invalid order; can use up to order 4 for any XYZ component," From ecdd01cde021575cec36980ca58b097702020e41 Mon Sep 17 00:00:00 2001 From: leila-pujal Date: Thu, 2 May 2024 12:50:54 -0400 Subject: [PATCH 4/4] Add tests for CBasis transform argument --- tests/test_libcint.py | 125 ++++++++++++++++++++++++++++++------------ 1 file changed, 91 insertions(+), 34 deletions(-) diff --git a/tests/test_libcint.py b/tests/test_libcint.py index e3656cca..0a8ff5d2 100644 --- a/tests/test_libcint.py +++ b/tests/test_libcint.py @@ -169,6 +169,11 @@ def test_integral(basis, atsyms, atcoords, coord_type, integral): pytest.param("h2o_hf_ccpv5z_sph.fchk", ["O", "H", "H"], "Spherical", id="h2o_sph"), ] +TEST_COORD_TRANSFORM = [ + pytest.param(False, id="no-transform"), + pytest.param(True, id="transform"), +] + TEST_INTEGRALS_IODATA = [ pytest.param("overlap", id="Overlap"), pytest.param("kinetic_energy", id="KineticEnergy"), @@ -183,8 +188,9 @@ def test_integral(basis, atsyms, atcoords, coord_type, integral): @pytest.mark.skipif(os.path.exists(f"{os.path.dirname(os.path.abspath(__file__).split('tests')[0])}/build") == False, reason="Libcint build not found") @pytest.mark.parametrize("fname, elements, coord_type", TEST_SYSTEMS_IODATA) +@pytest.mark.parametrize("transform", TEST_COORD_TRANSFORM) @pytest.mark.parametrize("integral", TEST_INTEGRALS_IODATA) -def test_integral_iodata(fname, elements, coord_type, integral): +def test_integral_iodata(fname, elements, coord_type, integral, transform): pytest.importorskip("iodata") from iodata import load_one from gbasis.integrals.libcint import ELEMENTS, LIBCINT, CBasis @@ -197,22 +203,41 @@ def test_integral_iodata(fname, elements, coord_type, integral): lc_basis = CBasis(py_basis, elements, mol.atcoords, coord_type=coord_type) if integral == "overlap": - py_int = overlap_integral(py_basis) - npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) - lc_int = lc_basis.overlap_integral() - npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + if transform: + py_int = overlap_integral(py_basis, transform=mol.mo.coeffs.T) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + lc_int = lc_basis.overlap_integral(transform=mol.mo.coeffs.T) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + else: + py_int = overlap_integral(py_basis) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + lc_int = lc_basis.overlap_integral() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) elif integral == "kinetic_energy": - py_int = kinetic_energy_integral(py_basis) - npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) - lc_int = lc_basis.kinetic_energy_integral() - npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + if transform: + py_int = kinetic_energy_integral(py_basis, transform=mol.mo.coeffs.T) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + lc_int = lc_basis.kinetic_energy_integral(transform=mol.mo.coeffs.T) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + else: + py_int = kinetic_energy_integral(py_basis) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + lc_int = lc_basis.kinetic_energy_integral() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) elif integral == "nuclear_attraction": - py_int = nuclear_electron_attraction_integral(py_basis, mol.atcoords, mol.atnums) - npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) - lc_int = lc_basis.nuclear_attraction_integral() - npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + if transform: + py_int = nuclear_electron_attraction_integral(py_basis, mol.atcoords, + mol.atnums, transform=mol.mo.coeffs.T) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + lc_int = lc_basis.nuclear_attraction_integral(transform=mol.mo.coeffs.T) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + else: + py_int = nuclear_electron_attraction_integral(py_basis, mol.atcoords, mol.atnums) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + lc_int = lc_basis.nuclear_attraction_integral() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) elif integral == "angular_momentum": # py_int = angular_momentum_integral(py_basis) @@ -223,29 +248,55 @@ def test_integral_iodata(fname, elements, coord_type, integral): return elif integral == "momentum": - py_int = momentum_integral(py_basis) - npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, 3)) - lc_int = lc_basis.momentum_integral(origin=np.zeros(3)) - npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, 3)) + if transform: + py_int = momentum_integral(py_basis, transform=mol.mo.coeffs.T) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, 3)) + lc_int = lc_basis.momentum_integral(origin=np.zeros(3), transform=mol.mo.coeffs.T) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, 3)) + else: + py_int = momentum_integral(py_basis) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, 3)) + lc_int = lc_basis.momentum_integral(origin=np.zeros(3)) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, 3)) elif integral == "electron_repulsion": - py_int = electron_repulsion_integral(py_basis) - npt.assert_array_equal( - py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn) - ) - lc_int = lc_basis.electron_repulsion_integral() - npt.assert_array_equal( - lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn) - ) + if transform: + py_int = electron_repulsion_integral(py_basis, transform=mol.mo.coeffs.T) + npt.assert_array_equal( + py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn) + ) + lc_int = lc_basis.electron_repulsion_integral(transform=mol.mo.coeffs.T) + npt.assert_array_equal( + lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn) + ) + else: + py_int = electron_repulsion_integral(py_basis) + npt.assert_array_equal( + py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn) + ) + lc_int = lc_basis.electron_repulsion_integral() + npt.assert_array_equal( + lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn) + ) elif integral == "point_charge": charge_coords = np.asarray([[2.0, 2.0, 2.0], [-3.0, -3.0, -3.0], [-1.0, 2.0, -3.0]]) charges = np.asarray([1.0, 0.666, -3.1415926]) - for i in range(1, len(charges) + 1): - py_int = point_charge_integral(py_basis, charge_coords[:i], charges[:i]) - npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) - lc_int = lc_basis.point_charge_integral(charge_coords[:i], charges[:i]) - npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) + if transform: + for i in range(1, len(charges) + 1): + py_int = point_charge_integral(py_basis, charge_coords[:i], + charges[:i], transform=mol.mo.coeffs.T) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) + lc_int = lc_basis.point_charge_integral(charge_coords[:i], + charges[:i], transform=mol.mo.coeffs.T) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) + + else: + for i in range(1, len(charges) + 1): + py_int = point_charge_integral(py_basis, charge_coords[:i], charges[:i]) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) + lc_int = lc_basis.point_charge_integral(charge_coords[:i], charges[:i]) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) elif integral == "moment": origin = np.zeros(3) @@ -263,10 +314,16 @@ def test_integral_iodata(fname, elements, coord_type, integral): [0, 1, 1], ] ) - py_int = moment_integral(py_basis, origin, orders) - npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) - lc_int = lc_basis.moment_integral(orders, origin=origin) - npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) + if transform: + py_int = moment_integral(py_basis, origin, orders, transform=mol.mo.coeffs.T) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) + lc_int = lc_basis.moment_integral(orders, origin=origin, transform=mol.mo.coeffs.T) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) + else: + py_int = moment_integral(py_basis, origin, orders) + npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) + lc_int = lc_basis.moment_integral(orders, origin=origin) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) else: raise ValueError("Invalid integral name '{integral}' passed")