pyfock.Integrals

 7def mmd_nuc_mat_symm(basis, mol, slice=None):
 8    #Here the lists are converted to numpy arrays for better use with Numba.
 9    #Once these conversions are done we pass these to a Numba decorated
10    #function that uses prange, etc. to calculate the matrix efficiently.
11
12    # This function calculates the nuclear matrix for a given basis object.
13    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
14    # It is possible to only calculate a slice (block) of the complete matrix.
15    # slice is a 4 element list whose first and second elements give the range of the rows to be calculated.
16    # the third and fourth element give the range of columns to be calculated.
17    # slice = [start_row, end_row, start_col, end_col]
18    # The integrals are performed using the formulas
19
20    #We convert the required properties to numpy arrays as this is what Numba likes.
21    bfs_coords = np.array([basis.bfs_coords])
22    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
23    bfs_lmn = np.array([basis.bfs_lmn])
24    bfs_nprim = np.array([basis.bfs_nprim])
25    coordsBohrs = np.array([mol.coordsBohrs])
26    Z = np.array([mol.Zcharges])
27    natoms = mol.natoms
28        
29
30    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
31    #Numba won't be able to work with these efficiently.
32    #So, we convert them to a numpy 2d array by applying a trick,
33    #that the second dimension is that of the largest list. So that
34    #it can accomadate all the lists.
35    maxnprim = max(basis.bfs_nprim)
36    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
37    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
38    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
39    for i in range(basis.bfs_nao):
40        for j in range(basis.bfs_nprim[i]):
41            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
42            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
43            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
44        
45
46
47    if slice is None:
48        slice = [0, basis.bfs_nao, 0, basis.bfs_nao]
49        
50    #Limits for the calculation of overlap integrals
51    a = int(slice[0]) #row start index
52    b = int(slice[1]) #row end index
53    c = int(slice[2]) #column start index
54    d = int(slice[3]) #column end index
55
56    # print([a,b,c,d])
57    
58    V = mmd_nuc_mat_symm_internal(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts, a, b, c, d, Z[0], coordsBohrs[0], natoms)
59    
60    return V 

  7def nuc_mat_symm(basis, mol, slice=None, sqrt_ints4c2e_diag=None):
  8    """
  9    Compute the nuclear attraction matrix for a given basis and molecular geometry.
 10
 11    This function evaluates the one-electron nuclear attraction integrals 
 12    ⟨χ_i | V_nuc | χ_j⟩, where V_nuc is the Coulombic potential from the atomic nuclei.
 13    The integrals are computed between all basis functions defined in `basis` and 
 14    and the nuclei within the `mol` object.
 15    So in principle, it can be used to compute the nuclear potential matrix due to the
 16    nuclei in one molecule, in the basis of another molecule.
 17
 18    Efficient Numba-accelerated routines are used to improve performance. 
 19    A partial block of the matrix can be computed by specifying a `slice`.
 20
 21    Parameters
 22    ----------
 23    basis : object
 24        Basis set object with fields including:
 25        - bfs_coords: Cartesian coordinates of basis function centers.
 26        - bfs_coeffs: Contraction coefficients.
 27        - bfs_expnts: Gaussian exponents.
 28        - bfs_prim_norms: Normalization constants for primitives.
 29        - bfs_contr_prim_norms: Normalization constants for contracted functions.
 30        - bfs_lmn: Angular momentum quantum numbers (ℓ, m, n).
 31        - bfs_nprim: Number of primitives per basis function.
 32        - bfs_nao: Number of atomic orbitals (contracted basis functions).
 33
 34    mol : object
 35        Molecule object containing:
 36        - coordsBohrs: Cartesian coordinates of nuclei (in Bohr).
 37        - Zcharges: Nuclear charges.
 38        - natoms: Number of atoms.
 39
 40    slice : list of int, optional
 41        A 4-element list `[start_row, end_row, start_col, end_col]` defining 
 42        the sub-block of the matrix to compute. Rows and columns refer to AOs.
 43        If `None` (default), the entire nuclear attraction matrix is computed.
 44
 45    sqrt_ints4c2e_diag : ndarray, optional
 46        Optional Schwarz screening preconditioner.
 47        If provided, Schwarz screening is enabled to skip insignificant integrals.
 48        If `None` (default), all integrals are computed without screening.
 49
 50    Returns
 51    -------
 52    V : ndarray of shape (end_row - start_row, end_col - start_col)
 53        The computed (sub)matrix of nuclear attraction integrals.
 54
 55    Notes
 56    -----
 57    The integrals are evaluated using standard expressions over contracted Gaussians.
 58    If Schwarz screening is used (`sqrt_ints4c2e_diag` is not None), a threshold-based 
 59    pruning is applied using precomputed upper bounds.
 60
 61    Examples
 62    --------
 63    >>> V = nuc_mat_symm(basis, mol)
 64    >>> V_block = nuc_mat_symm(basis, mol, slice=[0, 10, 0, 10])
 65    >>> V_schwarz = nuc_mat_symm(basis, mol, sqrt_ints4c2e_diag=sqrt_diag)
 66    """
 67    #Here the lists are converted to numpy arrays for better use with Numba.
 68    #Once these conversions are done we pass these to a Numba decorated
 69    #function that uses prange, etc. to calculate the matrix efficiently.
 70
 71    # This function calculates the nuclear matrix for a given basis object.
 72    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
 73    # It is possible to only calculate a slice (block) of the complete matrix.
 74    # slice is a 4 element list whose first and second elements give the range of the rows to be calculated.
 75    # the third and fourth element give the range of columns to be calculated.
 76    # slice = [start_row, end_row, start_col, end_col]
 77    # The integrals are performed using the formulas
 78
 79    #We convert the required properties to numpy arrays as this is what Numba likes.
 80    bfs_coords = np.array([basis.bfs_coords])
 81    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
 82    bfs_lmn = np.array([basis.bfs_lmn])
 83    bfs_nprim = np.array([basis.bfs_nprim])
 84    coordsBohrs = np.array([mol.coordsBohrs])
 85    Z = np.array([mol.Zcharges])
 86    natoms = mol.natoms
 87        
 88
 89    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
 90    #Numba won't be able to work with these efficiently.
 91    #So, we convert them to a numpy 2d array by applying a trick,
 92    #that the second dimension is that of the largest list. So that
 93    #it can accomadate all the lists.
 94    maxnprim = max(basis.bfs_nprim)
 95    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
 96    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
 97    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
 98    for i in range(basis.bfs_nao):
 99        for j in range(basis.bfs_nprim[i]):
