# Gbasis publication examples#

This notebook provides the examples of the use of `Gbasis`

provided in the publication article.

## Install dependencies and download data#

If this this notebook is used in google colab it requires several dependencies and data. The next cell will install the required dependencies if they are not already installed. It will also download the data required for the examples if it is not already downloaded. Uncomment the cell below and execute it.

```
# ! python -m pip install numpy
# ! python -m pip install matplotlib
# ! python -m pip install scipy
# ! python -m pip install urllib3
# ! python -m pip install git+https://github.com/theochem/grid.git
# ! python -m pip install git+https://github.com/theochem/iodata.git
# ! python -m pip install git+https://github.com/theochem/gbasis.git
# import os
# from urllib.request import urlretrieve
# # download the required data files
# file_path_1 = "hydrogen_def2-svp.1.gbs"
# if not os.path.isfile(file_path_1):
# url = "https://github.com/theochem/gbasis/blob/master/notebooks/tutorial/hydrogen_def2-svp.1.gbs?raw=true"
# urlretrieve(url, file_path_1)
# file_path_2 = "C2H4_hf_ccpvdz.fchk"
# if not os.path.isfile(file_path_2):
# url = "https://github.com/theochem/gbasis/blob/master/notebooks/tutorial/C2H4_hf_ccpvdz.fchk?raw=true"
# urlretrieve(url, file_path_2)
# file_path_3 = "C2H4_hf_ccpvtz.fchk"
# if not os.path.isfile(file_path_3):
# url = "https://github.com/theochem/gbasis/blob/master/notebooks/tutorial/C2H4_hf_ccpvtz.fchk?raw=true"
# urlretrieve(url, file_path_3)
```

## A. Building Basis Functions#

Gbasis supports loading basis functions through two different ways:

Basis functions from basis set text files. Currently it supports Gaussian94 and NewChem type files. The

`coord_type`

argument allows the user to specify Cartesian (“c” or “cartesian”), spherical (“p” or “spherical”), or mixed coordinate (as a list) for each contraction.

```
import numpy as np
from gbasis.parsers import parse_gbs, make_contractions
# Define atomic symbols and coordinates (i.e., basis function centers)
atoms = ["H", "H"]
atcoords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
# Obtain basis functions from the basis set files
basis_dict = parse_gbs("hydrogen_def2-svp.1.gbs")
basis = make_contractions(
basis_dict, atoms, atcoords, coord_types="c")
# To obtain the total number of AOs we compute the cartesian components for each angular momentum
total_ao = 0
print(f"Number of generalized shells: {len(basis)}") # Output 6
for shell in basis:
total_ao += shell.angmom_components_cart.shape[0]
print("Total number of AOs: ", total_ao) # output 10
```

```
Number of generalized shells: 6
Total number of AOs: 10
```

Quantum chemistry calculations through interfaces with IOData and PySCF packages. IOData provides unified access to various computational chemistry file formats (e.g. formatted checkpoint file, molden, and wfn/wfx wavefunction files). Basis set data stored in

`iodata`

object is converted to a list of`GeneralizedContractionShell`

objects, as exemplified below.

```
from iodata import load_one
from gbasis.wrappers import from_iodata
# Load ethylene fchk
mol = load_one("C2H4_hf_ccpvdz.fchk")
basis = from_iodata(mol)
print(f"Number of generalized shells: {len(basis)}") # output 24
# To obtain the total number of AOs we check for each shell its angular momentum and coordinate type
total_ao = 0
for shell in basis:
if shell.coord_type == "cartesian":
total_ao += shell.angmom_components_cart.shape[0]
elif shell.coord_type == "spherical":
total_ao += len(shell.angmom_components_sph)
print("Total number of AOs: ", total_ao) # output 48
```

```
Number of generalized shells: 24
Total number of AOs: 48
```

Each shell is an object of the base class `GeneralizedContractionShell`

. This object contains all the information to perform calculations for the different quantities available in the package. For the first shell in the basis we can access them:

