GBasis Tutorial - Integrals#
gbasis supports computation of 1- and 2-electron integrals. See the using the gbasis.integrals module for the complete list of integrals supported.
This notebook provides examples for computing various integrals for the helium hydride ion (\(\text{HeH}^{+}\)) using the 6-31G basis set, and 6-311G when two basis sets are needed. When available reference values are included from the calculation’s output file for comparison.
Building Basis Set#
Here, the basis set is built from a formatted checkpoint file (of a HF/6-31G calculation for the helium hydride ion) which is loaded using the IOData package. For a detailed tutorial on building basis sets, see the Building Basis Set tutorial.
from iodata import load_one
from gbasis.wrappers import from_iodata
# load FCHK file and build basis set in atomic orbitals (AO) basis
# Source:
mol = load_one("HeHp.fchk")
ao_basis = from_iodata(mol)
print(f"Loaded: {mol.title}")
print(f"Basis: {mol.obasis_name}")
print(f"Atomic numbers: {mol.atnums}")
print(f"Molecular charge: {mol.charge}")
print(f"Atomic coordinates:\n{mol.atcoords}")
print(f"Number of contraction shells: {len(ao_basis)}")
Loaded: HeH+ ion experimental R
Basis: 6-31g
Atomic numbers: [1 2]
Molecular charge: 1.0
Atomic coordinates:
[[ 3.12913779e-17 0.00000000e+00 -9.72579045e-01]
[ 3.10359322e-17 0.00000000e+00 4.86289523e-01]]
Number of contraction shells: 4
1. One-Electron Integrals#
1.1. Overlap Integrals#
Given a list of generalized contraction shells, the overlap matrix is computed by overlap_integral function. This result is the AO overlap integrals, unless the optional transform matrix is provided, which transforms the basis functions into a different basis set (e.g., from the atomic orbital basis to the molecular orbital basis) before computing the overlap integrals.
# From Gaussian16 Output File:
# Overlap matrix (S) of atomic orbitals:
# 1 0.100000D+01
# 2 0.658292D+00 0.100000D+01
# 3 0.300303D+00 0.334248D+00 0.100000D+01
# 4 0.536094D+00 0.746564D+00 0.634148D+00 0.100000D+01
import numpy as np
from gbasis.integrals.overlap import overlap_integral
# compute overlap of atomic orbitals (AO)
int1e_overlap_ao = overlap_integral(ao_basis)
print("Overlap matrix (S) of atomic orbitals:")
print(int1e_overlap_ao, end="\n\n")
# check whether overlap matrix is identity matrix (i.e., atomic orbitals are orthonormal)
is_orthogonal = np.allclose(int1e_overlap_ao, np.eye(int1e_overlap_ao.shape[0]))
print(f"Are AOs orthogonal? {is_orthogonal}")
Overlap matrix (S) of atomic orbitals:
[[1. 0.65829197 0.30030322 0.53609401]
[0.65829197 1. 0.33424849 0.74656436]
[0.30030322 0.33424849 1. 0.63414774]
[0.53609401 0.74656436 0.63414774 1. ]]
Are AOs orthogonal? False
Example: Overlap Integrals in MO Basis#
The transformation matrix from AO to MO basis is the transposed molecular orbital coefficient matrix. Providing this transformation as the transform argument to the overlap_integral function, results in the overlap integrals in the MO basis (i.e., transformed basis).
# get MO coefficients from IOData object
mo_coeffs = mol.mo.coeffs
# compute overlap of molecular orbitals (MO)
int1e_overlap_mo = overlap_integral(ao_basis, transform=mo_coeffs.T)
print("Overlap matrix (S) of molecular orbitals:")
print(int1e_overlap_mo, end="\n\n")
# check whether overlap matrix is identity matrix (i.e., molecular orbitals are orthonormal)
is_orthogonal = np.allclose(int1e_overlap_mo, np.eye(int1e_overlap_mo.shape[0]))
print(f"Are MOs orthogonal? {is_orthogonal}")
Overlap matrix (S) of molecular orbitals:
[[ 1.00000000e+00 2.38057226e-09 -1.60352739e-09 1.18604386e-09]
[ 2.38057211e-09 1.00000000e+00 -2.64133088e-09 1.11497621e-09]
[-1.60352734e-09 -2.64133095e-09 9.99999999e-01 -2.03584584e-09]
[ 1.18604397e-09 1.11497609e-09 -2.03584581e-09 1.00000000e+00]]
Are MOs orthogonal? True
1.2. Overlap Integrals Between Two Different Basis Sets#
The overlap_integral_asymmetric function computes the overlap integrals between two different basis sets, each denoted by a list of generalized contraction shells. The following example showcase how to use this feature to compute the overlap integrals between the 6-31G and 6-311G basis sets:
from gbasis.parsers import parse_gbs, make_contractions
from gbasis.integrals.overlap_asymm import overlap_integral_asymmetric
# create an 6-311G basis set for the helium hydride ion in spherical coordinates
basis_dict = parse_gbs("HeH_6-311g.gbs")
ao_basis_new = make_contractions(basis_dict, ["H", "He"], mol.atcoords, "p")
print(f"Number of shells in 6-311G basis: {len(ao_basis_new)}")
print(f"Number of shells in 6-31G basis: {len(ao_basis)}", end="\n\n")
# compute overlap of two different basis sets
int1e_overlap_basis = overlap_integral_asymmetric(ao_basis_new, ao_basis)
print(f"Shape of overlap matrix: {int1e_overlap_basis.shape}")
print("Overlap matrix (S) of atomic orbitals between old and new basis:")
print(int1e_overlap_basis)
Number of shells in 6-311G basis: 6
Number of shells in 6-31G basis: 4
Shape of overlap matrix: (6, 4)
Overlap matrix (S) of atomic orbitals between old and new basis:
[[0.96216463 0.49156933 0.22083097 0.40624233]
[0.86422289 0.91312074 0.37721246 0.7169519 ]
[0.52114949 0.96290073 0.28035313 0.69362244]
[0.15725291 0.17901201 0.90190087 0.37318667]
[0.40988749 0.46228859 0.92666879 0.81268342]
[0.53245911 0.78732494 0.57381839 0.9927251 ]]
Example: Project MO Coefficients from One Basis Set to Another#
Given a set of HF molecular orbitals (MO) \(\psi_{i}\) and basis functions (AOs) \(\phi_{\mu}\), they are related by the MO coefficient matrix \(C_{\mu i}\) by:
where \(C_{\mu i}\) can be seen as the transformation matrix from AOs to MOs. These MO coefficients \(\{C_{\mu i}\}\) are the elements of the overlap matrix between the MO basis \(\{\psi\}\) and the AO basis \(\{\phi\}\).
Given a different set of basis functions \(\chi_{\nu}\), the MO coefficients in the new basis can be expressed (projected) as:
where \(\phi_{\mu}^{\text{basis 2}}\) are the basis functions in the new basis set. This can be seen as the overlap matrix between the (transformed) MO basis \(\{\psi\}\) and the new basis set \(\{\phi^{\text{basis 2}}\}\).
In the following example we will use the overlap_integral_asymmetric method to project the MO coefficient matrix of \(\text{HeH}^{+}\) from the 6-31G basis set into the 6-311G basis set. For this we will use the MO coefficient matrix from the 6-31G basis set as the transformation matrix of this basis set into the MO basis.
# project MO coefficients to new basis
mo_coeffs_new = overlap_integral_asymmetric(ao_basis_new, ao_basis, transform_two=mo_coeffs.T)
print(f"Shape of MO coefficients: {mo_coeffs_new.shape}") # (nmo, nbasis)
print("MO coefficients projected in the 6-311G basis:")
print(mo_coeffs_new, end="\n\n")
print("Overlap matrix (S) of projected MOs (in the 6-311G basis):")
int1e_ovlp_mo_new = overlap_integral(ao_basis_new, mo_coeffs_new.T)
print(int1e_ovlp_mo_new, end="\n\n")
is_orthogonal = np.allclose(int1e_ovlp_mo_new, np.eye(int1e_ovlp_mo_new.shape[0])) # check if overlap is I
print(f"Projected MOs are orthogonal: {is_orthogonal}")
Shape of MO coefficients: (6, 4)
MO coefficients projected in the 6-311G basis:
[[ 0.4989715 0.52890845 0.61951534 0.2238233 ]
[ 0.68472982 0.6781266 0.00950138 0.16895291]
[ 0.54023005 0.64980473 -0.47179196 0.11732507]
[ 0.71577532 -0.26824295 0.01070959 -0.55172329]
[ 0.94958022 -0.21177303 -0.05871831 -0.03034387]
[ 0.8306327 0.14545388 -0.31205797 0.43115268]]
Overlap matrix (S) of projected MOs (in the 6-311G basis):
[[10.47059879 3.68860314 -0.96422155 1.2625921 ]
[ 3.68860314 2.62984775 -0.24828188 1.1829303 ]
[-0.96422155 -0.24828188 0.57472127 -0.13061703]
[ 1.2625921 1.1829303 -0.13061703 0.69815484]]
Projected MOs are orthogonal: False
# project MO coefficients to new basis
mo_coeffs_new = overlap_integral_asymmetric(ao_basis, ao_basis, transform_one=mo_coeffs.T)
print(f"Shape of MO coefficients: {mo_coeffs_new.shape}") # (nmo, nbasis)
print("MO coefficients projected in the 6-311G basis:")
print(mo_coeffs_new, end="\n\n")
print("Overlap matrix (S) of projected MOs (in the 6-311G basis):")
int1e_ovlp_mo_new = overlap_integral(ao_basis, mo_coeffs_new.T)
print(int1e_ovlp_mo_new, end="\n\n")
is_orthogonal = np.allclose(int1e_ovlp_mo_new, np.eye(int1e_ovlp_mo_new.shape[0])) # check if overlap is I
print(f"Projected coefficients are orthogonal: {is_orthogonal}")
Shape of MO coefficients: (4, 4)
MO coefficients projected in the 6-311G basis:
[[ 0.60290436 0.62310221 0.90324527 0.87056796]
[ 0.60630882 0.68647993 -0.2732068 0.07995231]
[ 0.47343187 -0.3462874 -0.01841626 -0.26148619]
[ 0.21156138 0.14343635 -0.3304041 0.40907705]]
Overlap matrix (S) of projected MOs (in the 6-311G basis):
[[2.29991242 1.53077726 0.3317762 1.46786878]
[1.53077726 1.4545096 0.39176399 1.27224729]
[0.3317762 0.39176399 0.49099634 0.46349437]
[1.46786878 1.27224729 0.46349437 1.23594644]]
Projected coefficients are orthogonal: False
1.3. Multipole Moment Integrals#
The moment_integral function computes the multipole moment integrals for two pairs of basis functions in AO basis, unless transform argument is provided. The center and order of multipole moments are specified by the moment_center and moment_orders arguments.
Example: Compute Dipole Moment#
The following example uses the moment_integral to compute molecular dipole moment. Given the atomic coordinates \(\{\mathbf{R}\}\) for \(N\) atoms, and the center of mass by \(\mathbf{R_c}\), the molecular dipole moment is given by:
from gbasis.integrals.moment import moment_integral
# extract the 1-RDM from the IOData object
rdm = mol.one_rdms['scf']
# set the orders of the moment integrals
p_ord = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
# compute center of nuclear charge of the molecule
mol_cm = mol.atcorenums @ mol.atcoords / np.sum(mol.atcorenums)
# calculate the nuclear moment of the molecule
mol_np_111 = mol.atcorenums @ (mol.atcoords - mol_cm)
print(f"{'Nuclear Dipole Moment of the molecule: ':>40s}{mol_np_111}")
# compute moment integrals of a set of molecular orbitals
int1e_mol_ep_111 = moment_integral(ao_basis, mol_cm, p_ord)
# compute electric moments of the molecule \mu_k = sum rdm_ij * p_ijk
mol_ep_111 = np.einsum("ij,ijk->k", rdm, int1e_mol_ep_111)
# gbasis electric moment of the molecule
print(f"{'Electric Dipole Moment of the molecule: ':>40s}{mol_ep_111}")
# calculate the total moment of the molecule
mol_tp_111 = mol_np_111 - mol_ep_111
print(f"{'Dipole Moment of the molecule: ':>40s}{mol_tp_111}")
print(f"{'Reference value: ':>40s}{mol.moments[(1,'c')]}")
Nuclear Dipole Moment of the molecule: [-1.23259516e-32 0.00000000e+00 0.00000000e+00]
Electric Dipole Moment of the molecule: [-9.81527656e-20 0.00000000e+00 5.60557428e-01]
Dipole Moment of the molecule: [ 9.81527656e-20 0.00000000e+00 -5.60557428e-01]
Reference value: [ 7.13633937e-18 0.00000000e+00 -5.60557427e-01]
1.4. Integrals over Differential Operator#
The gbasis.integrals module supports the computation of integrals over arbitrary differential operators.
Examples of these integrals are the kinetic energy integrals, and the momentum integrals.
1.5. Kinetic Energy Integrals#
The kinetic_energy_integral computes the kinetic energy integrals between pairs of basis functions in AO basis, unless transform argument is provided.
# From Gaussian16 Output File:
# Kinetic energy integral matrix in AO basis:
# 1 0.139568D+01
# 2 0.259735D+00 0.241917D+00
# 3 0.105657D+00 0.115030D+00 0.292565D+01
# 4 0.236735D+00 0.199566D+00 0.467227D+00 0.446946D+00
from gbasis.integrals.kinetic_energy import kinetic_energy_integral
# compute kinetic energy integrals in AO basis
int1e_kin_ao = kinetic_energy_integral(ao_basis)
print(f"Shape of kinetic energy integral matrix in AO basis: {int1e_kin_ao.shape}", end="\n\n")
print(f"Kinetic energy integral matrix in AO basis:")
print(int1e_kin_ao, end="\n\n")
# compute kinetic energy integrals in MO basis
int1e_kin_mo = kinetic_energy_integral(ao_basis, transform=mo_coeffs.T)
print(f"Shape of Kinetic Energy Integral matrix in MO basis: {int1e_kin_mo.shape}", end="\n\n")
print(f"Kinetic energy integral matrix in MO basis:")
print(int1e_kin_mo, end="\n\n")
Shape of kinetic energy integral matrix in AO basis: (4, 4)
Kinetic energy integral matrix in AO basis:
[[1.39567838 0.259735 0.10565728 0.23673475]
[0.259735 0.24191664 0.1150301 0.19956551]
[0.10565728 0.1150301 2.92565486 0.46722685]
[0.23673475 0.19956551 0.46722685 0.446946 ]]
Shape of Kinetic Energy Integral matrix in MO basis: (4, 4)
Kinetic energy integral matrix in MO basis:
[[ 1.46668936 -0.49993274 0.43649829 -1.52368766]
[-0.49993274 0.7947671 0.50202239 0.83209984]
[ 0.43649829 0.50202239 1.73387267 0.16632084]
[-1.52368766 0.83209984 0.16632084 3.39093136]]
1.6. Momentum Integrals#
The momentum_integral function computes the momentum integrals between pairs of basis functions in AO basis, unless transform argument is provided. This matrix corresponds to the unitary transformation from the position space to the momentum space representation. In the example below, the momentum vector for each AO and MO orbital is computed.
from gbasis.integrals.momentum import momentum_integral
# compute momentum integrals in AO basis
int1_p_ao = momentum_integral(ao_basis)
print(f"Shape of Momentum Integral matrix in AO basis: {int1_p_ao.shape}")
for i in range(len(ao_basis)):
print(f"Momentum for shell # {i+1}: {np.array2string(int1_p_ao[i][i], precision=2)}")
# compute momentum integrals in MO basis
int1_p_mo = momentum_integral(ao_basis)
round
print(f"\nShape of Momentum Integral matrix in MO basis: {int1_p_mo.shape}")
for i in range(len(ao_basis)):
print(f"Momentum for shell # {i+1}: {np.array2string(int1_p_mo[i][i], precision=2)}")
Shape of Momentum Integral matrix in AO basis: (4, 4, 3)
Momentum for shell # 1: [0.+0.0e+00j 0.+0.0e+00j 0.+1.6e-16j]
Momentum for shell # 2: [0.+0.00e+00j 0.+0.00e+00j 0.-3.58e-17j]
Momentum for shell # 3: [0.+0.00e+00j 0.+0.00e+00j 0.+1.65e-16j]
Momentum for shell # 4: [0.+0.00e+00j 0.+0.00e+00j 0.+3.31e-17j]
Shape of Momentum Integral matrix in MO basis: (4, 4, 3)
Momentum for shell # 1: [0.+0.0e+00j 0.+0.0e+00j 0.+1.6e-16j]
Momentum for shell # 2: [0.+0.00e+00j 0.+0.00e+00j 0.-3.58e-17j]
Momentum for shell # 3: [0.+0.00e+00j 0.+0.00e+00j 0.+1.65e-16j]
Momentum for shell # 4: [0.+0.00e+00j 0.+0.00e+00j 0.+3.31e-17j]
1.7. Angular Momentum Integrals#
The angular_momentum_integral function computes the angular momentum integral matrix for the AO and MO basis functions. The example below shows how to compute the angular momentum vector and norm.
from gbasis.integrals.angular_momentum import angular_momentum_integral
# compute angular momentum integrals in AO basis
int1_l_ao = angular_momentum_integral(ao_basis)
print("|{:^14}|{:^9}|{:^36}|".format("Contraction #", "Shell L", "Angular Momentum (x,y,z)"))
for i, ao_b in enumerate(ao_basis):
print(f"|{i+1:^14}|{ao_b.angmom:^9}| {np.array2string(int1_l_ao[i][i], precision=0)} |")
print()
# compute angular momentum integrals in MO basis
int1_l_mo = angular_momentum_integral(ao_basis, transform=mo_coeffs.T)
print("|{:^14}|{:^36}|".format("MO #", "Angular Momentum (x,y,z)"))
for i in range(len(ao_basis)):
print(f"|{i+1:^14}| {np.array2string(int1_l_mo[i][i], precision=0)} |")
|Contraction # | Shell L | Angular Momentum (x,y,z) |
| 1 | 0 | [0.+0.e+00j 0.-5.e-33j 0.+0.e+00j] |
| 2 | 0 | [0.+0.e+00j 0.+1.e-33j 0.+0.e+00j] |
| 3 | 0 | [0.+0.e+00j 0.-5.e-33j 0.+0.e+00j] |
| 4 | 0 | [0.+0.e+00j 0.-1.e-33j 0.+0.e+00j] |
| MO # | Angular Momentum (x,y,z) |
| 1 | [0.+0.e+00j 0.+5.e-18j 0.+0.e+00j] |
| 2 | [0.+0.e+00j 0.-2.e-17j 0.+0.e+00j] |
| 3 | [0.+0.e+00j 0.-2.e-19j 0.+0.e+00j] |
| 4 | [0.+0.e+00j 0.-1.e-17j 0.+0.e+00j] |
1.8. Nuclear Attraction Integrals#
The nuclear_electron_attraction function computes the nuclear attraction integrals to a set of nuclei of \(\{Z_A\}\) located at \(\{\mathbf{R}_{A}\}\) for pairs of AO and MO basis functions.
# From Gaussian16 Output File:
# Reference values:
# 1 2 3 4
# 1 0.300634D+01
# 2 0.159287D+01 0.168097D+01
# 3 0.118941D+01 0.124482D+01 0.549025D+01
# 4 0.149451D+01 0.159886D+01 0.244432D+01 0.235135D+01
from gbasis.integrals.nuclear_electron_attraction import nuclear_electron_attraction_integral
# compute nuclear-electron attraction integral in AO basis
int1e_nuc_ao = nuclear_electron_attraction_integral(ao_basis, mol.atcoords, mol.atcorenums)
print("Nuclear-Electron Attraction Integral of AOs (6-31G):")
print(int1e_nuc_ao, end="\n\n")
# compute nuclear-electron attraction integral in MO basis
int1e_nuc_mo = nuclear_electron_attraction_integral(
ao_basis, mol.atcoords, mol.atcorenums, transform=mo_coeffs.T
)
print("Nuclear-Electron Attraction Integral of MOs (6-31G):")
print(int1e_nuc_mo)
Nuclear-Electron Attraction Integral of AOs (6-31G):
[[-3.00634227 -1.59287383 -1.18941322 -1.49450521]
[-1.59287383 -1.68097365 -1.24481507 -1.59886117]
[-1.18941322 -1.24481507 -5.49024852 -2.44432492]
[-1.49450521 -1.59886117 -2.44432492 -2.35135264]]
Nuclear-Electron Attraction Integral of MOs (6-31G):
[[-4.11469726 0.64655654 -0.48797471 1.78211817]
[ 0.64655654 -2.01806794 -0.72857371 -1.09274852]
[-0.48797471 -0.72857371 -2.34608045 -0.12666144]
[ 1.78211817 -1.09274852 -0.12666144 -3.74760623]]
1.9. Point Charge Integrals#
The point_charge_integral function computes the point charge integral matrix for pairs of AO and MO basis functions. The represents the projection of the interaction between a (number) point charge(s) and an electron density onto a basis set.
In the example below, the nuclear attraction integrals (in the previous example) is reproduced by setting the point charge position and charge to the coordinates of the nucleus and their corresponding atomic number.
from gbasis.integrals.point_charge import point_charge_integral
# compute integral for interaction between these point charges and the set of atomic orbitals
int1e_pc_ao = point_charge_integral(ao_basis, mol.atcoords, mol.atcorenums)
print(f"Shape of Point Charge Integral matrix in AO basis: {int1e_pc_ao.shape}")
print("NOrbital, NOrbital, NAtom => Need to sum over NAtom", end="\n\n")
print("Nuclear-Electron Attraction Integral of AOs (6-31G):")
print(int1e_pc_ao.sum(axis=2), end="\n\n")
# compute integral for interaction between these point charges and the set of molecular orbitals
int1e_pc_mo = point_charge_integral(ao_basis, mol.atcoords, mol.atcorenums, transform=mo_coeffs.T)
print("Nuclear-Electron Attraction Integral of MOs (6-31G):")
print(int1e_pc_mo.sum(axis=2), end="\n\n")
Shape of Point Charge Integral matrix in AO basis: (4, 4, 2)
NOrbital, NOrbital, NAtom => Need to sum over NAtom
Nuclear-Electron Attraction Integral of AOs (6-31G):
[[-3.00634227 -1.59287383 -1.18941322 -1.49450521]
[-1.59287383 -1.68097365 -1.24481507 -1.59886117]
[-1.18941322 -1.24481507 -5.49024852 -2.44432492]
[-1.49450521 -1.59886117 -2.44432492 -2.35135264]]
Nuclear-Electron Attraction Integral of MOs (6-31G):
[[-4.11469726 0.64655654 -0.48797471 1.78211817]
[ 0.64655654 -2.01806794 -0.72857371 -1.09274852]
[-0.48797471 -0.72857371 -2.34608045 -0.12666144]
[ 1.78211817 -1.09274852 -0.12666144 -3.74760623]]
Example: Computing Core Hamiltonian Matrix#
The core Hamiltonian matrix is defined as:
where \(T_{\mu\nu}\) are the kinetic energy integral and \(V_{\mu\nu}\) are the nuclear attraction integral for basis functions \(\mu\) and \(\nu\). In the example below we will compute the core Hamiltonian matrix for the AO basis functions using the previously computed kinetic energy and nuclear attraction integral matrices.
# From Gaussian16 Output File (line 277 of HeHp.log):
# Core Hamiltonian matrix (H_core) in AO basis:
# 1 -0.161066D+01
# 2 -0.133314D+01 -0.143906D+01
# 3 -0.108376D+01 -0.112978D+01 -0.256459D+01
# 4 -0.125777D+01 -0.139930D+01 -0.197710D+01 -0.190441D+01
int1e_core = int1e_nuc_ao + int1e_kin_ao
print(f"Core Hamiltonian matrix (H_core) in AO basis:")
print(int1e_core, end="\n\n")
Core Hamiltonian matrix (H_core) in AO basis:
[[-1.6106639 -1.33313883 -1.08375595 -1.25777046]
[-1.33313883 -1.43905701 -1.12978497 -1.39929566]
[-1.08375595 -1.12978497 -2.56459366 -1.97709807]
[-1.25777046 -1.39929566 -1.97709807 -1.90440664]]
2. Two-Electron Repulsion Integrals#
The electron_repulsion function compute the electron-electron repulsion integrals in AO or MO basis for a pir of basis functions. In the example below, these integrals will be used to compute the Coulomb J energy.
from gbasis.integrals.electron_repulsion import electron_repulsion_integral
# compute electron-electron repulsion integral in AO basis
int2e_ao = electron_repulsion_integral(ao_basis, notation="chemist")
print(f"Shape of Electron-Electron Repulsion Integrals: {int2e_ao.shape} (#AO, #AO, #AO, #AO)")
# compute electron-electron repulsion integral in MO basis
int2e_mo = electron_repulsion_integral(ao_basis, notation="chemist", transform=mo_coeffs.T)
print(f"Shape of Electron-Electron Repulsion Integrals: {int2e_mo.shape} (#MO, #MO, #MO, #MO)")
Shape of Electron-Electron Repulsion Integrals: (4, 4, 4, 4) (#AO, #AO, #AO, #AO)
Shape of Electron-Electron Repulsion Integrals: (4, 4, 4, 4) (#MO, #MO, #MO, #MO)