100            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
101            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
102            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
103        
104
105
106    if slice is None:
107        slice = [0, basis.bfs_nao, 0, basis.bfs_nao]
108        
109    #Limits for the calculation of overlap integrals
110    a = int(slice[0]) #row start index
111    b = int(slice[1]) #row end index
112    c = int(slice[2]) #column start index
113    d = int(slice[3]) #column end index
114
115    # print([a,b,c,d])
116
117    if sqrt_ints4c2e_diag is None:
118        #Create dummy array
119        sqrt_ints4c2e_diag = np.zeros((1,1), dtype=np.float64)
120        isSchwarz = False
121    else:
122        isSchwarz = True
123    
124    V = nuc_mat_symm_internal(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts, a, b, c, d, Z[0], coordsBohrs[0], natoms, sqrt_ints4c2e_diag, isSchwarz)
125    
126    return V 

 6def kin_mat_symm(basis, slice=None):
 7    """
 8    Compute the kinetic energy matrix for a given basis set.
 9
10    This function evaluates the one-electron kinetic energy integrals 
11    ⟨χ_i | -½∇² | χ_j⟩ for a set of contracted Gaussian basis functions defined in `basis`. 
12    It uses an optimized backend with improved memory layout, early termination, 
13    and partial vectorization for enhanced performance.
14
15    Block-wise computation is supported through the optional `slice` parameter.
16
17    Parameters
18    ----------
19    basis : object
20        Basis set object containing properties such as:
21        - bfs_coords: Cartesian coordinates of basis function centers.
22        - bfs_coeffs: Contraction coefficients.
23        - bfs_expnts: Gaussian exponents.
24        - bfs_prim_norms: Normalization constants for primitives.
25        - bfs_contr_prim_norms: Normalization constants for contracted functions.
26        - bfs_lmn: Angular momentum quantum numbers (ℓ, m, n).
27        - bfs_nprim: Number of primitives per basis function.
28        - bfs_nao: Number of atomic orbitals (contracted basis functions).
29
30    slice : list of int, optional
31        A 4-element list `[start_row, end_row, start_col, end_col]` specifying the 
32        matrix block to compute. Rows and columns refer to AOs. If `None` (default), 
33        the entire kinetic energy matrix is calculated.
34
35    Returns
36    -------
37    T : ndarray of shape (end_row - start_row, end_col - start_col)
38        The computed (sub)matrix of kinetic energy integrals.
39
40    Notes
41    -----
42    Internally, this function:
43    - Pads coefficient/exponent arrays for uniform shape and Numba compatibility.
44    - Reduces redundant operations and loops.
45    - Utilizes an optimized `kin_mat_symm_internal_optimized` backend.
46
47    Examples
48    --------
49    >>> T = kin_mat_symm(basis)
50    >>> T_block = kin_mat_symm(basis, slice=[0, 10, 0, 10])
51    """
52    
53    # Convert to numpy arrays (same as before)
54    bfs_coords = np.array(basis.bfs_coords)  # Remove unnecessary list wrapping
55    bfs_contr_prim_norms = np.array(basis.bfs_contr_prim_norms)
56    bfs_lmn = np.array(basis.bfs_lmn)
57    bfs_nprim = np.array(basis.bfs_nprim)
58    
59    # Optimize the coefficient/exponent arrays
60    maxnprim = max(basis.bfs_nprim)
61    bfs_coeffs = np.zeros((basis.bfs_nao, maxnprim))
62    bfs_expnts = np.zeros((basis.bfs_nao, maxnprim))
63    bfs_prim_norms = np.zeros((basis.bfs_nao, maxnprim))
64    
65    # Vectorized assignment (faster than nested loops)
66    for i in range(basis.bfs_nao):
67        n = basis.bfs_nprim[i]
68        bfs_coeffs[i, :n] = basis.bfs_coeffs[i]
69        bfs_expnts[i, :n] = basis.bfs_expnts[i]
70        bfs_prim_norms[i, :n] = basis.bfs_prim_norms[i]
71    
72    if slice is None:
73        slice = [0, basis.bfs_nao, 0, basis.bfs_nao]
74    
75    a, b, c, d = map(int, slice)
76    
77    T = kin_mat_symm_internal_optimized(
78        bfs_coords, bfs_contr_prim_norms, bfs_lmn, bfs_nprim, 
79        bfs_coeffs, bfs_prim_norms, bfs_expnts, a, b, c, d
80    )
81    
82    return T

 7def overlap_mat_symm(basis, slice=None):
 8    """
 9    Compute the overlap matrix for a given basis set using symmetry-aware integrals.
10
11    This function evaluates the overlap integrals ⟨χ_i | χ_j⟩ between all pairs 
12    of basis functions defined in the `basis` object. It supports partial evaluation 
13    of the overlap matrix by specifying a slice.
14
15    The integrals are computed using an efficient Numba-accelerated backend that 
16    benefits from parallelization via `prange` and preprocessed NumPy arrays.
17
18    Parameters
19    ----------
20    basis : object
21        A basis set object that contains information about basis functions, such as:
22        - bfs_coords: Cartesian coordinates of the basis function centers.
23        - bfs_coeffs: Contraction coefficients.
24        - bfs_expnts: Gaussian exponents.
25        - bfs_prim_norms: Primitive normalization constants.
26        - bfs_contr_prim_norms: Contraction normalization factors.
27        - bfs_lmn: Angular momentum quantum numbers (ℓ, m, n).
28        - bfs_nprim: Number of primitives per basis function.
29        - bfs_nao: Total number of atomic orbitals.
30    
31    slice : list of int, optional
32        A 4-element list specifying a sub-block of the matrix to compute:
33        [start_row, end_row, start_col, end_col].
34        If None (default), the full overlap matrix is computed.
35
36    Returns
37    -------
38    S : ndarray of shape (end_row - start_row, end_col - start_col)
39        The computed (sub)matrix of overlap integrals.
40
41    Notes
42    -----
43    This function is optimized for performance using preallocated NumPy arrays
44    and avoids nested Python lists which are not supported efficiently by Numba.
45
46    Examples
47    --------
48    >>> S = overlap_mat_symm(basis)
49    >>> S_block = overlap_mat_symm(basis, slice=[0, 10, 0, 10])
50    """
51    # Here the lists are converted to numpy arrays for better use with Numba.
52    # Once these conversions are done we pass these to a Numba decorated
53    # function that uses prange, etc. to calculate the matrix efficiently.
54
55    # This function calculates the overlap matrix for a given basis object.
56    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
57    # It is possible to only calculate a slice (block) of the complete matrix.
58    # slice is a 4 element list whose first and second elements give the range of the rows to be calculated.
59    # the third and fourth element give the range of columns to be calculated.
60    # slice = [start_row, end_row, start_col, end_col]
61    # The integrals are performed using the formulas
62
63    #We convert the required properties to numpy arrays as this is what Numba likes.
64    bfs_coords = np.array([basis.bfs_coords])
65    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
66    bfs_lmn = np.array([basis.bfs_lmn])
67    bfs_nprim = np.array([basis.bfs_nprim])
68        
69
70    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
71    #Numba won't be able to work with these efficiently.
72    #So, we convert them to a numpy 2d array by applying a trick,
73    #that the second dimension is that of the largest list. So that
74    #it can accomadate all the lists.
75    maxnprim = max(basis.bfs_nprim)
76    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
77    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
78    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
79    for i in range(basis.bfs_nao):
80        for j in range(basis.bfs_nprim[i]):
81            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
82            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
83            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
84        
85
86    if slice is None:
87        slice = [0, basis.bfs_nao, 0, basis.bfs_nao]
88        
89    #Limits for the calculation of overlap integrals
90    a = int(slice[0]) #row start index
91    b = int(slice[1]) #row end index
92    c = int(slice[2]) #column start index
93    d = int(slice[3]) #column end index
94        
95    S = overlap_mat_symm_internal(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts, a, b, c, d)
96    return S

  7def conv_4c2e_symm(basis, slice=None):
  8    """
  9    Compute four-center two-electron (4c2e) electron repulsion integrals (ERIs)
 10    using the conventional and slow (analytical) formula-based method with full symmetry exploitation.
 11
 12    This function evaluates integrals of the form (A B | C D), where
 13    A, B, C, D are basis functions from the same primary basis set.
 14
 15    The "conv" variant uses explicit analytical integral formulas and
 16    nested loops over primitive Gaussians, following the derivations in:
 17
 18        J. Chem. Educ. 2018, 95, 9, 1572–1578
 19        https://pubs.acs.org/doi/full/10.1021/acs.jchemed.8b00255
 20
 21    Compared to Rys quadrature, this conventional approach is much more
 22    computationally expensive but exact for the given formula set, making
 23    it useful for validation and special-purpose computations.
 24
 25    Symmetries exploited
 26    --------------------
 27    For 4c2e integrals, the full 8-fold permutational symmetry is used:
 28
 29        (A B | C D) = (B A | C D) = (A B | D C) = (B A | D C)
 30                    = (C D | A B) = (D C | A B) = (C D | B A) = (D C | B A)
 31
 32    This reduces the number of independent integrals from:
 33        N_bf^4   →   N_bf*(N_bf+1)/2 * N_bf*(N_bf+1)/2
 34    when computing the full tensor.
 35
 36    Parameters
 37    ----------
 38    basis : object
 39        Basis set object containing:
 40        - bfs_coords : Cartesian coordinates of basis function centers
 41        - bfs_coeffs : Contraction coefficients
 42        - bfs_expnts : Gaussian exponents
 43        - bfs_prim_norms : Primitive normalization constants
 44        - bfs_contr_prim_norms : Contraction normalization factors
 45        - bfs_lmn : Angular momentum quantum numbers (ℓ, m, n)
 46        - bfs_nprim : Number of primitives per basis function
 47        - bfs_nao : Total number of atomic orbitals
 48
 49    slice : list of int, optional
 50        An 8-element list specifying a sub-block of integrals to compute:
 51        [start_A, end_A, start_B, end_B, start_C, end_C, start_D, end_D]
 52        If None (default), computes the full (Nbf, Nbf, Nbf, Nbf) tensor.
 53
 54        **Note:** When slices are used, symmetry exploitation is restricted
 55        to permutations that lie entirely within the specified slice.
 56
 57    Returns
 58    -------
 59    ints4c2e : ndarray
 60        The computed 4-center 2-electron integrals for the requested range.
 61
 62    Notes
 63    -----
 64    - All basis set data are pre-packed into NumPy arrays for Numba acceleration.
 65    - Uses explicit primitive Gaussian formula evaluation without numerical quadrature.
 66    - Symmetry exploitation is maximal only for full-range computations.
 67    - Conventional evaluation scales as O(N_bf^4) without symmetry, but
 68      symmetry reduces the number of computations drastically.
 69
 70    Examples
 71    --------
 72    >>> eri_full = conv_4c2e_symm(basis)
 73    >>> eri_block = conv_4c2e_symm(basis, slice=[0,5, 0,5, 0,5, 0,5])
 74    """
 75    # Here the lists are converted to numpy arrays for better use with Numba.
 76    # Once these conversions are done we pass these to a Numba decorated
 77    # function that uses prange, etc. to calculate the 4c2e integrals efficiently.
 78    # This function calculates the 4c2e electron-electron ERIs for a given basis object.
 79    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
 80    
 81    # The integrals are performed using the formulas https://pubs.acs.org/doi/full/10.1021/acs.jchemed.8b00255
 82
 83    # It is possible to only calculate a slice (block/subset) of the complete set of integrals.
 84    # slice is an 8 element list whose first and second elements give the range of the A functions to be calculated.
 85    # and so on.
 86
 87    #We convert the required properties to numpy arrays as this is what Numba likes.
 88    bfs_coords = np.array([basis.bfs_coords])
 89    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
 90    bfs_lmn = np.array([basis.bfs_lmn])
 91    bfs_nprim = np.array([basis.bfs_nprim])
 92        
 93
 94    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
 95    #Numba won't be able to work with these efficiently.
 96    #So, we convert them to a numpy 2d array by applying a trick,
 97    #that the second dimension is that of the largest list. So that
 98    #it can accomadate all the lists.
 99    maxnprim = max(basis.bfs_nprim)