```
# For shell 1
shell = basis[0]
print(f"Generalized contraction Shell = 1")
print("Coordinates (Bohr):")
print(shell.coord)
print("Shape:",shell.coord.shape, "--> x,y,z coordinates")
print("Angular momentum: ", shell.angmom)
print("Exponents primitives:")
print(shell.exps)
print("Shape:", shell.exps.shape, "--> (K,) K=# primitives")
print("Contraction coefficients")
print(shell.coeffs)
print("Shape:", shell.coeffs.shape, "--> (K, M) \
K=# of primitives, M=# number segmented contraction shells")
print("Normalization constants:")
print(shell.norm_cont)
print("Shape:", shell.norm_cont.shape, "--> (M, L), \
M=# number segmented contraction shells, \
L=# different angular momentum components")
```

```
Generalized contraction Shell = 1
Coordinates (Bohr):
[-1.41331444e-15 -8.50766406e-17 1.24804461e+00]
Shape: (3,) --> x,y,z coordinates
Angular momentum: 0
Exponents primitives:
[6.665e+03 1.000e+03 2.280e+02 6.471e+01 2.106e+01 7.495e+00 2.797e+00]
Shape: (7,) --> (K,) K=# primitives
Contraction coefficients
[[0.00069352]
[0.0053415 ]
[0.02713667]
[0.10199239]
[0.27550864]
[0.45108643]
[0.28756574]]
Shape: (7, 1) --> (K, M) K=# of primitives, M=# number segmented contraction shells
Normalization constants:
[[1.]]
Shape: (1, 1) --> (M, L), M=# number segmented contraction shells, L=# different angular momentum components
```

We can retrieve information from the molecular system as well as the wavefunction (i.e molecular orbtials) stored inside `iodata`

object in the `mol`

variables. With molecular orbitals coefficients and occupations it is possible to construct the density matrix in terms of atomic orbitals. This will be used to compute various quantities.

```
# Get atomic numbers, atomic charges,
# atomic masses and Cartesian coordinates
atnums, atcharges = mol.atnums, mol.atcorenums
atcoords, atmasses = mol.atcoords, mol.atmasses
# Get molecular orbitals coefficients and occupations
mo_coeffs, mo_occs = mol.mo.coeffs, mol.mo.occs
# Calculate density matrix
dm = np.dot(mo_coeffs * mo_occs, mo_coeffs.T)
```

## B. Computing Integrals#

The `gbasis.integrals`

module supports various 1- and 2-electron integrals as well as an interface to `Libcint`

library.

## 1. One-Electron Integrals#

### 1.1. Overlap Integrals#

The overlap integral is computed by `overlap_integral`

. This returns a matrix of size (AOs, AOs). The molecular overlap can also be generated if the transformation matrix (i.e molecular coefficients) is passed through the `transform`

attribute. To speed up overlap integral calculation a screening based on the exponents and the distance between
their centers of the basis functions. To use this functionality the molecule neeeds to be loaded and the argument `overlap`

set to `True`

(default is `False`

)

```
import numpy as np
from gbasis.integrals.overlap import overlap_integral
# compute overlap integrals in AO and MO basis
olp_ao = overlap_integral(basis)
olp_mo = overlap_integral(basis, transform=mo_coeffs.T)
# check whether overlap integrals are orthonormal
print("Is AO Overlap Normalized?", np.allclose(np.diag(olp_ao),np.ones(total_ao)))
print("Is AO Overlap Orthogonal?", np.allclose(olp_ao, np.eye(total_ao)))
print("Is MO Overlap Orthonormal?", np.allclose(olp_mo, np.eye(total_ao), atol=1e-7))
# compute overlap integrals in AO and MO basis with screening
olp_ao_screening = overlap_integral(basis, tol_screen=1.0e-5)
olp_mo_screening = overlap_integral(basis, tol_screen=1.0e-5, transform=mo_coeffs.T)
print("Is Overlap equal to Overlap screening?", np.allclose(olp_ao, olp_ao_screening, atol=1e-4, rtol=1e-4))
```

