Calculating integrals over GTOs#
# Force the local gqcpy to be imported
import sys
sys.path.insert(0, '../../build/gqcpy/')
import gqcpy
import numpy as np
np.set_printoptions(precision=8, linewidth=120)
In this example, we’ll go over the high- and low-level machinery that GQCP offers in order to calculate integrals over Cartesian GTOs. We assume that you’re familiar with the mathematical concepts of a scalar basis, a shell, a basis function, and a primitive. If you’re not, here’s a very succinct overview:
Spin-orbitals are expanded in an underlying scalar basis.
In order to compactify scalar bases, we use shells that group primitives according to their angular momentum.
Therefore, in every shell, a set of basis functions (of the same angular momentum) is implicitly defined.
Every basis function is defined as a contraction (i.e. a linear combination, where the coefficients are called contraction coefficients) of primitives.
Here, the primitives are Cartesian GTOs.
We’ll start off by setting up a small molecular system and a scalar basis. Since we don’t want to manually read in a basis set, we’ll use RSpinOrbitalBasis’s functionality to provide us with a scalar basis.
nuclei = [gqcpy.Nucleus(1, 0.0, 0.0, 0.0), gqcpy.Nucleus(1, 0.0, 0.0, 1.0)]
molecule = gqcpy.Molecule(nuclei)
spin_orbital_basis = gqcpy.RSpinOrbitalBasis_d(molecule, "STO-3G")
scalar_basis = spin_orbital_basis.scalarBasis()
shell_set = scalar_basis.shellSet() # A shell set is just a collection of shells.
Primitive engines#
In GQCP, we define an engine to be a computational object that is able to calculate integrals over shells, while a primitive engine is defined to be a computational entity that can calculate integrals over primitives.
In this example, we’re taking the expansion of contracted GTOs in terms of their primitives for granted (using the implementations provided by the combination of FunctionalOneElectronIntegralEngine and IntegralCalculator), and we’ll be focusing on calculating integrals over primitives.
Overlap integrals#
Let’s get straight into it. Using the McMurchie-Davidson integral scheme, we can calculate the overlap integral over two primitives as follows.
def overlap_function_1D(a, K, i, b, L, j):
# Negative Cartesian exponents should be ignored: the correct value for the corresponding integral is 0.
if (i < 0) or (j < 0):
return 0.0
# Use the McMurchie-Davidson recursion to calculate the overlap integral.
p = a + b
E = gqcpy.McMurchieDavidsonCoefficient(K, a, L, b)
return np.power(np.pi / p, 0.5) * E(i, j, 0)
def overlap_function(left, right):
primitive_integral = 1.0
# The overlap integral is separable in its three Cartesian components.
for direction in [gqcpy.CartesianDirection.x, gqcpy.CartesianDirection.y, gqcpy.CartesianDirection.z]:
i = left.cartesianExponents().value(direction)
j = right.cartesianExponents().value(direction)
a = left.gaussianExponent()
b = right.gaussianExponent()
K = left.center()[direction]
L = right.center()[direction]
primitive_integral *= overlap_function_1D(a, K, i, b, L, j)
return primitive_integral
We’ll wrap this function in to a FunctionalOneElectronPrimitiveIntegralEngine, and then supply it as the primitive engine that should be used in a FunctionalOneElectronIntegralEngine. We’re doing this in order to use GQCP’s internal handling of the shells and contractions through the IntegralCalculator.calculate call.
primitive_overlap_engine = gqcpy.FunctionalOneElectronPrimitiveIntegralEngine_d(overlap_function)
overlap_engine = gqcpy.FunctionalOneElectronIntegralEngine_d(primitive_overlap_engine)
S = gqcpy.IntegralCalculator.calculate(overlap_engine, shell_set, shell_set)
print(S)
[[0.99999999 0.79658829]
[0.79658829 0.99999999]]
We can verify our results by letting the spin-orbital basis quantize the overlap operator.
S_ref = spin_orbital_basis.quantize(gqcpy.OverlapOperator()).parameters()
print(S_ref)
[[0.99999999 0.79658829]
[0.79658829 0.99999999]]
Kinetic integrals#
def kinetic_function_1D(a, K, i, b, L, j):
# The kinetic 1D integral is a sum of three 1D overlap integrals.
return -2.0 * np.power(b, 2) * overlap_function_1D(a, K, i, b, L, j + 2) + \
b * (2 * j + 1) * overlap_function_1D(a, K, i, b, L, j) - \
0.5 * j * (j - 1) * overlap_function_1D(a, K, i, b, L, j - 2)
def kinetic_function(left, right):
# Prepare some variables
i = left.cartesianExponents().value(gqcpy.CartesianDirection.x)
k = left.cartesianExponents().value(gqcpy.CartesianDirection.y)
m = left.cartesianExponents().value(gqcpy.CartesianDirection.z)
j = right.cartesianExponents().value(gqcpy.CartesianDirection.x)
l = right.cartesianExponents().value(gqcpy.CartesianDirection.y)
n = right.cartesianExponents().value(gqcpy.CartesianDirection.z)
a = left.gaussianExponent()
b = right.gaussianExponent()
K_x = left.center()[gqcpy.CartesianDirection.x]
K_y = left.center()[gqcpy.CartesianDirection.y]
K_z = left.center()[gqcpy.CartesianDirection.z]
L_x = right.center()[gqcpy.CartesianDirection.x]
L_y = right.center()[gqcpy.CartesianDirection.y]
L_z = right.center()[gqcpy.CartesianDirection.z]
# The 3D kinetic energy integral is a sum of three contributions (dx^2, dy^2, dz^2).
return kinetic_function_1D(a, K_x, i, b, L_x, j) * \
overlap_function_1D(a, K_y, k, b, L_y, l) * \
overlap_function_1D(a, K_z, m, b, L_z, n) + \
overlap_function_1D(a, K_x, i, b, L_x, j) * \
kinetic_function_1D(a, K_y, k, b, L_y, l) * \
overlap_function_1D(a, K_z, m, b, L_z, n) + \
overlap_function_1D(a, K_x, i, b, L_x, j) * \
overlap_function_1D(a, K_y, k, b, L_y, l) * \
kinetic_function_1D(a, K_z, m, b, L_z, n);
Same thing for the kinetic integrals. We wrap the kinetic_function in a functional primitive engine and proceed analoguously.
primitive_kinetic_engine = gqcpy.FunctionalOneElectronPrimitiveIntegralEngine_d(kinetic_function)
kinetic_engine = gqcpy.FunctionalOneElectronIntegralEngine_d(primitive_kinetic_engine)
T = gqcpy.IntegralCalculator.calculate(kinetic_engine, shell_set, shell_set)
print(T)
[[0.76003188 0.38325367]
[0.38325367 0.76003188]]
As a comparison, here are the GQCP integrals.
T_ref = spin_orbital_basis.quantize(gqcpy.KineticOperator()).parameters()
print(T_ref)
[[0.76003188 0.38325367]
[0.38325367 0.76003188]]
Nuclear attraction integrals#
A point of caution here. The type of functions that FunctionalOneElectronPrimitiveIntegralEngine_d accepts only has two CartesianGTO arguments. However, the nuclear attraction function needs to know which nucleic contributions have to be calculated. Since this is just a tutorial/example, we have opted to declare nuclei as a global variable.
def nuclear_attraction_function(left, right):
# Prepare some variables.
i = left.cartesianExponents().value(gqcpy.CartesianDirection.x)
k = left.cartesianExponents().value(gqcpy.CartesianDirection.y)
m = left.cartesianExponents().value(gqcpy.CartesianDirection.z)
j = right.cartesianExponents().value(gqcpy.CartesianDirection.x)
l = right.cartesianExponents().value(gqcpy.CartesianDirection.y)
n = right.cartesianExponents().value(gqcpy.CartesianDirection.z)
a = left.gaussianExponent()
b = right.gaussianExponent()
K_x = left.center()[gqcpy.CartesianDirection.x]
K_y = left.center()[gqcpy.CartesianDirection.y]
K_z = left.center()[gqcpy.CartesianDirection.z]
L_x = right.center()[gqcpy.CartesianDirection.x]
L_y = right.center()[gqcpy.CartesianDirection.y]
L_z = right.center()[gqcpy.CartesianDirection.z]
# Prepare the McMurchie-Davidson coefficients.
E_x = gqcpy.McMurchieDavidsonCoefficient(K_x, a, L_x, b)
E_y = gqcpy.McMurchieDavidsonCoefficient(K_y, a, L_y, b)
E_z = gqcpy.McMurchieDavidsonCoefficient(K_z, a, L_z, b)
p = a + b;
P = np.array([E_x.centerOfMass(), E_y.centerOfMass(), E_z.centerOfMass()])
# Calculate the contributions from every nuclear center.
total_integral = 0.0
for nucleus in nuclei:
integral = 0.0
C = nucleus.position()
R = gqcpy.HermiteCoulombIntegral(p, P, C)
for t in range(0, i + j + 1):
for u in range(0, k + l + 1):
for v in range(0, m + n + 1):
# Add the contribution to the integral. The prefactor will be applied at the end.
integral += E_x(i, j, t) * E_y(k, l, u) * E_z(m, n, v) * R(0, t, u, v)
total_integral += (-nucleus.charge()) * integral
return 2 * np.pi / p * total_integral
You’re getting the hang of this, right? Wrap the nuclear_attraction_function in the functional primitive engine and then provide that as input for the one-electron engine.
primitive_nuclear_engine = gqcpy.FunctionalOneElectronPrimitiveIntegralEngine_d(nuclear_attraction_function)
nuclear_engine = gqcpy.FunctionalOneElectronIntegralEngine_d(primitive_nuclear_engine)
V = gqcpy.IntegralCalculator.calculate(nuclear_engine, shell_set, shell_set)
print(V)
[[-2.03852055 -1.60241665]
[-1.60241665 -2.03852055]]
V_ref = spin_orbital_basis.quantize(gqcpy.NuclearAttractionOperator(molecule.nuclearFramework())).parameters()
print(V_ref)
[[-2.03852055 -1.60241665]
[-1.60241665 -2.03852055]]
Coulomb repulsion integrals#
The two-electron Coulomb repulsion integrals aren’t very different. They just require a little more code since the integrals are calculated over four basis functions instead of the previous two.
def coulomb_repulsion_function(left1, left2, right1, right2):
# Prepare some variables. Those with an extra underscore represent the 'primed' indices in the notes.
i = left1.cartesianExponents().value(gqcpy.CartesianDirection.x)
k = left1.cartesianExponents().value(gqcpy.CartesianDirection.y)
m = left1.cartesianExponents().value(gqcpy.CartesianDirection.z)
j = left2.cartesianExponents().value(gqcpy.CartesianDirection.x)
l = left2.cartesianExponents().value(gqcpy.CartesianDirection.y)
n = left2.cartesianExponents().value(gqcpy.CartesianDirection.z)
i_ = right1.cartesianExponents().value(gqcpy.CartesianDirection.x)
k_ = right1.cartesianExponents().value(gqcpy.CartesianDirection.y)
m_ = right1.cartesianExponents().value(gqcpy.CartesianDirection.z)
j_ = right2.cartesianExponents().value(gqcpy.CartesianDirection.x)
l_ = right2.cartesianExponents().value(gqcpy.CartesianDirection.y)
n_ = right2.cartesianExponents().value(gqcpy.CartesianDirection.z)
a = left1.gaussianExponent()
b = left2.gaussianExponent()
c = right1.gaussianExponent()
d = right2.gaussianExponent()
K_x = left1.center()[gqcpy.CartesianDirection.x]
K_y = left1.center()[gqcpy.CartesianDirection.y]
K_z = left1.center()[gqcpy.CartesianDirection.z]
L_x = left2.center()[gqcpy.CartesianDirection.x]
L_y = left2.center()[gqcpy.CartesianDirection.y]
L_z = left2.center()[gqcpy.CartesianDirection.z]
M_x = right1.center()[gqcpy.CartesianDirection.x]
M_y = right1.center()[gqcpy.CartesianDirection.y]
M_z = right1.center()[gqcpy.CartesianDirection.z]
N_x = right2.center()[gqcpy.CartesianDirection.x]
N_y = right2.center()[gqcpy.CartesianDirection.y]
N_z = right2.center()[gqcpy.CartesianDirection.z]
# Prepare the McMurchie-Davidson coefficients.
E_x = gqcpy.McMurchieDavidsonCoefficient(K_x, a, L_x, b)
E_y = gqcpy.McMurchieDavidsonCoefficient(K_y, a, L_y, b)
E_z = gqcpy.McMurchieDavidsonCoefficient(K_z, a, L_z, b)
E_x_ = gqcpy.McMurchieDavidsonCoefficient(M_x, c, N_x, d)
E_y_ = gqcpy.McMurchieDavidsonCoefficient(M_y, c, N_y, d)
E_z_ = gqcpy.McMurchieDavidsonCoefficient(M_z, c, N_z, d)
# Prepare the Hermite Coulomb integral.
p = a + b;
q = c + d;
alpha = p * q / (p + q);
P = np.array([E_x.centerOfMass(), E_y.centerOfMass(), E_z.centerOfMass()])
Q = np.array([E_x_.centerOfMass(), E_y_.centerOfMass(), E_z_.centerOfMass()])
R = gqcpy.HermiteCoulombIntegral(alpha, P, Q);
# Calculate the Coulomb repulsion integrals over the primitives.
integral = 0.0
for t in range(0, i + j + 1):
for u in range(0, k + l + 1):
for v in range(0, m + n + 1):
for tau in range(0, i_ + j_ + 1):
for mu in range(0, k_ + l_ + 1):
for nu in range(0, m_ + n_ + 1):
# Add the contribution to the integral. The prefactor will be applied at the end.
integral += E_x(i, j, t) * E_y(k, l, u) * E_z(m, n, v) * \
E_x_(i_, j_, tau) * E_y_(k_, l_, mu) * E_z_(m_, n_, nu) * \
np.power(-1, tau + mu + nu) * R(0, t + tau, u + mu, v + nu)
return 2 * np.power(np.pi, 2.5) / (p * q * np.sqrt(p + q)) * integral;
Instead of using a one-electron engine, we’ll have to use a two-electron engine for the Coulomb repulsion integrals.
primitive_coulomb_engine = gqcpy.FunctionalTwoElectronPrimitiveIntegralEngine_d(coulomb_repulsion_function)
coulomb_engine = gqcpy.FunctionalTwoElectronIntegralEngine_d(primitive_coulomb_engine)
g = gqcpy.IntegralCalculator.calculate(coulomb_engine, shell_set, shell_set)
print(g)
[[[[0.77460593 0.56886143]
[0.56886143 0.65017746]]
[[0.56886143 0.45590151]
[0.45590151 0.56886143]]]
[[[0.56886143 0.45590151]
[0.45590151 0.56886143]]
[[0.65017746 0.56886143]
[0.56886143 0.77460593]]]]
g_ref = spin_orbital_basis.quantize(gqcpy.CoulombRepulsionOperator()).parameters()
print(g_ref)
[[[[0.77460593 0.56886143]
[0.56886143 0.65017746]]
[[0.56886143 0.45590151]
[0.45590151 0.56886143]]]
[[[0.56886143 0.45590151]
[0.45590151 0.56886143]]
[[0.65017746 0.56886143]
[0.56886143 0.77460593]]]]