100    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
101    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
102    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
103    for i in range(basis.bfs_nao):
104        for j in range(basis.bfs_nprim[i]):
105            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
106            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
107            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
108        
109
110
111    
112    if slice is None:
113        slice = [0, basis.bfs_nao, 0, basis.bfs_nao, 0, basis.bfs_nao, 0, basis.bfs_nao]
114        
115    #Limits for the calculation of 4c2e integrals
116    indx_startA = int(slice[0])
117    indx_endA = int(slice[1])
118    indx_startB = int(slice[2])
119    indx_endB = int(slice[3])
120    indx_startC = int(slice[4])
121    indx_endC = int(slice[5])
122    indx_startD = int(slice[6])
123    indx_endD = int(slice[7])
124
125    ints4c2e = conv_4c2e_symm_internal(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts,\
126        indx_startA,indx_endA, indx_startB, indx_endB, indx_startC, indx_endC, indx_startD,indx_endD)
127
128    return ints4c2e

 7def mmd_4c2e_symm(basis, slice=None):
 8    # Here the lists are converted to numpy arrays for better use with Numba.
 9    # Once these conversions are done we pass these to a Numba decorated
10    # function that uses prange, etc. to calculate the 4c2e integrals efficiently.
11    # This function calculates the 4c2e electron-electron ERIs for a given basis object.
12    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
13    
14    # The integrals are performed using the formulas https://pubs.acs.org/doi/full/10.1021/acs.jchemed.8b00255
15
16    # It is possible to only calculate a slice (block/subset) of the complete set of integrals.
17    # slice is an 8 element list whose first and second elements give the range of the A functions to be calculated.
18    # and so on.
19
20    #We convert the required properties to numpy arrays as this is what Numba likes.
21    bfs_coords = np.array([basis.bfs_coords])
22    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
23    bfs_lmn = np.array([basis.bfs_lmn])
24    bfs_nprim = np.array([basis.bfs_nprim])
25        
26
27    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
28    #Numba won't be able to work with these efficiently.
29    #So, we convert them to a numpy 2d array by applying a trick,
30    #that the second dimension is that of the largest list. So that
31    #it can accomadate all the lists.
32    maxnprim = max(basis.bfs_nprim)
33    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
34    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
35    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
36    for i in range(basis.bfs_nao):
37        for j in range(basis.bfs_nprim[i]):
38            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
39            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
40            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
41        
42
43
44    
45    if slice is None:
46        slice = [0, basis.bfs_nao, 0, basis.bfs_nao, 0, basis.bfs_nao, 0, basis.bfs_nao]
47        
48    #Limits for the calculation of 4c2e integrals
49    indx_startA = int(slice[0])
50    indx_endA = int(slice[1])
51    indx_startB = int(slice[2])
52    indx_endB = int(slice[3])
53    indx_startC = int(slice[4])
54    indx_endC = int(slice[5])
55    indx_startD = int(slice[6])
56    indx_endD = int(slice[7])
57
58    ints4c2e = mmd_4c2e_symm_internal2(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts,\
59        indx_startA,indx_endA, indx_startB, indx_endB, indx_startC, indx_endC, indx_startD,indx_endD)
60
61    return ints4c2e

  8def rys_4c2e_symm(basis, slice=None):
  9    """
 10    Compute four-center two-electron (4c2e) electron repulsion integrals (ERIs) 
 11    using the Rys quadrature method with exploitation of 8-fold permutational 
 12    symmetry.
 13
 14    This function evaluates integrals of the form (A B | C D), where A, B, C, D 
 15    are basis functions from the same primary basis set. It uses Numba-accelerated 
 16    backends and symmetry-aware optimizations to reduce computational cost.
 17
 18    Symmetries exploited
 19    --------------------
 20    The 4c2e integrals obey the following permutational symmetry relations:
 21
 22        (A B | C D) = (B A | C D)   # exchange within bra
 23                    = (A B | D C)  # exchange within ket
 24                    = (B A | D C)  
 25                    = (C D | A B)  # bra ↔ ket exchange
 26                    = (D C | A B)
 27                    = (C D | B A)
 28                    = (D C | B A)
 29
 30    These 8 equivalent permutations mean that only a subset of integrals 
 31    needs to be explicitly computed; the rest can be obtained by symmetry.
 32    
 33    This reduces the number of independent integrals from:
 34        N_bf^4   →   N_bf*(N_bf+1)/2 * N_bf*(N_bf+1)/2
 35    when computing the full tensor.
 36
 37    Parameters
 38    ----------
 39    basis : object
 40        Primary basis set object containing:
 41        - bfs_coords : Cartesian coordinates of basis function centers.
 42        - bfs_coeffs : Contraction coefficients.
 43        - bfs_expnts : Gaussian exponents.
 44        - bfs_prim_norms : Primitive normalization constants.
 45        - bfs_contr_prim_norms : Contraction normalization factors.
 46        - bfs_lmn : Angular momentum quantum numbers (ℓ, m, n).
 47        - bfs_nprim : Number of primitives per basis function.
 48        - bfs_shell_index : Index of the shell each basis function belongs to.
 49        - bfs_nao : Total number of atomic orbitals.
 50
 51    slice : list of int, optional
 52        An 8-element list specifying a sub-block of integrals to compute:
 53        [start_A, end_A, start_B, end_B, start_C, end_C, start_D, end_D]
 54        If None (default), computes the full (Nbf, Nbf, Nbf, Nbf) block.
 55
 56        **Note:** When slices are used, not all 8-fold symmetries may be 
 57        available because the requested block may not contain all 
 58        symmetric permutations. In such cases, only the applicable 
 59        symmetries are used.
 60
 61    Returns
 62    -------
 63    ints4c2e : ndarray
 64        The computed 4-center 2-electron integrals for the requested range.
 65
 66    Notes
 67    -----
 68    - Precomputes and stores basis set data in NumPy arrays for Numba efficiency.
 69    - Exploits all possible symmetry permutations when the full tensor is 
 70      computed (no slice) to reduce redundant work.
 71    - If a slice is specified, symmetry exploitation is limited to the 
 72      permutations that fall within the slice's index ranges.
 73
 74    Examples
 75    --------
 76    >>> eri_full = rys_4c2e_symm(basis)
 77    >>> eri_block = rys_4c2e_symm(basis, slice=[0,5, 0,5, 0,5, 0,5])
 78    """
 79    # Here the lists are converted to numpy arrays for better use with Numba.
 80    # Once these conversions are done we pass these to a Numba decorated
 81    # function that uses prange, etc. to calculate the 4c2e integrals efficiently.
 82    # This function calculates the 4c2e electron-electron ERIs for a given basis object.
 83    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
 84
 85    # It is possible to only calculate a slice (block/subset) of the complete set of integrals.
 86    # slice is an 8 element list whose first and second elements give the range of the A functions to be calculated.
 87    # and so on.
 88
 89    #We convert the required properties to numpy arrays as this is what Numba likes.
 90    bfs_coords = np.array([basis.bfs_coords])
 91    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
 92    bfs_lmn = np.array([basis.bfs_lmn])
 93    bfs_nprim = np.array([basis.bfs_nprim])
 94        
 95
 96    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
 97    #Numba won't be able to work with these efficiently.
 98    #So, we convert them to a numpy 2d array by applying a trick,
 99    #that the second dimension is that of the largest list. So that
