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:

1. 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

2. 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

2. Density Related properties

Gbasis supports evaluating the electron density and arbitrary-order derivatives

from gbasis.evals.density import evaluate_density

# compute electronic density
dens = evaluate_density(dm, basis, grid.points)
# integrate density and compare to the number of electrons in the system
print("Total number of electrons: ", grid.integrate(dens)) # Total number of electrons ethylene = 16

Total number of electrons:  15.99989097350276

In the special case of the gradient, Hessian, and the Laplacian of the electron density, Gbasis has high-level functions implementing the calculation and combination of the necessary derivatives orders. Up to second order, the user has the option to choose between the arbitrary-order-derivative implementation, or the explicitly coded 1st and 2nd derivatives.

from gbasis.evals.density import evaluate_density_laplacian

# compute Laplacian of electron density
# general arbitrary-order-derivative implementation
lap_gen = evaluate_density_laplacian(dm, basis, grid.points, deriv_type='general')
print(grid.integrate(lap_gen)) # Expected value: 0
# Direct implementation of 1s and 2nd order derivatives
lap_dir = evaluate_density_laplacian(dm, basis, grid.points, deriv_type='direct')
print(grid.integrate(lap_dir)) # Expected value: 0
print("Are General and Direct implementation the same?", 
      np.allclose(lap_gen, lap_dir))

-0.0049650299404096635
-0.004965029940581106
Are General and Direct implementation the same? True

3. Density Matrix Related properties

GBasis also supports the evaluation of the one-electron reduced density matrix and arbitrary order derivatives on a set of grid points. Supported related properties are the stress tensor, the Ehrenfest force and the Ehrenfest Hessian. Below we show an example on how to use the Ehrenfest force. Other properties are accessed similarly with their corresponding high level function.

from gbasis.evals.stress_tensor import evaluate_ehrenfest_force

#center of mass ethylene
com = (np.mean(atmasses[:,None] * atcoords, axis=0) / np.sum(atmasses))[None,:]
print(f"Center of mass of ethylene: {com}", end="\n\n")
alpha, beta = 1, 0
print(f"Ehrenfest force (alpha = {alpha} and beta = {beta}):")
ef_com = evaluate_ehrenfest_force(dm, 
                                  basis, 
                                  com, 
                                  alpha=alpha, 
                                  beta=beta)
print(f"Shape of a Ehrenfest force at one point: {ef_com[0].shape}")
print(f"Ehrenfest force at center of mass:\n{ef_com}")
print(f"Ehrenfest force at center of mass is 0: {np.allclose(ef_com, 0)}")
print(f"(should be zero, it is Lorentz point)", end="\n\n")

Center of mass of ethylene: [[ 5.03790449e-17  5.93301676e-18 -5.93301676e-18]]

Ehrenfest force (alpha = 1 and beta = 0):
Shape of a Ehrenfest force at one point: (3,)
Ehrenfest force at center of mass:
[[-8.59906292e-17 -3.88840561e-17  8.51183409e-18]]
Ehrenfest force at center of mass is 0: True
(should be zero, it is Lorentz point)

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