```
Is AO Overlap Normalized? True
Is AO Overlap Orthogonal? False
Is MO Overlap Orthonormal? True
Is Overlap equal to Overlap screening? 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 `CC-pVDZ`

and `CC-pTDZ`

basis sets:

```
from gbasis.parsers import parse_gbs, make_contractions
from gbasis.integrals.overlap_asymm import overlap_integral_asymmetric
# load ethylene with basis set cc-pTDZ
mol_extended_basis = load_one("C2H4_hf_ccpvtz.fchk")
basis_extended_basis = from_iodata(mol_extended_basis)
print(f"Number of shells in CC-pVDZ basis: {len(basis_extended_basis)}")
print(f"Number of shells in CC-pTDZ basis: {len(basis)}", end="\n\n")
# compute overlap of two different basis sets
olp_asym = overlap_integral_asymmetric(basis_extended_basis, basis)
print(f"Shape of overlap matrix: {olp_asym.shape}")
```

```
Number of shells in CC-pVDZ basis: 44
Number of shells in CC-pTDZ basis: 24
Shape of overlap matrix: (116, 48)
```

### 1.3 Integral over arbitrary differential operator#

The `gbasis.integrals`

module supports the computation of integrals over arbitrary differential operators.

### 1.3.1 Kinetic energy#

The `kinetic_energy_integral`

computes the kinetic energy integrals between pairs of basis functions in AO basis, unless `transform`

argument is provided.

```
from gbasis.integrals.kinetic_energy import kinetic_energy_integral
# compute kinetic energy integrals in AO basis
k_int1e = kinetic_energy_integral(basis)
print("Shape kinetic energy integral: ", k_int1e.shape, "(#AO, #AO)")
kin_e = np.trace(dm.dot(k_int1e))
print("Kinetic energy (Hartree):", kin_e)
```

```
Shape kinetic energy integral: (48, 48) (#AO, #AO)
Kinetic energy (Hartree): 77.9285467708042
```

### 1.4 Nuclear electron attraction integral#

The `nuclear_electron_attraction`

function computes the nuclear attraction integrals to a set of nuclei of \(\{Z_C\}\) located at \(\{\mathbf{R}_{C}\}\) for pairs of AO or MO basis functions. Its computation is build upon the The `point_charge_integral`

.

```
from gbasis.integrals.nuclear_electron_attraction import \
nuclear_electron_attraction_integral
# compute nuclear-electron attraction integrals in AO basis
nuc_int1e = nuclear_electron_attraction_integral(
basis, atcoords, atnums)
print("Shape Nuclear-electron integral: ", nuc_int1e.shape, "(#AO, #AO)")
ne_e = np.trace(dm.dot(nuc_int1e))
print("Nuclear-electron energy (Hartree):", ne_e)
```

```
Shape Nuclear-electron integral: (48, 48) (#AO, #AO)
Nuclear-electron energy (Hartree): -248.2770425530744
```

## 2. Two-Electron Repulsion Integrals#

The `electron_repulsion`

function compute the electron-electron repulsion integrals in AO or MO basis for a pair of basis functions. This integrals can be used to compute Coulomb(J) and Exchange(K)

```
from gbasis.integrals.electron_repulsion import electron_repulsion_integral
#Compute e-e repulsion integral in MO basis, shape=(#MO, #MO, #MO, #MO)
int2e_mo = electron_repulsion_integral(basis, transform=mo_coeffs.T, notation='chemist')
print('Shape e-e repulsion integrals:',int2e_mo.shape,'(#MO, #MO, #MO, #MO)')
j_coul = 0
k_ex = 0
# Mask only occupied Molecular orbitals
occ_mo = mo_occs[mo_occs > 0].shape[0]
for i in range(occ_mo):
for j in range(occ_mo):
j_coul += 2 * int2e_mo[i,i,j,j]
k_ex += int2e_mo[i,j,i,j]
print("Coulomb energy (Hartree):", j_coul)
print("Exchange energy (Hartree):", k_ex)
```

```
Shape e-e repulsion integrals: (48, 48, 48, 48) (#MO, #MO, #MO, #MO)
Coulomb energy (Hartree): 70.53719756188356
Exchange energy (Hartree): 11.756060443132313
```

```
# Compute Nucleus-Nucleus repulsion
rab = np.triu(np.linalg.norm(atcoords[:, None]- atcoords, axis=-1))
at_charges = np.triu(atnums[:, None] * atnums)[np.where(rab > 0)]
nn_e = np.sum(at_charges / rab[rab > 0])
# Combine all terms to obtain total energy at Restricted HF level
e = nn_e + ne_e + kin_e + j_coul - k_ex
print(f"Total energy - Gbasis (Hartree): , {e: 5.11f}")
# SCF energy is stored in iodata object
print("Total energy - Gaussian16 (Hartre): ", mol.energy ) #-78.0401652
```

```
Total energy - Gbasis (Hartree): , -78.04016522086
Total energy - Gaussian16 (Hartre): -78.04016529602737
```

## 3. Libcint interface#

Gbasis includes a python interface to `Libcint`

library. To access the interface a `CBasis`

objects needs to be initialized, which needs as arguments `basis`

, list of atom elements and coordinates in atomic units (Bohr). Different from native functionality in Gbasis, Libcint does not work with mix type basis, meaning that all shells passed have to be in either cartesian or spherical coordinates. To showcase its functionality below we compute again the same energy terms. Note: To execute the code Libcint must be installed.

### IMPORTANT:#

This block deals with Libcint installation. This needs to be uncommented if you are executing this notebook on Google colab and want to use Libcint interface. Do not uncoment if you are using this notebook locally, Libcint should be installed using bash in a terminal following `README.md`

instructions. If Python version is not 3.10 this will most likely break.

```
# %%shell
# git clone https://github.com/theochem/gbasis.git
# cp -r gbasis/tools /usr/local/lib/python3.10/dist-packages/gbasis/
# sudo apt-get install -y sbcl
# cd /usr/local/lib/python3.10/dist-packages/gbasis/; ./tools/install_libcint.sh
# mv /usr/local/lib/python3.10/dist-packages/gbasis/gbasis/integrals/include/ /usr/local/lib/python3.10/dist-packages/gbasis/integrals
# mv /usr/local/lib/python3.10/dist-packages/gbasis/gbasis/integrals/lib/ /usr/local/lib/python3.10/dist-packages/gbasis/integrals/
```

```
from gbasis.integrals.libcint import CBasis
# Define atomic elements matching atom numbers
atoms_ethylene = ["C", "H", "H", "C", "H", "H"]
cbasis = CBasis(basis, atoms_ethylene, atcoords, coord_type="spherical")
kin_int1e_c= cbasis.kinetic_energy_integral()
kin_e_c = np.trace(dm.dot(kin_int1e_c))
nuc_int1e_c = cbasis.nuclear_attraction_integral()
ne_e_c = np.trace(dm.dot(nuc_int1e_c))
int2e_mo_c = cbasis.electron_repulsion_integral(transform=mo_coeffs.T, notation='chemist')
j_coul_c = 0
k_ex_c = 0
# Mask only occupied Molecular orbitals
occ_mo = mo_occs[mo_occs > 0].shape[0]
occ_mo = mol.mo.occs[mol.mo.occs > 0].shape[0]
for i in range(occ_mo):
for j in range(occ_mo):
j_coul_c += 2 * int2e_mo_c[i,i,j,j]
k_ex_c += int2e_mo_c[i,j,i,j]
print("Kinetic energy (CBasis/Hartree)):", kin_e_c)
print("Nuclear-electron energy (CBasis/Hartree)):", ne_e_c)
print("Coulomb energy (CBasis/Hartree)):", j_coul_c)
print("Exchange energy (CBasis/Hartree)):", k_ex_c)
# Combine all terms to obtain total energy at Restricted HF level
e_cbasis = nn_e + ne_e_c + kin_e_c + j_coul_c - k_ex_c
print(f"Total energy - Gbasis (Hartree): , {e: 5.11f}")
print(f"Total energy - Cbasis (Hartree): , {e_cbasis: 5.11f}")
# SCF energy is stored in iodata object
print("Total energy - Gaussian16 (Hartree): ", mol.energy ) #-78.0401652
```

```
Kinetic energy (CBasis/Hartree)): 77.92854683833846
Nuclear-electron energy (CBasis/Hartree)): -248.27704265438183
Coulomb energy (CBasis/Hartree)): 70.5371975684723
Exchange energy (CBasis/Hartree)): 11.756060454118147
Total energy - Gbasis (Hartree): , -78.04016522086
Total energy - Cbasis (Hartree): , -78.04016525903
Total energy - Gaussian16 (Hartree): -78.04016529602737
```

## C. Evaluations#

## 1. Basis functions#

The `gbasis.evals`

module supports evaluating functions that are expanded in
Gaussian basis set on a set of points.For the following examples, we use the `Grid`

library to
generate the Becke-Lebedev molecular grid:

```
from grid.molgrid import MolGrid
grid = MolGrid.from_preset(
atnums=mol.atnums, atcoords=atcoords, preset="fine")
print("Number of grid points = ", grid.size)
```

```
Number of grid points = 16796
```

The following example show how to evaluate atomic and molecular orbitals and their derivatives

```
from gbasis.evals.eval import evaluate_basis
from gbasis.evals.eval_deriv import evaluate_deriv_basis
# Evaluate the MOs on the grid points
basis_mo = evaluate_basis(basis, grid.points, transform=mo_coeffs.T)
print("Shape MOs evaluated in a Molecular grid: ")
print(basis_mo.shape, "(#MOs, #Grid points)")
# Integration MO
# Mask only occupied Molecular orbitals
occ_mo = mo_occs[mo_occs > 0].shape[0]
for i in range(occ_mo): # Only occupied MOs
eval_mo = basis_mo[i] * basis_mo[i]
print(f"MO {i+1} integrated squared: ", grid.integrate(eval_mo))
# evaluate 6th derivative w.r.t. y for MOs on the grid points
deriv_6_mo = evaluate_deriv_basis(basis, grid.points,
orders=np.array([0,6,0]), transform=mo_coeffs.T)
```

```
Shape MOs evaluated in a Molecular grid:
(48, 16796) (#MOs, #Grid points)
MO 1 integrated squared: 1.0000027687589634
MO 2 integrated squared: 1.0000027654426193
MO 3 integrated squared: 0.9999978011521944
MO 4 integrated squared: 0.9999860542308123
MO 5 integrated squared: 0.9999893708000488
MO 6 integrated squared: 0.9999931468746739
MO 7 integrated squared: 0.99997589380213
MO 8 integrated squared: 0.9999976856899412
```

## 4. Electrostatic potential#

An high-level function for computing the molecular electrostatic potential is also provided. It is more suitable to use a cubic grid for this computation, which is typically used for visualizing molecular surface properties.

```
from grid import UniformGrid
from gbasis.evals.electrostatic_potential import electrostatic_potential
# construct a cubic grid
cube_grid = UniformGrid.from_molecule(atnums,
atcoords,
spacing=0.75,
extension=5.0,
rotate=True)
# compute electrostatic potential on the grid point
esp = electrostatic_potential(basis=basis,
one_density_matrix=dm,
points=cube_grid.points,
nuclear_coords=atcoords,
nuclear_charges=atcharges)
print("Shape electrostatic potential output: ", esp.shape, "#Grid points")
# generate a cube
cube_grid.generate_cube(fname="ethylene.cube",
data=esp,
atcoords=atcoords,
atnums=atnums,
pseudo_numbers=atcharges)
```

```
Shape electrostatic potential output: (5040,) #Grid points
```