100    #it can accomadate all the lists.
101    maxnprim = max(basis.bfs_nprim)
102    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
103    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
104    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
105    shell_indices = np.array([basis.bfs_shell_index], dtype=np.uint16)[0]
106    for i in range(basis.bfs_nao):
107        for j in range(basis.bfs_nprim[i]):
108            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
109            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
110            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
111        
112
113
114    
115    if slice is None:
116        slice = [0, basis.bfs_nao, 0, basis.bfs_nao, 0, basis.bfs_nao, 0, basis.bfs_nao]
117        
118    #Limits for the calculation of 4c2e integrals
119    indx_startA = int(slice[0])
120    indx_endA = int(slice[1])
121    indx_startB = int(slice[2])
122    indx_endB = int(slice[3])
123    indx_startC = int(slice[4])
124    indx_endC = int(slice[5])
125    indx_startD = int(slice[6])
126    indx_endD = int(slice[7])
127
128
129    ints4c2e = rys_4c2e_symm_internal(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts, shell_indices,\
130        indx_startA,indx_endA, indx_startB, indx_endB, indx_startC, indx_endC, indx_startD,indx_endD)
131
132    return ints4c2e

416def rys_4c2e_symm_old(basis, slice=None):
417    # Here the lists are converted to numpy arrays for better use with Numba.
418    # Once these conversions are done we pass these to a Numba decorated
419    # function that uses prange, etc. to calculate the 4c2e integrals efficiently.
420    # This function calculates the 4c2e electron-electron ERIs for a given basis object.
421    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
422
423    # It is possible to only calculate a slice (block/subset) of the complete set of integrals.
424    # slice is an 8 element list whose first and second elements give the range of the A functions to be calculated.
425    # and so on.
426
427    #We convert the required properties to numpy arrays as this is what Numba likes.
428    bfs_coords = np.array([basis.bfs_coords])
429    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
430    bfs_lmn = np.array([basis.bfs_lmn])
431    bfs_nprim = np.array([basis.bfs_nprim])
432        
433
434    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
435    #Numba won't be able to work with these efficiently.
436    #So, we convert them to a numpy 2d array by applying a trick,
437    #that the second dimension is that of the largest list. So that
438    #it can accomadate all the lists.
439    maxnprim = max(basis.bfs_nprim)
440    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
441    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
442    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
443    for i in range(basis.bfs_nao):
444        for j in range(basis.bfs_nprim[i]):
445            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
446            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
447            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
448        
449
450
451    
452    if slice is None:
453        slice = [0, basis.bfs_nao, 0, basis.bfs_nao, 0, basis.bfs_nao, 0, basis.bfs_nao]
454        
455    #Limits for the calculation of 4c2e integrals
456    indx_startA = int(slice[0])
457    indx_endA = int(slice[1])
458    indx_startB = int(slice[2])
459    indx_endB = int(slice[3])
460    indx_startC = int(slice[4])
461    indx_endC = int(slice[5])
462    indx_startD = int(slice[6])
463    indx_endD = int(slice[7])
464
465    ints4c2e = rys_4c2e_symm_internal_old(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts,\
466        indx_startA,indx_endA, indx_startB, indx_endB, indx_startC, indx_endC, indx_startD,indx_endD)
467
468    return ints4c2e

  7def conv_3c2e_symm(basis, auxbasis, slice=None):
  8    """
  9    Compute three-center two-electron (3c2e) electron repulsion integrals (ERIs)
 10    using the conventional and slow (analytical) formula-based method with symmetry exploitation.
 11
 12    This function evaluates integrals of the form (A B | C), where:
 13        - A and B are primary basis functions from `basis`
 14        - C is an auxiliary basis function from `auxbasis`
 15
 16    The "conv" variant uses explicit analytical integral formulas and nested
 17    loops over primitive Gaussians, following the derivations in:
 18
 19        J. Chem. Educ. 2018, 95, 9, 1572–1578
 20        https://pubs.acs.org/doi/full/10.1021/acs.jchemed.8b00255
 21
 22    Compared to the Rys quadrature method, this conventional approach is
 23    significantly more computationally expensive, but can be useful for
 24    validation or when Rys is not applicable.
 25
 26    Symmetries exploited
 27    --------------------
 28    For 3c2e integrals, the following bra symmetry is used:
 29
 30        (A B | C) = (B A | C)
 31
 32    This reduces the number of computed integrals from N_bf² * N_aux
 33    to N_bf*(N_bf+1)/2 * N_aux when the full tensor is computed.
 34
 35    Parameters
 36    ----------
 37    basis : object
 38        Primary basis set object containing:
 39        - bfs_coords : Cartesian coordinates of basis function centers
 40        - bfs_coeffs : Contraction coefficients
 41        - bfs_expnts : Gaussian exponents
 42        - bfs_prim_norms : Primitive normalization constants
 43        - bfs_contr_prim_norms : Contraction normalization factors
 44        - bfs_lmn : Angular momentum quantum numbers (ℓ, m, n)
 45        - bfs_nprim : Number of primitives per basis function
 46        - bfs_nao : Total number of atomic orbitals
 47
 48    auxbasis : object
 49        Auxiliary basis set object with the same attributes as `basis`.
 50
 51    slice : list of int, optional
 52        A 6-element list specifying a sub-block of integrals to compute:
 53        [start_A, end_A, start_B, end_B, start_C, end_C]
 54        If None (default), computes the full (Nbf, Nbf, Naux) tensor.
 55
 56        **Note:** When slices are used, AB symmetry exploitation is limited
 57        to permutations that lie entirely within the specified slice.
 58
 59    Returns
 60    -------
 61    ints3c2e : ndarray
 62        The computed 3-center 2-electron integrals for the requested range.
 63        Shape: (Nbf, Nbf, Nauxbf) or 
 64        (end_A - start_A, end_B - start_B, end_C - start_C) 
 65        if slice is given.
 66
 67    Notes
 68    -----
 69    - All basis set data are pre-packed into NumPy arrays for Numba acceleration.
 70    - Uses explicit primitive Gaussian formula evaluation without numerical quadrature.
 71    - Symmetry exploitation is maximal only for full-range computations.
 72    - Conventional evaluation scales poorly compared to Rys quadrature, but is
 73      exact for the given formula set.
 74
 75    Examples
 76    --------
 77    >>> eri_full = conv_3c2e_symm(basis, auxbasis)
 78    >>> eri_block = conv_3c2e_symm(basis, auxbasis, slice=[0,5, 0,5, 0,10])
 79    """
 80    # Here the lists are converted to numpy arrays for better use with Numba.
 81    # Once these conversions are done we pass these to a Numba decorated
 82    # function that uses prange, etc. to calculate the 3c2e integrals efficiently.
 83    # This function calculates the 3c2e electron-electron ERIs for a given basis object and auxbasis object.
 84    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
 85    
 86    # The integrals are performed using the formulas https://pubs.acs.org/doi/full/10.1021/acs.jchemed.8b00255
 87
 88    # It is possible to only calculate a slice (block/subset) of the complete set of integrals.
 89    # slice is an 6 element list whose first and second elements give the range of the A functions to be calculated.
 90    # and so on.
 91    # slice = [indx_startA, indx_endA, indx_startB, indx_endB, indx_startC, indx_endC]
 92
 93    #We convert the required properties to numpy arrays as this is what Numba likes.
 94    bfs_coords = np.array([basis.bfs_coords])
 95    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
 96    bfs_lmn = np.array([basis.bfs_lmn])
 97    bfs_nprim = np.array([basis.bfs_nprim])
 98
 99    #We convert the required properties to numpy arrays as this is what Numba likes.
