diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d243677..7705735b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -56,6 +56,7 @@ set_target_properties(gbasis_ext PROPERTIES # Expose Python and NumPy include directories to C extension module. target_include_directories(gbasis_ext PRIVATE + ${cint_SOURCE_DIR}/include ${Python_INCLUDE_DIRS} ${Python_NumPy_INCLUDE_DIRS}) # Link C extension module with libcint/qcint. diff --git a/gbasis/integrals/libcint.py b/gbasis/integrals/libcint.py index ee2ad0c7..64a0f1d2 100644 --- a/gbasis/integrals/libcint.py +++ b/gbasis/integrals/libcint.py @@ -2,13 +2,16 @@ Python C-API bindings for ``libcint`` GTO integrals library. """ + from gbasis.integrals.lib import libcint_bindings from contextlib import contextmanager -from ctypes import CDLL, POINTER, Structure, cdll, byref, c_int, c_double, c_void_p +# to Remove 1 +from ctypes import CDLL, POINTER, Structure, cdll, byref, c_int, c_double -from itertools import chain +# to Remove 2 +# from itertools import chain from operator import attrgetter @@ -20,7 +23,8 @@ from scipy.special import factorial -from gbasis.utils import factorial2 +# to Remove 3 +# from gbasis.utils import factorial2 __all__ = [ @@ -194,6 +198,7 @@ def from_param(cls, obj): class PairData(Structure): r"""``libcint`` ``PairData`` class.""" + _fields_ = [ ("rij", c_double * 3), ("eij", c_double), @@ -203,6 +208,7 @@ class PairData(Structure): class CINTOpt(Structure): r"""``libcint`` ``CINTOpt`` class.""" + _fields_ = [ ("index_xyz_array", POINTER(POINTER(c_int))), ("non0ctr", POINTER(POINTER(c_int))), @@ -220,6 +226,7 @@ class _LibCInt: """ import platform as _platform + _lib_dir = Path(__file__).parent / "lib" _system = _platform.system() if _system == "Darwin": @@ -309,7 +316,12 @@ def __getattr__(self, attr): # opt POINTER(CINTOpt), # cache - ndptr(enable_null=True, dtype=c_double, ndim=1, flags=("C_CONTIGUOUS",)), + ndptr( + enable_null=True, + dtype=c_double, + ndim=1, + flags=("C_CONTIGUOUS",), + ), ) cfunc.restype = c_int @@ -396,10 +408,12 @@ class CBasis: Methods ------- - make_int1e(self, func_name, components=tuple(), constant=None, is_complex=False, origin=False, inv_origin=False) - Make an instance-bound 1-electron integral method from a ``libcint`` function. - make_int2e(self, func_name, components=tuple(), constant=None, is_complex=False, origin=False, inv_origin=False) - Make an instance-bound 2-electron integral method from a ``libcint`` function. + make_int1e(self, func_name, components=tuple(), constant=None, + is_complex=False, origin=False, inv_origin=False) + Make an instance-bound 1-electron integral method from a ``libcint`` function. + make_int2e(self, func_name, components=tuple(), constant=None, + is_complex=False, origin=False, inv_origin=False) + Make an instance-bound 2-electron integral method from a ``libcint`` function. overlap(self) Compute the overlap integrals. kinetic_energy(self) @@ -575,7 +589,7 @@ def __init__(self, basis, atnums, atcoords, coord_type="spherical"): "int1e_ipovlp", components=(3,), constant=-1j, is_complex=True, origin=True ) # self._amom = self.make_int1e( - # "int1e_rxp", components=(3,), constant=-1j, is_complex=True, origin=True + # "int1e_rxp", components=(3,), constant=-1j, is_complex=True, origin=True # ) self._d_ovlp = self.make_int1e("int1e_ipovlp", components=(3,)) self._d_kin = self.make_int1e("int1e_ipkin", components=(3,)) @@ -585,29 +599,29 @@ def __init__(self, basis, atnums, atcoords, coord_type="spherical"): self._moments = {} # Order 1: int1e_r has 3 components (x, y, z) _r = self.make_int1e("int1e_r", components=(3,), origin=True) - self._moments[(1,0,0)] = lambda **kw: _r(**kw)[..., 0] - self._moments[(0,1,0)] = lambda **kw: _r(**kw)[..., 1] - self._moments[(0,0,1)] = lambda **kw: _r(**kw)[..., 2] + self._moments[(1, 0, 0)] = lambda **kw: _r(**kw)[..., 0] + self._moments[(0, 1, 0)] = lambda **kw: _r(**kw)[..., 1] + self._moments[(0, 0, 1)] = lambda **kw: _r(**kw)[..., 2] # Order 2: int1e_rr has 9 components _rr = self.make_int1e("int1e_rr", components=(9,), origin=True) - self._moments[(2,0,0)] = lambda **kw: _rr(**kw)[..., 0] - self._moments[(1,1,0)] = lambda **kw: _rr(**kw)[..., 1] - self._moments[(1,0,1)] = lambda **kw: _rr(**kw)[..., 2] - self._moments[(0,2,0)] = lambda **kw: _rr(**kw)[..., 4] - self._moments[(0,1,1)] = lambda **kw: _rr(**kw)[..., 5] - self._moments[(0,0,2)] = lambda **kw: _rr(**kw)[..., 8] + self._moments[(2, 0, 0)] = lambda **kw: _rr(**kw)[..., 0] + self._moments[(1, 1, 0)] = lambda **kw: _rr(**kw)[..., 1] + self._moments[(1, 0, 1)] = lambda **kw: _rr(**kw)[..., 2] + self._moments[(0, 2, 0)] = lambda **kw: _rr(**kw)[..., 4] + self._moments[(0, 1, 1)] = lambda **kw: _rr(**kw)[..., 5] + self._moments[(0, 0, 2)] = lambda **kw: _rr(**kw)[..., 8] # Order 3: int1e_rrr has 27 components _rrr = self.make_int1e("int1e_rrr", components=(27,), origin=True) - self._moments[(3,0,0)] = lambda **kw: _rrr(**kw)[..., 0] - self._moments[(0,3,0)] = lambda **kw: _rrr(**kw)[..., 13] - self._moments[(0,0,3)] = lambda **kw: _rrr(**kw)[..., 26] - self._moments[(2,1,0)] = lambda **kw: _rrr(**kw)[..., 3] - self._moments[(2,0,1)] = lambda **kw: _rrr(**kw)[..., 6] - self._moments[(1,2,0)] = lambda **kw: _rrr(**kw)[..., 1] - self._moments[(0,2,1)] = lambda **kw: _rrr(**kw)[..., 14] - self._moments[(1,0,2)] = lambda **kw: _rrr(**kw)[..., 2] - self._moments[(0,1,2)] = lambda **kw: _rrr(**kw)[..., 17] - self._moments[(1,1,1)] = lambda **kw: _rrr(**kw)[..., 4] + self._moments[(3, 0, 0)] = lambda **kw: _rrr(**kw)[..., 0] + self._moments[(0, 3, 0)] = lambda **kw: _rrr(**kw)[..., 13] + self._moments[(0, 0, 3)] = lambda **kw: _rrr(**kw)[..., 26] + self._moments[(2, 1, 0)] = lambda **kw: _rrr(**kw)[..., 3] + self._moments[(2, 0, 1)] = lambda **kw: _rrr(**kw)[..., 6] + self._moments[(1, 2, 0)] = lambda **kw: _rrr(**kw)[..., 1] + self._moments[(0, 2, 1)] = lambda **kw: _rrr(**kw)[..., 14] + self._moments[(1, 0, 2)] = lambda **kw: _rrr(**kw)[..., 2] + self._moments[(0, 1, 2)] = lambda **kw: _rrr(**kw)[..., 17] + self._moments[(1, 1, 1)] = lambda **kw: _rrr(**kw)[..., 4] # Set proper value for inverse sqrt of overlap integral # for cartesian basis sets @@ -1088,7 +1102,7 @@ def nuclear_attraction_integral(self, notation="physicist", transform=None): """ return self._nuc(notation=notation, transform=transform) - def overlap(self): + def overlap(self, transform=None): r""" Compute the overlap integrals. @@ -1099,20 +1113,39 @@ def overlap(self): .. math:: S_{ij} = \langle \phi_i | \phi_j \rangle + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + Returns ------- out : np.ndarray(Nbasis, Nbasis, dtype=float) Overlap integral array. """ - out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order='F') - libcint_bindings.overlap_integral_shellloop( - out, self.natm, self.atm, self.nbas, - self.bas, self.env, self._offs, self.nbfn + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.overlap_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 - def kinetic_energy(self): + def kinetic_energy(self, transform=None): r""" Compute the kinetic energy integrals. @@ -1123,20 +1156,40 @@ def kinetic_energy(self): .. math:: T_{ij} = \langle \phi_i | -\frac{1}{2}\nabla^2 | \phi_j \rangle + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + Returns ------- out : np.ndarray(Nbasis, Nbasis, dtype=float) Kinetic energy integral array. """ - out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order='F') - libcint_bindings.kinetic_integral_shellloop( - out, self.natm, self.atm, self.nbas, - self.bas, self.env, self._offs, self.nbfn + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.kinetic_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, ) + + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 - def nuclear_attraction(self): + def nuclear_attraction(self, transform=None): r""" Compute the nuclear attraction integrals. @@ -1149,21 +1202,40 @@ def nuclear_attraction(self): where :math:`Z_A` is the nuclear charge and :math:`\mathbf{R}_A` is the position of atom :math:`A`. - + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. Returns ------- out : np.ndarray(Nbasis, Nbasis, dtype=float) Nuclear attraction integral array. """ - out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order='F') - libcint_bindings.nuclear_integral_shellloop( - out, self.natm, self.atm, self.nbas, - self.bas, self.env, self._offs, self.nbfn + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.nuclear_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + - def momentum(self): + def momentum(self, origin=None, transform=None): r""" Compute the momentum integrals. @@ -1174,28 +1246,29 @@ def momentum(self): .. math:: p_{ij} = \langle \phi_i | -i\nabla | \phi_j \rangle + Parameters + ---------- + origin : np.ndarray(3, dtype=float), default=[0, 0, 0] + Origin about which to evaluate integrals. + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + Returns ------- - out : np.ndarray(Nbasis, Nbasis, dtype=float) + out : np.ndarray(Nbasis, Nbasis, 3, dtype=complex) Momentum integral array. Notes ----- - Returns the raw single-component output from ``int1e_ipovlp_sph`` - without the :math:`-i` scaling factor. The full 3-component momentum - integral (x, y, z) with proper scaling is available via - ``momentum_integral()``. - + Returns the full 3-component complex momentum integral (x, y, z) + with proper :math:`-i` scaling. Equivalent to ``momentum_integral()``. """ + if origin is None: + origin = np.zeros(3) + return self._mom(origin=origin, transform=transform) - out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order='F') - libcint_bindings.momentum_integral_shellloop( - out, self.natm, self.atm, self.nbas, - self.bas, self.env, self._offs, self.nbfn - ) - return out - - def rinv(self): + def rinv(self, inv_origin=None, transform=None): r""" Compute the :math:`1/\left|\mathbf{r} - \mathbf{R}_\text{inv}\right|` integrals. @@ -1205,21 +1278,46 @@ def rinv(self): .. math:: V_{ij} = \langle \phi_i | \frac{1}{|\mathbf{r} - \mathbf{R}_\text{inv}|} | \phi_j \rangle - - Returns - ------- - out : np.ndarray(Nbasis, Nbasis, dtype=float) - 1/r integral array. - + + Parameters + ---------- + inv_origin : np.ndarray(3, dtype=float), optional + Origin for 1/|r - R| operator. + Default is [0, 0, 0]. + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + 1/r integral array. + """ - out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order='F') - libcint_bindings.rinv_integral_shellloop( - out, self.natm, self.atm, self.nbas, - self.bas, self.env, self._offs, self.nbfn + if inv_origin is None: + inv_origin = np.zeros(3) + self.env[4:7] = inv_origin + + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.rinv_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 - def dipole(self): + + def dipole(self, transform=None): r""" Compute the dipole moment integrals. @@ -1230,6 +1328,12 @@ def dipole(self): .. math:: \mu_{ij} = \langle \phi_i | \mathbf{r} | \phi_j \rangle + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + Returns ------- out : np.ndarray(Nbasis, Nbasis, dtype=float) @@ -1242,14 +1346,28 @@ def dipole(self): available via ``moment_integral()``. """ - out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order='F') - libcint_bindings.dipole_integral_shellloop( - out, self.natm, self.atm, self.nbas, - self.bas, self.env, self._offs, self.nbfn + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.dipole_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + - def quadrupole(self): + def quadrupole(self, transform=None): r""" Compute the quadrupole moment integrals. @@ -1260,6 +1378,13 @@ def quadrupole(self): .. math:: Q_{ij} = \langle \phi_i | \mathbf{r}\mathbf{r} | \phi_j \rangle + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + Returns ------- out : np.ndarray(Nbasis, Nbasis, dtype=float) @@ -1272,14 +1397,28 @@ def quadrupole(self): available via ``moment_integral()``. """ - out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order='F') - libcint_bindings.quadrupole_integral_shellloop( - out, self.natm, self.atm, self.nbas, - self.bas, self.env, self._offs, self.nbfn + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.quadrupole_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 - def octupole(self): + + def octupole(self, transform=None): r""" Compute the octupole moment integrals. @@ -1289,6 +1428,12 @@ def octupole(self): .. math:: O_{ij} = \langle \phi_i | \mathbf{r}\mathbf{r}\mathbf{r} | \phi_j \rangle + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. Returns ------- @@ -1302,14 +1447,450 @@ def octupole(self): available via ``moment_integral()``. """ - out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order='F') - libcint_bindings.octupole_integral_shellloop( - out, self.natm, self.atm, self.nbas, - self.bas, self.env, self._offs, self.nbfn + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.octupole_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + def gradient_kinetic(self, transform=None): + r""" + Compute the gradient of kinetic energy integrals (i∇ kinetic). + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + Gradient kinetic integral array. + """ + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.ipkin_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + def gradient_nuclear(self, transform=None): + r""" + Compute the gradient of nuclear attraction integrals (i∇ nuclear). + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + Gradient nuclear integral array. + """ + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.ipnuc_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + def gradient_rinv(self, inv_origin=None, transform=None): + r""" + Compute the gradient of 1/r integrals (i∇ rinv). + + Parameters + ---------- + inv_origin : np.ndarray(3, dtype=float), optional + Origin for 1/|r - R| operator. Default is [0, 0, 0]. + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + Gradient 1/r integral array. + """ + if inv_origin is None: + inv_origin = np.zeros(3) + self.env[4:7] = inv_origin + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.iprinv_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + def ia01p(self, transform=None): + r""" + Compute the GIAO paramagnetic shielding integrals (ia01p). + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + GIAO ia01p integral array. + """ + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.ia01p_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + def ircxp(self, transform=None): + r""" + Compute the GIAO angular momentum integrals (ircxp). + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + GIAO ircxp integral array. + """ + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.ircxp_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + + def iking(self, transform=None): + r""" + Compute the GIAO kinetic energy integrals (igkin). + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + GIAO igkin integral array. + """ + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.igkin_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + def iovlpg(self, transform=None): + r""" + Compute the GIAO overlap gradient integrals (igovlp). + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + GIAO igovlp integral array. + """ + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.igovlp_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + def inucg(self, transform=None): + r""" + Compute the GIAO nuclear attraction integrals (ignuc). + + Parameters + ---------- + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, dtype=float) + GIAO ignuc integral array. + """ + out = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.ignuc_integral_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, ) + # Apply permutation + out = out[self._permutations, :][:, self._permutations] + # 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 + + def point_charge(self, point_coords, point_charges): + r""" + Compute the point charge integrals. + + The point charge integral represents the electrostatic potential due to + a set of point charges at given coordinates. For each pair of basis + functions :math:`\phi_i` and :math:`\phi_j`, it is defined as: + + .. math:: + V_{ij}^{(n)} = -q_n \langle \phi_i | \frac{1}{|\mathbf{r} - \mathbf{R}_n|} | \phi_j \rangle + + Parameters + ---------- + point_coords : np.ndarray(N, 3, dtype=float) + Coordinates of point charges. + point_charges : np.ndarray(N, dtype=float) + Charges of point charges. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, N, dtype=float) + Point charge integral array. + + """ + out = np.zeros((self.nbfn, self.nbfn, len(point_charges)), dtype=c_double, order="F") + for icharge, (coord, charge) in enumerate(zip(point_coords, point_charges)): + # Set inv_origin in env for this charge + self.env[4:7] = coord + val = np.zeros((self.nbfn, self.nbfn), dtype=c_double, order="F") + libcint_bindings.rinv_integral_array( + val, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + val *= -charge + out[:, :, icharge] = val + return out + + def moment(self, orders, origin=None): + r""" + Compute the moment integrals. + + The moment integral represents the expectation value of the position + operator raised to a given order between basis functions :math:`\phi_i` + and :math:`\phi_j`. + + Parameters + ---------- + orders : np.ndarray(N, 3, dtype=int) + Moment orders :math:`\left[x, y, z\right]` to evaluate. + origin : np.ndarray(3, dtype=float), default=[0, 0, 0] + Origin about which to evaluate integrals. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, N, dtype=float) + Moment integral array. + + Notes + ----- + Uses the C shell-loop bindings for dipole, quadrupole, and octupole + integrals internally. Supports up to 3rd order moments. + + """ + if origin is None: + origin = np.zeros(3) + out = np.zeros((self.nbfn, self.nbfn, len(orders)), dtype=np.float64) + for i, order in enumerate(orders): + self.env[1:4] = origin + if sum(order) == 0: + out[:, :, i] = self.overlap() + else: + out[:, :, i] = self._moments[tuple(order)](origin=origin) return out - + def electron_repulsion(self, notation="physicist", transform=None): + r""" + Compute the electron repulsion integrals. + + The two-electron repulsion integral between basis functions + :math:`\phi_i`, :math:`\phi_j`, :math:`\phi_k`, and :math:`\phi_l` + is defined as: + + .. math:: + g_{ijkl} = \langle \phi_i \phi_j | \frac{1}{r_{12}} | \phi_k \phi_l \rangle + + Parameters + ---------- + notation : ("physicist" | "chemist"), default="physicist" + Axis order convention. + transform : np.ndarray(K, K_cont), optional + Transformation matrix from AO to MO basis. + Default is no transformation. + + Returns + ------- + out : np.ndarray(Nbasis, Nbasis, Nbasis, Nbasis, dtype=float) + Electron repulsion integral array. + + """ + if notation not in ("physicist", "chemist"): + raise ValueError("``notation`` must be one of 'physicist' or 'chemist'") + + out = np.zeros((self.nbfn, self.nbfn, self.nbfn, self.nbfn), dtype=c_double) + libcint_bindings.eri_array( + out, + self.natm, + self.atm, + self.nbas, + self.bas, + self.env, + self._offs, + self.nbfn, + ) + + # Apply permutation + out = out[self._permutations] + out = out[:, self._permutations] + out = out[:, :, self._permutations] + out = out[:, :, :, self._permutations] + # Apply notation + if notation == "chemist": + out = out.transpose(0, 2, 1, 3) + # 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 def electron_repulsion_integral(self, notation="physicist", transform=None): r""" diff --git a/gbasis/integrals/src/libcint_wrap.c b/gbasis/integrals/src/libcint_wrap.c index 9d639b48..e7373b1f 100644 --- a/gbasis/integrals/src/libcint_wrap.c +++ b/gbasis/integrals/src/libcint_wrap.c @@ -37,7 +37,9 @@ #include #include #include - +/* CINTOpt forward declaration */ +typedef struct CINTOpt CINTOpt; +extern void CINTdel_optimizer(CINTOpt **); /* Forward declarations for libcint spherical integral functions. * Signature: (out, dims, shls, atm, natm, bas, nbas, env, opt, cache) @@ -54,60 +56,122 @@ */ /* Forward declarations — 1-electron integrals */ -extern int int1e_ovlp_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); -extern int int1e_kin_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); -extern int int1e_nuc_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); -extern int int1e_ipovlp_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); -extern int int1e_cg_irxp_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); -extern int int1e_rinv_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); -extern int int1e_r_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); -extern int int1e_rr_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); -extern int int1e_rrr_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); +extern int int1e_ovlp_sph(double *out, int *dims, int *shls, int *atm, int natm, + int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_kin_sph(double *out, int *dims, int *shls, int *atm, int natm, + int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_nuc_sph(double *out, int *dims, int *shls, int *atm, int natm, + int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_ipovlp_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, + void *opt, double *cache); +extern int int1e_cg_irxp_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, + void *opt, double *cache); +extern int int1e_rinv_sph(double *out, int *dims, int *shls, int *atm, int natm, + int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_r_sph(double *out, int *dims, int *shls, int *atm, int natm, + int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_rr_sph(double *out, int *dims, int *shls, int *atm, int natm, + int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_rrr_sph(double *out, int *dims, int *shls, int *atm, int natm, + int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_ipkin_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_ipnuc_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_iprinv_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, + void *opt, double *cache); +extern int int1e_ia01p_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_igkin_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, void *opt, + double *cache); +extern int int1e_igovlp_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, + void *opt, double *cache); +extern int int1e_ignuc_sph(double *out, int *dims, int *shls, int *atm, + int natm, int *bas, int nbas, double *env, void *opt, + double *cache); + +/* Optimizer forward declarations */ +extern void int1e_ovlp_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_kin_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_nuc_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_ipovlp_optimizer(CINTOpt **, int *, int, int *, int, + double *); +extern void int1e_cg_irxp_optimizer(CINTOpt **, int *, int, int *, int, + double *); +extern void int1e_rinv_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_r_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_rr_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_rrr_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void cint2e_sph_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void CINTall_1e_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_ipkin_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_ipnuc_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_iprinv_optimizer(CINTOpt **, int *, int, int *, int, + double *); +extern void int1e_ia01p_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_igkin_optimizer(CINTOpt **, int *, int, int *, int, double *); +extern void int1e_igovlp_optimizer(CINTOpt **, int *, int, int *, int, + double *); +extern void int1e_ignuc_optimizer(CINTOpt **, int *, int, int *, int, double *); /* Forward declaration — 2-electron integral */ -extern int int2e_sph(double *out, int *dims, int *shls, int *atm, int natm, int *bas, int nbas, double *env, void *opt, double *cache); +extern int int2e_sph(double *out, int *dims, int *shls, int *atm, int natm, + int *bas, int nbas, double *env, void *opt, double *cache); +/* Optimizer macro using token pasting */ +#define MAKE_OPTIMIZER(func) c##func##_optimizer /* - * DEFINE_INTEGRAL_INT1e(func_name, libcint_func) + * DEFINE_INT1E_ARRAY_FN(func_name, libcint_func) * Generates a Python/C API wrapper for a 1-electron libcint integral. * Accepts NumPy arrays directly — uses PyArray_GETPTR1 (NumPy C-API). */ /* Macro for 1-electron wrappers */ -#define DEFINE_INTEGRAL_INT1e(func_name, libcint_func) \ -static PyObject * \ -func_name(PyObject *self, PyObject *args) \ -{ \ - PyArrayObject *out_arr, *dims_arr, *shls_arr, *atm_arr, *bas_arr, *env_arr; \ +#define DEFINE_INT1E_ARRAY_FN(func_name, libcint_func) \ + static PyObject *func_name(PyObject *self, PyObject *args) { \ + PyArrayObject *out_arr, *dims_arr, *shls_arr, *atm_arr, *bas_arr, \ + *env_arr; \ int natm, nbas; \ - if (!PyArg_ParseTuple(args, "O!O!O!O!iO!iO!", \ - &PyArray_Type, &out_arr, \ - &PyArray_Type, &dims_arr, \ - &PyArray_Type, &shls_arr, \ - &PyArray_Type, &atm_arr, &natm, \ - &PyArray_Type, &bas_arr, &nbas, \ - &PyArray_Type, &env_arr)) \ - return NULL; \ - double *out = (double *)PyArray_GETPTR1(out_arr, 0); \ - int *dims = (int *) PyArray_GETPTR1(dims_arr, 0); \ - int *shls = (int *) PyArray_GETPTR1(shls_arr, 0); \ - int *atm = (int *) PyArray_GETPTR2(atm_arr, 0, 0); \ - int *bas = (int *) PyArray_GETPTR2(bas_arr, 0, 0); \ - double *env = (double *)PyArray_GETPTR1(env_arr, 0); \ - int result = libcint_func(out, dims, shls, atm, natm, \ - bas, nbas, env, NULL, NULL); \ + if (!PyArg_ParseTuple(args, "O!O!O!O!iO!iO!", &PyArray_Type, &out_arr, \ + &PyArray_Type, &dims_arr, &PyArray_Type, &shls_arr, \ + &PyArray_Type, &atm_arr, &natm, &PyArray_Type, \ + &bas_arr, &nbas, &PyArray_Type, &env_arr)) \ + return NULL; \ + double *out = (double *)PyArray_GETPTR1(out_arr, 0); \ + int *dims = (int *)PyArray_GETPTR1(dims_arr, 0); \ + int *shls = (int *)PyArray_GETPTR1(shls_arr, 0); \ + int *atm = (int *)PyArray_GETPTR2(atm_arr, 0, 0); \ + int *bas = (int *)PyArray_GETPTR2(bas_arr, 0, 0); \ + double *env = (double *)PyArray_GETPTR1(env_arr, 0); \ + int result = \ + libcint_func(out, dims, shls, atm, natm, bas, nbas, env, NULL, NULL); \ return PyLong_FromLong(result); \ -} + } -DEFINE_INTEGRAL_INT1e(overlap_sph, int1e_ovlp_sph) -DEFINE_INTEGRAL_INT1e(kinetic_sph, int1e_kin_sph) -DEFINE_INTEGRAL_INT1e(nuclear_sph, int1e_nuc_sph) -DEFINE_INTEGRAL_INT1e(momentum_sph, int1e_ipovlp_sph) -DEFINE_INTEGRAL_INT1e(angular_momentum_sph, int1e_cg_irxp_sph) -DEFINE_INTEGRAL_INT1e(rinv_sph, int1e_rinv_sph) -DEFINE_INTEGRAL_INT1e(dipole_sph, int1e_r_sph) -DEFINE_INTEGRAL_INT1e(quadrupole_sph, int1e_rr_sph) -DEFINE_INTEGRAL_INT1e(octupole_sph, int1e_rrr_sph) +DEFINE_INT1E_ARRAY_FN(overlap_sph, int1e_ovlp_sph) +DEFINE_INT1E_ARRAY_FN(kinetic_sph, int1e_kin_sph) +DEFINE_INT1E_ARRAY_FN(nuclear_sph, int1e_nuc_sph) +DEFINE_INT1E_ARRAY_FN(momentum_sph, int1e_ipovlp_sph) +DEFINE_INT1E_ARRAY_FN(angular_momentum_sph, int1e_cg_irxp_sph) +DEFINE_INT1E_ARRAY_FN(rinv_sph, int1e_rinv_sph) +DEFINE_INT1E_ARRAY_FN(dipole_sph, int1e_r_sph) +DEFINE_INT1E_ARRAY_FN(quadrupole_sph, int1e_rr_sph) +DEFINE_INT1E_ARRAY_FN(octupole_sph, int1e_rrr_sph) /* 2-electron wrapper */ @@ -117,131 +181,306 @@ DEFINE_INTEGRAL_INT1e(octupole_sph, int1e_rrr_sph) * Same signature as 1-electron but uses int2e_sph (4-center integral). */ -static PyObject * -electron_repulsion_sph(PyObject *self, PyObject *args) -{ - PyArrayObject *out_arr, *dims_arr, *shls_arr, *atm_arr, *bas_arr, *env_arr; - int natm, nbas; - if (!PyArg_ParseTuple(args, "O!O!O!O!iO!iO!", - &PyArray_Type, &out_arr, - &PyArray_Type, &dims_arr, - &PyArray_Type, &shls_arr, - &PyArray_Type, &atm_arr, &natm, - &PyArray_Type, &bas_arr, &nbas, - &PyArray_Type, &env_arr)) - return NULL; - double *out = (double *)PyArray_GETPTR1(out_arr, 0); - int *dims = (int *) PyArray_GETPTR1(dims_arr, 0); - int *shls = (int *) PyArray_GETPTR1(shls_arr, 0); - int *atm = (int *) PyArray_GETPTR2(atm_arr, 0, 0); - int *bas = (int *) PyArray_GETPTR2(bas_arr, 0, 0); - double *env = (double *)PyArray_GETPTR1(env_arr, 0); - int result = int2e_sph(out, dims, shls, atm, natm, - bas, nbas, env, NULL, NULL); - return PyLong_FromLong(result); +static PyObject *electron_repulsion_sph(PyObject *self, PyObject *args) { + PyArrayObject *out_arr, *dims_arr, *shls_arr, *atm_arr, *bas_arr, *env_arr; + int natm, nbas; + if (!PyArg_ParseTuple(args, "O!O!O!O!iO!iO!", &PyArray_Type, &out_arr, + &PyArray_Type, &dims_arr, &PyArray_Type, &shls_arr, + &PyArray_Type, &atm_arr, &natm, &PyArray_Type, &bas_arr, + &nbas, &PyArray_Type, &env_arr)) + return NULL; + double *out = (double *)PyArray_GETPTR1(out_arr, 0); + int *dims = (int *)PyArray_GETPTR1(dims_arr, 0); + int *shls = (int *)PyArray_GETPTR1(shls_arr, 0); + int *atm = (int *)PyArray_GETPTR2(atm_arr, 0, 0); + int *bas = (int *)PyArray_GETPTR2(bas_arr, 0, 0); + double *env = (double *)PyArray_GETPTR1(env_arr, 0); + int result = + int2e_sph(out, dims, shls, atm, natm, bas, nbas, env, NULL, NULL); + return PyLong_FromLong(result); } +/* Macro for integrals WITH optimizer support */ +#define DEFINE_INT1E_LOOP_FN_OPT(func_name, libcint_func, opt_func) \ + static PyObject *func_name(PyObject *self, PyObject *args) { \ + PyArrayObject *out_arr, *atm_arr, *bas_arr, *env_arr, *offs_arr; \ + int natm, nbas, nbfn; \ + if (!PyArg_ParseTuple(args, "O!iO!iO!O!O!i", &PyArray_Type, &out_arr, \ + &natm, &PyArray_Type, &atm_arr, &nbas, \ + &PyArray_Type, &bas_arr, &PyArray_Type, &env_arr, \ + &PyArray_Type, &offs_arr, &nbfn)) \ + return NULL; \ + double *out = (double *)PyArray_DATA(out_arr); \ + int *atm = (int *)PyArray_GETPTR2(atm_arr, 0, 0); \ + int *bas = (int *)PyArray_GETPTR2(bas_arr, 0, 0); \ + double *env = (double *)PyArray_GETPTR1(env_arr, 0); \ + int *offs = (int *)PyArray_GETPTR1(offs_arr, 0); \ + int max_off = 0; \ + for (int i = 0; i < nbas; i++) { \ + if (offs[i] > max_off) \ + max_off = offs[i]; \ + } \ + size_t buf_size = (size_t)max_off * max_off; \ + double *buf = calloc(buf_size, sizeof(double)); \ + if (!buf) { \ + PyErr_NoMemory(); \ + return NULL; \ + } \ + CINTOpt *opt = NULL; \ + opt_func##_optimizer(&opt, atm, natm, bas, nbas, env); \ + int shls[2]; \ + int ipos = 0; \ + for (int ishl = 0; ishl < nbas; ishl++) { \ + shls[0] = ishl; \ + int p_off = offs[ishl]; \ + int jpos = 0; \ + for (int jshl = 0; jshl <= ishl; jshl++) { \ + shls[1] = jshl; \ + int q_off = offs[jshl]; \ + libcint_func(buf, NULL, shls, atm, natm, bas, nbas, env, opt, NULL); \ + for (int p = 0; p < p_off; p++) { \ + for (int q = 0; q < q_off; q++) { \ + double val = buf[p + q * p_off]; \ + out[(ipos + p) * nbfn + (jpos + q)] = val; \ + out[(jpos + q) * nbfn + (ipos + p)] = val; \ + } \ + } \ + memset(buf, 0, buf_size * sizeof(double)); \ + jpos += q_off; \ + } \ + ipos += p_off; \ + } \ + CINTdel_optimizer(&opt); \ + free(buf); \ + Py_RETURN_NONE; \ + } + /* - * DEFINE_SHELLLOOP_INT1e(func_name, libcint_func) + * DEFINE_INT1E_LOOP_FN(func_name, libcint_func) * Generates a C shell-loop wrapper for a 1-electron libcint integral that * returns the FULL integral array over all shells (not just one shell pair). * Loops over shells I, J; calls libcint_func per shell pair; fills the * symmetric output matrix using PyArray_GETPTR (NumPy C-API). */ -#define DEFINE_SHELLLOOP_INT1e(func_name, libcint_func) \ -static PyObject * \ -func_name(PyObject *self, PyObject *args) \ -{ \ - PyArrayObject *out_arr, *atm_arr, *bas_arr, *env_arr, *offs_arr; \ - int natm, nbas, nbfn; \ - if (!PyArg_ParseTuple(args, "O!iO!iO!O!O!i", \ - &PyArray_Type, &out_arr, \ - &natm, \ - &PyArray_Type, &atm_arr, \ - &nbas, \ - &PyArray_Type, &bas_arr, \ - &PyArray_Type, &env_arr, \ - &PyArray_Type, &offs_arr, \ - &nbfn)) \ - return NULL; \ - double *out = (double *)PyArray_DATA(out_arr); \ - int *atm = (int *) PyArray_GETPTR2(atm_arr, 0, 0); \ - int *bas = (int *) PyArray_GETPTR2(bas_arr, 0, 0); \ - double *env = (double *)PyArray_GETPTR1(env_arr, 0); \ - int *offs = (int *) PyArray_GETPTR1(offs_arr, 0); \ - int shls[2]; \ - double buf[10000] = {0}; \ - int ipos = 0; \ - for (int ishl = 0; ishl < nbas; ishl++) { \ - shls[0] = ishl; \ - int p_off = offs[ishl]; \ - int jpos = 0; \ - for (int jshl = 0; jshl <= ishl; jshl++) { \ - shls[1] = jshl; \ - int q_off = offs[jshl]; \ - libcint_func(buf, NULL, shls, atm, natm, bas, nbas, env, NULL, NULL); \ - for (int p = 0; p < p_off; p++) { \ - for (int q = 0; q < q_off; q++) { \ - double val = buf[p + q * p_off]; \ - out[(ipos+p) * nbfn + (jpos+q)] = val; \ - out[(jpos+q) * nbfn + (ipos+p)] = val; \ - } \ - } \ - memset(buf, 0, sizeof(buf)); \ - jpos += q_off; \ - } \ - ipos += p_off; \ - } \ - Py_RETURN_NONE; \ -} +#define DEFINE_INT1E_LOOP_FN(func_name, libcint_func, opt_func) \ + static PyObject *func_name(PyObject *self, PyObject *args) { \ + PyArrayObject *out_arr, *atm_arr, *bas_arr, *env_arr, *offs_arr; \ + int natm, nbas, nbfn; \ + if (!PyArg_ParseTuple(args, "O!iO!iO!O!O!i", &PyArray_Type, &out_arr, \ + &natm, &PyArray_Type, &atm_arr, &nbas, \ + &PyArray_Type, &bas_arr, &PyArray_Type, &env_arr, \ + &PyArray_Type, &offs_arr, &nbfn)) \ + return NULL; \ + double *out = (double *)PyArray_DATA(out_arr); \ + int *atm = (int *)PyArray_GETPTR2(atm_arr, 0, 0); \ + int *bas = (int *)PyArray_GETPTR2(bas_arr, 0, 0); \ + double *env = (double *)PyArray_GETPTR1(env_arr, 0); \ + int *offs = (int *)PyArray_GETPTR1(offs_arr, 0); \ + int shls[2]; \ + int max_off = 0; \ + for (int i = 0; i < nbas; i++) { \ + if (offs[i] > max_off) \ + max_off = offs[i]; \ + } \ + size_t buf_size = (size_t)max_off * max_off * 27; \ + double *buf = calloc(buf_size, sizeof(double)); \ + if (!buf) { \ + PyErr_NoMemory(); \ + return NULL; \ + } \ + CINTOpt *opt = NULL; \ + opt_func##_optimizer(&opt, atm, natm, bas, nbas, env); \ + int ipos = 0; \ + for (int ishl = 0; ishl < nbas; ishl++) { \ + shls[0] = ishl; \ + int p_off = offs[ishl]; \ + int jpos = 0; \ + for (int jshl = 0; jshl <= ishl; jshl++) { \ + shls[1] = jshl; \ + int q_off = offs[jshl]; \ + libcint_func(buf, NULL, shls, atm, natm, bas, nbas, env, opt, NULL); \ + for (int p = 0; p < p_off; p++) { \ + for (int q = 0; q < q_off; q++) { \ + double val = buf[p + q * p_off]; \ + out[(ipos + p) * nbfn + (jpos + q)] = val; \ + out[(jpos + q) * nbfn + (ipos + p)] = val; \ + } \ + } \ + memset(buf, 0, buf_size * sizeof(double)); \ + jpos += q_off; \ + } \ + ipos += p_off; \ + } \ + CINTdel_optimizer(&opt); \ + free(buf); \ + Py_RETURN_NONE; \ + } /* Generate shell-loop wrappers for all 1-electron integrals using the macro */ -DEFINE_SHELLLOOP_INT1e(overlap_integral_shellloop, int1e_ovlp_sph) -DEFINE_SHELLLOOP_INT1e(kinetic_integral_shellloop, int1e_kin_sph) -DEFINE_SHELLLOOP_INT1e(nuclear_integral_shellloop, int1e_nuc_sph) -DEFINE_SHELLLOOP_INT1e(momentum_integral_shellloop, int1e_ipovlp_sph) -DEFINE_SHELLLOOP_INT1e(rinv_integral_shellloop, int1e_rinv_sph) -DEFINE_SHELLLOOP_INT1e(dipole_integral_shellloop, int1e_r_sph) -DEFINE_SHELLLOOP_INT1e(quadrupole_integral_shellloop, int1e_rr_sph) -DEFINE_SHELLLOOP_INT1e(octupole_integral_shellloop, int1e_rrr_sph) +DEFINE_INT1E_LOOP_FN_OPT(overlap_integral_array, int1e_ovlp_sph, int1e_ovlp) +DEFINE_INT1E_LOOP_FN_OPT(kinetic_integral_array, int1e_kin_sph, int1e_kin) +DEFINE_INT1E_LOOP_FN_OPT(nuclear_integral_array, int1e_nuc_sph, int1e_nuc) +DEFINE_INT1E_LOOP_FN(momentum_integral_array, int1e_ipovlp_sph, int1e_ipovlp) +DEFINE_INT1E_LOOP_FN(rinv_integral_array, int1e_rinv_sph, int1e_rinv) +DEFINE_INT1E_LOOP_FN(dipole_integral_array, int1e_r_sph, int1e_r) +DEFINE_INT1E_LOOP_FN(quadrupole_integral_array, int1e_rr_sph, int1e_rr) +DEFINE_INT1E_LOOP_FN(octupole_integral_array, int1e_rrr_sph, int1e_rrr) +DEFINE_INT1E_LOOP_FN(ipkin_integral_array, int1e_ipkin_sph, int1e_ipkin) +DEFINE_INT1E_LOOP_FN(ipnuc_integral_array, int1e_ipnuc_sph, int1e_ipnuc) +DEFINE_INT1E_LOOP_FN(iprinv_integral_array, int1e_iprinv_sph, int1e_iprinv) +DEFINE_INT1E_LOOP_FN(ia01p_integral_array, int1e_ia01p_sph, int1e_ia01p) +DEFINE_INT1E_LOOP_FN(ircxp_integral_array, int1e_cg_irxp_sph, int1e_cg_irxp) +DEFINE_INT1E_LOOP_FN(igkin_integral_array, int1e_igkin_sph, int1e_igkin) +DEFINE_INT1E_LOOP_FN(igovlp_integral_array, int1e_igovlp_sph, int1e_igovlp) +DEFINE_INT1E_LOOP_FN(ignuc_integral_array, int1e_ignuc_sph, int1e_ignuc) + +/* eri_array — array-based loop in C for 2-electron ERI (4 shells: I,J,K,L) */ +static PyObject *eri_array(PyObject *self, PyObject *args) { + PyArrayObject *out_arr, *atm_arr, *bas_arr, *env_arr, *offs_arr; + int natm, nbas, nbfn; -/* eri_shellloop — shell-by-shell loop in C for 2-electron ERI (4 shells: I,J,K,L) */ + if (!PyArg_ParseTuple(args, "O!iO!iO!O!O!i", &PyArray_Type, &out_arr, &natm, + &PyArray_Type, &atm_arr, &nbas, &PyArray_Type, &bas_arr, + &PyArray_Type, &env_arr, &PyArray_Type, &offs_arr, + &nbfn)) + return NULL; + double *out = (double *)PyArray_DATA(out_arr); + int *atm = (int *)PyArray_GETPTR2(atm_arr, 0, 0); + int *bas = (int *)PyArray_GETPTR2(bas_arr, 0, 0); + double *env = (double *)PyArray_GETPTR1(env_arr, 0); + int *offs = (int *)PyArray_GETPTR1(offs_arr, 0); + int shls[4]; + int max_off = 0; + for (int i = 0; i < nbas; i++) { + if (offs[i] > max_off) + max_off = offs[i]; + } + size_t buf_size = (size_t)max_off * max_off * max_off * max_off; + double *buf = calloc(buf_size, sizeof(double)); + if (!buf) { + PyErr_NoMemory(); + return NULL; + } + CINTOpt *opt = NULL; + cint2e_sph_optimizer(&opt, atm, natm, bas, nbas, env); + int ipos = 0; + for (int ishl = 0; ishl < nbas; ishl++) { + shls[0] = ishl; + int p_off = offs[ishl]; + int jpos = 0; + for (int jshl = 0; jshl <= ishl; jshl++) { + int ij = ((ishl + 1) * ishl) / 2 + jshl; + shls[1] = jshl; + int q_off = offs[jshl]; + int kpos = 0; + for (int kshl = 0; kshl < nbas; kshl++) { + shls[2] = kshl; + int r_off = offs[kshl]; + int lpos = 0; + for (int lshl = 0; lshl <= kshl; lshl++) { + int kl = ((kshl + 1) * kshl) / 2 + lshl; + if (ij < kl) { + lpos += offs[lshl]; + continue; + } + shls[3] = lshl; + int s_off = offs[lshl]; + int2e_sph(buf, NULL, shls, atm, natm, bas, nbas, env, opt, NULL); + for (int p = 0; p < p_off; p++) { + for (int q = 0; q < q_off; q++) { + for (int r = 0; r < r_off; r++) { + for (int s = 0; s < s_off; s++) { + double val = buf[p + p_off * (q + q_off * (r + r_off * s))]; + int i = ipos + p, j = jpos + q, k = kpos + r, l = lpos + s; + out[i * nbfn * nbfn * nbfn + k * nbfn * nbfn + j * nbfn + l] = + val; + out[i * nbfn * nbfn * nbfn + l * nbfn * nbfn + j * nbfn + k] = + val; + out[j * nbfn * nbfn * nbfn + k * nbfn * nbfn + i * nbfn + l] = + val; + out[j * nbfn * nbfn * nbfn + l * nbfn * nbfn + i * nbfn + k] = + val; + out[k * nbfn * nbfn * nbfn + i * nbfn * nbfn + l * nbfn + j] = + val; + out[k * nbfn * nbfn * nbfn + j * nbfn * nbfn + l * nbfn + i] = + val; + out[l * nbfn * nbfn * nbfn + i * nbfn * nbfn + k * nbfn + j] = + val; + out[l * nbfn * nbfn * nbfn + j * nbfn * nbfn + k * nbfn + i] = + val; + } + } + } + } + memset(buf, 0, buf_size * sizeof(double)); + lpos += s_off; + } + kpos += r_off; + } + jpos += q_off; + } + ipos += p_off; + } + CINTdel_optimizer(&opt); + free(buf); + Py_RETURN_NONE; +} static PyMethodDef LibcintMethods[] = { - {"overlap_sph", overlap_sph, METH_VARARGS, "Overlap integral"}, - {"kinetic_sph", kinetic_sph, METH_VARARGS, "Kinetic energy integral"}, - {"nuclear_sph", nuclear_sph, METH_VARARGS, "Nuclear attraction integral"}, - {"momentum_sph", momentum_sph, METH_VARARGS, "Momentum integral"}, - {"angular_momentum_sph", angular_momentum_sph, METH_VARARGS, "Angular momentum integral"}, - {"rinv_sph", rinv_sph, METH_VARARGS, "1/r integral"}, - {"dipole_sph", dipole_sph, METH_VARARGS, "Dipole moment integral"}, - {"quadrupole_sph", quadrupole_sph, METH_VARARGS, "Quadrupole moment integral"}, - {"octupole_sph", octupole_sph, METH_VARARGS, "Octupole moment integral"}, - {"electron_repulsion_sph", electron_repulsion_sph, METH_VARARGS, "Electron repulsion integral"}, - {"overlap_integral_shellloop", overlap_integral_shellloop, METH_VARARGS, "Overlap integral shell loop in C"}, - {"kinetic_integral_shellloop", kinetic_integral_shellloop, METH_VARARGS, "Kinetic integral shell loop in C"}, - {"nuclear_integral_shellloop", nuclear_integral_shellloop, METH_VARARGS, "Nuclear attraction integral shell loop in C"}, - {"momentum_integral_shellloop", momentum_integral_shellloop, METH_VARARGS, "Momentum integral shell loop in C"}, - {"rinv_integral_shellloop", rinv_integral_shellloop, METH_VARARGS, "1/r integral shell loop in C"}, - {"dipole_integral_shellloop", dipole_integral_shellloop, METH_VARARGS, "Dipole integral shell loop in C"}, - {"quadrupole_integral_shellloop", quadrupole_integral_shellloop, METH_VARARGS, "Quadrupole integral shell loop in C"}, - {"octupole_integral_shellloop", octupole_integral_shellloop, METH_VARARGS, "Octupole integral shell loop in C"}, - {NULL, NULL, 0, NULL} -}; + {"overlap_sph", overlap_sph, METH_VARARGS, "Overlap integral"}, + {"kinetic_sph", kinetic_sph, METH_VARARGS, "Kinetic energy integral"}, + {"nuclear_sph", nuclear_sph, METH_VARARGS, "Nuclear attraction integral"}, + {"momentum_sph", momentum_sph, METH_VARARGS, "Momentum integral"}, + {"angular_momentum_sph", angular_momentum_sph, METH_VARARGS, + "Angular momentum integral"}, + {"rinv_sph", rinv_sph, METH_VARARGS, "1/r integral"}, + {"dipole_sph", dipole_sph, METH_VARARGS, "Dipole moment integral"}, + {"quadrupole_sph", quadrupole_sph, METH_VARARGS, + "Quadrupole moment integral"}, + {"octupole_sph", octupole_sph, METH_VARARGS, "Octupole moment integral"}, + {"electron_repulsion_sph", electron_repulsion_sph, METH_VARARGS, + "Electron repulsion integral"}, + {"overlap_integral_array", overlap_integral_array, METH_VARARGS, + "Overlap integral array in C"}, + {"kinetic_integral_array", kinetic_integral_array, METH_VARARGS, + "Kinetic integral array in C"}, + {"nuclear_integral_array", nuclear_integral_array, METH_VARARGS, + "Nuclear attraction integral array in C"}, + {"momentum_integral_array", momentum_integral_array, METH_VARARGS, + "Momentum integral array in C"}, + {"rinv_integral_array", rinv_integral_array, METH_VARARGS, + "1/r integral array in C"}, + {"dipole_integral_array", dipole_integral_array, METH_VARARGS, + "Dipole integral array in C"}, + {"quadrupole_integral_array", quadrupole_integral_array, METH_VARARGS, + "Quadrupole integral array in C"}, + {"octupole_integral_array", octupole_integral_array, METH_VARARGS, + "Octupole integral array in C"}, + {"ipkin_integral_array", ipkin_integral_array, METH_VARARGS, + "Momentum kinetic integral array in C"}, + {"ipnuc_integral_array", ipnuc_integral_array, METH_VARARGS, + "Momentum nuclear attraction integral array in C"}, + {"iprinv_integral_array", iprinv_integral_array, METH_VARARGS, + "Momentum 1/r integral array in C"}, + {"ia01p_integral_array", ia01p_integral_array, METH_VARARGS, + "GIAO ia01p integral array in C"}, + {"ircxp_integral_array", ircxp_integral_array, METH_VARARGS, + "GIAO ircxp integral array in C"}, + {"igkin_integral_array", igkin_integral_array, METH_VARARGS, + "GIAO igkin integral array in C"}, + {"igovlp_integral_array", igovlp_integral_array, METH_VARARGS, + "GIAO igovlp integral array in C"}, + {"ignuc_integral_array", ignuc_integral_array, METH_VARARGS, + "GIAO ignuc integral array in C"}, + {"eri_array", eri_array, METH_VARARGS, "ERI 2-electron array in C"}, + {NULL, NULL, 0, NULL}}; static struct PyModuleDef libcintmodule = { - PyModuleDef_HEAD_INIT, - "libcint_bindings", - NULL, - -1, - LibcintMethods -}; - -PyMODINIT_FUNC -PyInit_libcint_bindings(void) -{ - import_array(); - return PyModule_Create(&libcintmodule); + PyModuleDef_HEAD_INIT, "libcint_bindings", NULL, -1, LibcintMethods}; + +PyMODINIT_FUNC PyInit_libcint_bindings(void) { + import_array(); + return PyModule_Create(&libcintmodule); } diff --git a/tests/test_libcint.py b/tests/test_libcint.py index 1d1deddb..14f42aca 100644 --- a/tests/test_libcint.py +++ b/tests/test_libcint.py @@ -68,7 +68,7 @@ @pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows") @pytest.mark.skipif( - len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.so*"))) == 0, + len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.*"))) == 0, reason="The libcint shared library object was not found", ) @pytest.mark.parametrize("integral", TEST_INTEGRALS) @@ -198,7 +198,7 @@ def test_integral(basis, atsyms, atcoords, coord_type, integral): ] @pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows") @pytest.mark.skipif( - len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.so*"))) == 0, + len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.*"))) == 0, reason="The libcint shared library object was not found", ) @pytest.mark.parametrize("fname, elements, coord_type", TEST_SYSTEMS_IODATA) @@ -343,3 +343,456 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform): raise ValueError("Invalid integral name '{integral}' passed") npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + + + + +# ───────────────────────────────────────────────────────────────────────── +# New test list for C shell-loop bindings introduced in PR-5 and PR-6 +# ───────────────────────────────────────────────────────────────────────── + +TEST_C_SHELLLOOP_INTEGRALS = [ + pytest.param("overlap", id="C-Overlap"), + pytest.param("kinetic_energy", id="C-KineticEnergy"), + pytest.param("nuclear_attraction", id="C-NuclearAttraction"), + pytest.param("rinv", id="C-Rinv"), + pytest.param("momentum", id="C-Momentum"), + pytest.param("dipole", id="C-Dipole"), + pytest.param("quadrupole", id="C-Quadrupole"), + pytest.param("octupole", id="C-Octupole"), + pytest.param("point_charge", id="C-PointCharge"), + pytest.param("moment", id="C-Moment"), + pytest.param("electron_repulsion", id="C-ElectronRepulsion"), +] + + +@pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows") +@pytest.mark.skipif( + len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.*"))) == 0, + reason="The libcint shared library object was not found", +) +@pytest.mark.parametrize("integral", TEST_C_SHELLLOOP_INTEGRALS) +@pytest.mark.parametrize("atsyms, atcoords", TEST_SYSTEMS) +@pytest.mark.parametrize("basis", TEST_BASIS_SETS) +def test_c_shellloop_integral(basis, atsyms, atcoords, integral): + r""" + Test the C shell-loop bindings (PR-5: 1-electron, PR-6: ERI) added to + ``gbasis.integrals.libcint.CBasis`` against the existing GBasis Python + integral implementations. + + These are the ``.overlap()``, ``.kinetic_energy()``, + ``.nuclear_attraction()``, ``.rinv()``, ``.dipole()``, ``.quadrupole()``, + ``.octupole()``, and ``.electron_repulsion()`` methods, which loop over + shells directly in C (as opposed to the ``*_integral()`` methods, which + loop over shells in Python and only call into C per shell pair). + + """ + from gbasis.integrals.libcint import ELEMENTS, CBasis + + atol, rtol = 1e-6, 1e-6 + + atcoords = atcoords / 0.5291772083 + + atnums = np.asarray([ELEMENTS.index(i) for i in atsyms], dtype=float) + + basis_dict = parse_nwchem(find_datafile(basis)) + + # C shell-loop bindings are implemented for spherical only + py_basis = make_contractions(basis_dict, atsyms, atcoords, coord_types="spherical") + + lc_basis = CBasis(py_basis, atsyms, atcoords, coord_type="spherical") + + if integral == "overlap": + py_int = overlap_integral(py_basis, screen_basis=False) + lc_int = lc_basis.overlap() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.overlap(transform=transform) + npt.assert_array_equal(lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int_t, py_int, atol=atol, rtol=rtol) + + elif integral == "kinetic_energy": + py_int = kinetic_energy_integral(py_basis, screen_basis=False) + lc_int = lc_basis.kinetic_energy() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.kinetic_energy(transform=transform) + npt.assert_array_equal(lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int_t, py_int, atol=atol, rtol=rtol) + + elif integral == "nuclear_attraction": + py_int = nuclear_electron_attraction_integral(py_basis, atcoords, atnums) + lc_int = lc_basis.nuclear_attraction() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.nuclear_attraction(transform=transform) + npt.assert_array_equal(lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int_t, py_int, atol=atol, rtol=rtol) + + elif integral == "rinv": + # Compare against the point_charge Python integral with a single + # unit charge at the origin, since rinv == 1/|r - origin| + origin = np.zeros(3) + py_int = point_charge_integral( + py_basis, origin.reshape(1, 3), np.asarray([-1.0]) + )[:, :, 0] + lc_int = lc_basis.rinv() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + # Test with inv_origin + lc_int_inv = lc_basis.rinv(inv_origin=np.zeros(3)) + npt.assert_allclose(lc_int_inv, py_int, atol=atol, rtol=rtol) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.rinv(transform=transform) + npt.assert_allclose(lc_int_t, py_int, atol=atol, rtol=rtol) + + elif integral == "dipole": + origin = np.zeros(3) + orders = np.asarray([[1, 0, 0]]) + py_int = moment_integral(py_basis, origin, orders, screen_basis=False)[:, :, 0] + lc_int = lc_basis.dipole() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.dipole(transform=transform) + npt.assert_array_equal(lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int_t, py_int, atol=atol, rtol=rtol) + + + elif integral == "quadrupole": + origin = np.zeros(3) + orders = np.asarray([[2, 0, 0]]) + py_int = moment_integral(py_basis, origin, orders, screen_basis=False)[:, :, 0] + lc_int = lc_basis.quadrupole() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.quadrupole(transform=transform) + npt.assert_array_equal(lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int_t, py_int, atol=atol, rtol=rtol) + + elif integral == "octupole": + origin = np.zeros(3) + orders = np.asarray([[3, 0, 0]]) + py_int = moment_integral(py_basis, origin, orders, screen_basis=False)[:, :, 0] + lc_int = lc_basis.octupole() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.octupole(transform=transform) + npt.assert_array_equal(lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int_t, py_int, atol=atol, rtol=rtol) + + elif integral == "momentum": + py_int = momentum_integral(py_basis, screen_basis=False) + lc_int = lc_basis.momentum(origin=np.zeros(3)) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, 3)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.momentum(origin=np.zeros(3), transform=transform) + npt.assert_array_equal(lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn, 3)) + npt.assert_allclose(lc_int_t, py_int, atol=atol, rtol=rtol) + + 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]) + lc_int = lc_basis.point_charge(charge_coords[:i], charges[:i]) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + + elif integral == "moment": + origin = np.zeros(3) + orders = np.asarray( + [ + [0, 0, 0], + [1, 0, 0], + [0, 1, 0], + [0, 0, 1], + [2, 0, 0], + [0, 2, 0], + [0, 0, 2], + [1, 1, 0], + [1, 0, 1], + [0, 1, 1], + [3, 0, 0], + [0, 3, 0], + [0, 0, 3], + ] + ) + py_int = moment_integral(py_basis, origin, orders, screen_basis=False) + lc_int = lc_basis.moment(orders, origin=origin) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) + npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) + + elif integral == "electron_repulsion": + py_int = electron_repulsion_integral_improved(py_basis) + lc_int = lc_basis.electron_repulsion(notation="physicist") + npt.assert_array_equal( + lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn) + ) + # ERI uses a looser tolerance, consistent with the existing + # electron_repulsion_integral test above + npt.assert_allclose(lc_int, py_int, atol=1e-4, rtol=1e-5) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.electron_repulsion(notation="physicist",transform=transform) + npt.assert_array_equal( + lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn, lc_basis.nbfn) + ) + npt.assert_allclose(lc_int_t, py_int, atol=1e-4, rtol=1e-5) + else: + raise ValueError(f"Invalid integral name '{integral}' passed") + + +@pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows") +@pytest.mark.skipif( + len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.*"))) == 0, + reason="The libcint shared library object was not found", +) +@pytest.mark.parametrize("atsyms, atcoords", TEST_SYSTEMS) +@pytest.mark.parametrize("basis", TEST_BASIS_SETS) + +def test_c_shellloop_matches_make_int1e(basis, atsyms, atcoords): + r""" + Cross-check the new C shell-loop bindings (``.overlap()``, + ``.kinetic_energy()``, ``.nuclear_attraction()``, ``.electron_repulsion()``) + directly against the existing ``make_int1e``/``make_int2e``-based methods + (``.overlap_integral()``, ``.kinetic_energy_integral()``, + ``.nuclear_attraction_integral()``, ``.electron_repulsion_integral()``) + on the *same* ``CBasis`` instance. + + This isolates the C shell-loop logic itself (PR-5/PR-6) from any + differences against the pure-Python GBasis implementation, since both + sides here come from libcint. + + """ + from gbasis.integrals.libcint import ELEMENTS, CBasis + + atcoords = atcoords / 0.5291772083 + + basis_dict = parse_nwchem(find_datafile(basis)) + + py_basis = make_contractions(basis_dict, atsyms, atcoords, coord_types="spherical") + + lc_basis = CBasis(py_basis, atsyms, atcoords, coord_type="spherical") + + # 1-electron integrals: C shell-loop vs. make_int1e shell-loop + npt.assert_allclose( + lc_basis.overlap(), lc_basis.overlap_integral(), atol=1e-10, rtol=1e-10 + ) + npt.assert_allclose( + lc_basis.kinetic_energy(), + lc_basis.kinetic_energy_integral(), + atol=1e-10, + rtol=1e-10, + ) + npt.assert_allclose( + lc_basis.nuclear_attraction(), + lc_basis.nuclear_attraction_integral(), + atol=1e-10, + rtol=1e-10, + ) + # overlap with transform + transform = np.eye(lc_basis.nbfn) + npt.assert_allclose( + lc_basis.overlap(transform=transform), + lc_basis.overlap_integral(transform=transform), + atol=1e-10, + rtol=1e-10, + ) + + # kinetic_energy with transform + npt.assert_allclose( + lc_basis.kinetic_energy(transform=transform), + lc_basis.kinetic_energy_integral(transform=transform), + atol=1e-10, rtol=1e-10, + ) + + # nuclear_attraction with transform + npt.assert_allclose( + lc_basis.nuclear_attraction(transform=transform), + lc_basis.nuclear_attraction_integral(transform=transform), + atol=1e-10, rtol=1e-10, + ) + + # rinv with inv_origin and transform + npt.assert_allclose( + lc_basis.rinv(inv_origin=np.zeros(3), transform=transform), + lc_basis.r_inv_integral(origin=np.zeros(3), transform=transform), + atol=1e-10, rtol=1e-10, + ) + + # momentum with transform + npt.assert_allclose( + lc_basis.momentum(origin=np.zeros(3), transform=transform), + lc_basis.momentum_integral(origin=np.zeros(3), transform=transform), + atol=1e-10, rtol=1e-10, + ) + + + # 2-electron ERI: C shell-loop vs. make_int2e shell-loop + npt.assert_allclose( + lc_basis.electron_repulsion(notation="chemist"), + lc_basis.electron_repulsion_integral(notation="chemist"), + atol=1e-8, + rtol=1e-8, + ) + + # electron_repulsion with transform + npt.assert_allclose( + lc_basis.electron_repulsion(notation="chemist", transform=transform), + lc_basis.electron_repulsion_integral(notation="chemist", transform=transform), + atol=1e-8, rtol=1e-8, + ) + +# ───────────────────────────────────────────────────────────────────────── +# Tests for gradient integral bindings +# ───────────────────────────────────────────────────────────────────────── + +TEST_GRADIENT_INTEGRALS = [ + pytest.param("gradient_kinetic", id="C-GradKinetic"), + pytest.param("gradient_nuclear", id="C-GradNuclear"), + pytest.param("gradient_rinv", id="C-GradRinv"), +] + + +@pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows") +@pytest.mark.skipif( + len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.*"))) == 0, + reason="The libcint shared library object was not found", +) +@pytest.mark.parametrize("integral", TEST_GRADIENT_INTEGRALS) +@pytest.mark.parametrize("atsyms, atcoords", TEST_SYSTEMS) +@pytest.mark.parametrize("basis", TEST_BASIS_SETS) +def test_c_gradient_integral(basis, atsyms, atcoords, integral): + r""" + Test the C shell-loop gradient integral bindings (PR-8) added to + ``gbasis.integrals.libcint.CBasis`` against the existing make_int1e + based implementations. + + These are the ``.gradient_kinetic()``, ``.gradient_nuclear()``, + and ``.gradient_rinv()`` methods which are the building blocks + for computing nuclear coordinate gradients. + """ + from gbasis.integrals.libcint import ELEMENTS, CBasis + + atcoords = atcoords / 0.5291772083 + + basis_dict = parse_nwchem(find_datafile(basis)) + + py_basis = make_contractions(basis_dict, atsyms, atcoords, coord_types="spherical") + lc_basis = CBasis(py_basis, atsyms, atcoords, coord_type="spherical") + + if integral == "gradient_kinetic": + # Compare C shell-loop against make_int1e path on same CBasis instance + py_int = lc_basis._d_kin() + lc_int = lc_basis.gradient_kinetic() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int[..., 0], atol=1e-10, rtol=1e-10) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.gradient_kinetic(transform=transform) + npt.assert_allclose(lc_int_t, py_int[..., 0], atol=1e-10, rtol=1e-10) + + elif integral == "gradient_nuclear": + py_int = lc_basis._d_nuc() + lc_int = lc_basis.gradient_nuclear() + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int[..., 0], atol=1e-10, rtol=1e-10) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.gradient_nuclear(transform=transform) + npt.assert_allclose(lc_int_t, py_int[..., 0], atol=1e-10, rtol=1e-10) + + elif integral == "gradient_rinv": + py_int = lc_basis._d_rinv(inv_origin=np.zeros(3)) + lc_int = lc_basis.gradient_rinv(inv_origin=np.zeros(3)) + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + npt.assert_allclose(lc_int, py_int[..., 0], atol=1e-10, rtol=1e-10) + # Test with transform + transform = np.eye(lc_basis.nbfn) + lc_int_t = lc_basis.gradient_rinv(inv_origin=np.zeros(3), transform=transform) + npt.assert_allclose(lc_int_t, py_int[..., 0], atol=1e-10, rtol=1e-10) + else: + raise ValueError(f"Invalid integral name '{integral}' passed") + + +# ───────────────────────────────────────────────────────────────────────── +# Tests for GIAO/magnetic integral bindings +# ───────────────────────────────────────────────────────────────────────── + +TEST_GIAO_INTEGRALS = [ + pytest.param("ia01p", id="C-ia01p"), + pytest.param("ircxp", id="C-ircxp"), + pytest.param("iking", id="C-iking"), + pytest.param("iovlpg", id="C-iovlpg"), + pytest.param("inucg", id="C-inucg"), +] + + +@pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows") +@pytest.mark.skipif( + len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.*"))) == 0, + reason="The libcint shared library object was not found", +) +@pytest.mark.parametrize("integral", TEST_GIAO_INTEGRALS) +@pytest.mark.parametrize("atsyms, atcoords", TEST_SYSTEMS) +@pytest.mark.parametrize("basis", TEST_BASIS_SETS) +def test_c_giao_integral(basis, atsyms, atcoords, integral): + r""" + Test the GIAO/magnetic integral bindings (PR-8) added to + ``gbasis.integrals.libcint.CBasis``. + + These are the ``.ia01p()``, ``.ircxp()``, ``.iking()``, + ``.iovlpg()``, and ``.inucg()`` methods which are building + blocks for NMR/magnetic property calculations. + + Since GBasis has no Python reference implementation for GIAO + integrals, we verify shape and that results are finite and + non-trivially zero for multi-atom systems. + """ + from gbasis.integrals.libcint import ELEMENTS, CBasis + + atcoords = atcoords / 0.5291772083 + + basis_dict = parse_nwchem(find_datafile(basis)) + py_basis = make_contractions(basis_dict, atsyms, atcoords, coord_types="spherical") + lc_basis = CBasis(py_basis, atsyms, atcoords, coord_type="spherical") + + if integral == "ia01p": + lc_int = lc_basis.ia01p() + elif integral == "ircxp": + lc_int = lc_basis.ircxp() + elif integral == "iking": + lc_int = lc_basis.iking() + elif integral == "iovlpg": + lc_int = lc_basis.iovlpg() + elif integral == "inucg": + lc_int = lc_basis.inucg() + else: + raise ValueError(f"Invalid integral name '{integral}' passed") + + # Shape check + npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) + + # Finiteness check — no NaN or Inf + assert np.all(np.isfinite(lc_int)), f"{integral} contains NaN or Inf" + # Test with transform + transform = np.eye(lc_basis.nbfn) + func = getattr(lc_basis, integral) + lc_int_t = func(transform=transform) + npt.assert_array_equal(lc_int_t.shape, (lc_basis.nbfn, lc_basis.nbfn)) + assert np.all(np.isfinite(lc_int_t)), f"{integral} with transform contains NaN or Inf" \ No newline at end of file