100    aux_bfs_coords = np.array([auxbasis.bfs_coords])
101    aux_bfs_contr_prim_norms = np.array([auxbasis.bfs_contr_prim_norms])
102    aux_bfs_lmn = np.array([auxbasis.bfs_lmn])
103    aux_bfs_nprim = np.array([auxbasis.bfs_nprim])
104        
105
106    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
107    #Numba won't be able to work with these efficiently.
108    #So, we convert them to a numpy 2d array by applying a trick,
109    #that the second dimension is that of the largest list. So that
110    #it can accomadate all the lists.
111    maxnprim = max(basis.bfs_nprim)
112    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
113    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
114    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
115    for i in range(basis.bfs_nao):
116        for j in range(basis.bfs_nprim[i]):
117            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
118            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
119            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
120
121    maxnprimaux = max(auxbasis.bfs_nprim)
122    aux_bfs_coeffs = np.zeros([auxbasis.bfs_nao, maxnprimaux])
123    aux_bfs_expnts = np.zeros([auxbasis.bfs_nao, maxnprimaux])
124    aux_bfs_prim_norms = np.zeros([auxbasis.bfs_nao, maxnprimaux])
125    for i in range(auxbasis.bfs_nao):
126        for j in range(auxbasis.bfs_nprim[i]):
127            aux_bfs_coeffs[i,j] = auxbasis.bfs_coeffs[i][j]
128            aux_bfs_expnts[i,j] = auxbasis.bfs_expnts[i][j]
129            aux_bfs_prim_norms[i,j] = auxbasis.bfs_prim_norms[i][j]
130            
131
132
133        
134    if slice is None:
135        slice = [0, basis.bfs_nao, 0, basis.bfs_nao, 0, auxbasis.bfs_nao]
136        
137    #Limits for the calculation of 4c2e integrals
138    indx_startA = int(slice[0])
139    indx_endA = int(slice[1])
140    indx_startB = int(slice[2])
141    indx_endB = int(slice[3])
142    indx_startC = int(slice[4])
143    indx_endC = int(slice[5])
144
145    ints3c2e = conv_3c2e_symm_internal(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts,aux_bfs_coords[0], aux_bfs_contr_prim_norms[0], aux_bfs_lmn[0], aux_bfs_nprim[0], aux_bfs_coeffs, aux_bfs_prim_norms, aux_bfs_expnts,indx_startA, indx_endA, indx_startB, indx_endB, indx_startC, indx_endC)
146    return ints3c2e

 11def rys_3c2e_symm(basis, auxbasis, slice=None, schwarz=False, schwarz_threshold=1e-9):
 12    """
 13    Compute three-center two-electron (3c2e) electron repulsion integrals using
 14    the Rys quadrature method with symmetry considerations.
 15
 16    This function evaluates integrals of the form (A B | C), where A and B are
 17    basis functions from a primary basis set and C is from an auxiliary basis set.
 18    Symmetry in the first two indices is exploited ((A B | C) = (B A | C)),
 19    so only the upper-triangular part of the (A,B) block is computed and the (B A | C)
 20    integrals are formed by symmetry.
 21    The implementation uses Numba-accelerated backends with symmetry-aware
 22    optimizations and optional Schwarz screening to skip negligible contributions.
 23    Symmetry is utilized to only compute N_{bf}*(N_{bf}+1)/2*N_{auxbf}
 24
 25    Parameters
 26    ----------
 27    basis : object
 28        Primary basis set object containing information about atomic orbitals, such as:
 29        - bfs_coords : Cartesian coordinates of basis function centers.
 30        - bfs_coeffs : Contraction coefficients.
 31        - bfs_expnts : Gaussian exponents.
 32        - bfs_prim_norms : Primitive normalization constants.
 33        - bfs_contr_prim_norms : Contraction normalization factors.
 34        - bfs_lmn : Angular momentum quantum numbers (ℓ, m, n).
 35        - bfs_nprim : Number of primitives per basis function.
 36        - bfs_nao : Total number of atomic orbitals.
 37
 38    auxbasis : object
 39        Auxiliary basis set object with the same attributes as `basis`, 
 40        but typically used for resolution-of-the-identity (RI) expansions.
 41
 42    slice : list of int, optional
 43        A 6-element list specifying a subset of integrals to compute:
 44        [start_A, end_A, start_B, end_B, start_C, end_C]
 45        If None (default), computes all N_{bf}*N_{bf}*N_{auxbf} integrals.
 46
 47    schwarz : bool, optional
 48        If True, applies Schwarz screening to skip calculations where the 
 49        product of integral bounds is below `schwarz_threshold`.
 50
 51    schwarz_threshold : float, optional
 52        The threshold for Schwarz screening (default is 1e-9).
 53
 54    Returns
 55    -------
 56    ints3c2e : ndarray of shape 
 57        (Nbf, Nbf, Nauxbf) or 
 58        (end_A - start_A, end_B - start_B, end_C - start_C) 
 59        if slice is given.
 60        The computed 3-center 2-electron integrals.
 61
 62    Notes
 63    -----
 64    - Uses preallocated NumPy arrays to store primitive data for efficient Numba processing.
 65    - Handles irregular contraction patterns by padding primitive arrays to the 
 66      size of the largest contraction.
 67    - If Schwarz screening is enabled, precomputes diagonal two-electron integrals 
 68      and uses their square roots for screening.
 69    - Symmetry relations are exploited to avoid redundant calculations.
 70
 71    Examples
 72    --------
 73    >>> ints = rys_3c2e_symm(basis, auxbasis)
 74    >>> ints_block = rys_3c2e_symm(basis, auxbasis, slice=[0, 5, 0, 5, 0, 10])
 75    >>> ints_screened = rys_3c2e_symm(basis, auxbasis, schwarz=True, schwarz_threshold=1e-10)
 76    """
 77    # Here the lists are converted to numpy arrays for better use with Numba.
 78    # Once these conversions are done we pass these to a Numba decorated
 79    # function that uses prange, etc. to calculate the 3c2e integrals efficiently.
 80    # This function calculates the 3c2e electron-electron ERIs for a given basis object and auxbasis object.
 81    # The basis object holds the information of basis functions like: exponents, coeffs, etc.
 82    
 83
 84    # It is possible to only calculate a slice (block/subset) of the complete set of integrals.
 85    # slice is a 6 element list whose first and second elements give the range of the A functions to be calculated.
 86    # and so on.
 87    # slice = [indx_startA, indx_endA, indx_startB, indx_endB, indx_startC, indx_endC]
 88
 89    #We convert the required properties to numpy arrays as this is what Numba likes.
 90    bfs_coords = np.array([basis.bfs_coords])
 91    bfs_contr_prim_norms = np.array([basis.bfs_contr_prim_norms])
 92    bfs_lmn = np.array([basis.bfs_lmn])
 93    bfs_nprim = np.array([basis.bfs_nprim])
 94
 95    #We convert the required properties to numpy arrays as this is what Numba likes.
 96    aux_bfs_coords = np.array([auxbasis.bfs_coords])
 97    aux_bfs_contr_prim_norms = np.array([auxbasis.bfs_contr_prim_norms])
 98    aux_bfs_lmn = np.array([auxbasis.bfs_lmn])
 99    aux_bfs_nprim = np.array([auxbasis.bfs_nprim])
100        
101
102    #The remaining properties like bfs_coeffs are a list of lists of unequal sizes.
103    #Numba won't be able to work with these efficiently.
104    #So, we convert them to a numpy 2d array by applying a trick,
105    #that the second dimension is that of the largest list. So that
106    #it can accomadate all the lists.
107    maxnprim = max(basis.bfs_nprim)
108    bfs_coeffs = np.zeros([basis.bfs_nao, maxnprim])
109    bfs_expnts = np.zeros([basis.bfs_nao, maxnprim])
110    bfs_prim_norms = np.zeros([basis.bfs_nao, maxnprim])
111    for i in range(basis.bfs_nao):
112        for j in range(basis.bfs_nprim[i]):
113            bfs_coeffs[i,j] = basis.bfs_coeffs[i][j]
114            bfs_expnts[i,j] = basis.bfs_expnts[i][j]
115            bfs_prim_norms[i,j] = basis.bfs_prim_norms[i][j]
116
117    maxnprimaux = max(auxbasis.bfs_nprim)
118    aux_bfs_coeffs = np.zeros([auxbasis.bfs_nao, maxnprimaux])
119    aux_bfs_expnts = np.zeros([auxbasis.bfs_nao, maxnprimaux])
120    aux_bfs_prim_norms = np.zeros([auxbasis.bfs_nao, maxnprimaux])
121    for i in range(auxbasis.bfs_nao):
122        for j in range(auxbasis.bfs_nprim[i]):
123            aux_bfs_coeffs[i,j] = auxbasis.bfs_coeffs[i][j]
124            aux_bfs_expnts[i,j] = auxbasis.bfs_expnts[i][j]
125            aux_bfs_prim_norms[i,j] = auxbasis.bfs_prim_norms[i][j]
126            
127
128
129        
130    if slice is None:
131        slice = [0, basis.bfs_nao, 0, basis.bfs_nao, 0, auxbasis.bfs_nao]
132        
133    #Limits for the calculation of 4c2e integrals
134    indx_startA = int(slice[0])
135    indx_endA = int(slice[1])
136    indx_startB = int(slice[2])
137    indx_endB = int(slice[3])
138    indx_startC = int(slice[4])
139    indx_endC = int(slice[5])
140
141    if schwarz:
142        ints4c2e_diag = eri_4c2e_diag(basis)
143        ints2c2e = rys_2c2e_symm(auxbasis)
144        sqrt_ints4c2e_diag = np.sqrt(np.abs(ints4c2e_diag))
145        sqrt_diag_ints2c2e = np.sqrt(np.abs(np.diag(ints2c2e)))
146        print('Prelims calc done for Schwarz screening!')
147    else:
148        #Create dummy array
149        sqrt_ints4c2e_diag = np.zeros((1,1), dtype=np.float64)
150        sqrt_diag_ints2c2e = np.zeros((1), dtype=np.float64)
151
152    ints3c2e = rys_3c2e_symm_internal(bfs_coords[0], bfs_contr_prim_norms[0], bfs_lmn[0], bfs_nprim[0], bfs_coeffs, bfs_prim_norms, bfs_expnts,aux_bfs_coords[0], aux_bfs_contr_prim_norms[0], aux_bfs_lmn[0], aux_bfs_nprim[0], aux_bfs_coeffs, aux_bfs_prim_norms, aux_bfs_expnts,indx_startA, indx_endA, indx_startB, indx_endB, indx_startC, indx_endC, schwarz, sqrt_ints4c2e_diag, sqrt_diag_ints2c2e, schwarz_threshold)
153    return ints3c2e