becke-multicenter-integration 0.0.2

Description:

beckemulticenterintegration 0.0.2

Numerical Molecular Integrals on Multicenter Grids (Becke’s Integration Scheme)
This python package computes molecular integrals numerically for arbitrary
basis functions using Becke’s multicenter grids [1].

The following matrix elements are available:


overlap (a|b)
kinetic energy (a|T|b)
nuclear attraction energy: sum_I -Z(I) * (a|1/rI|b)
dipole operator (a|e*r|b)
electron repulsive integrals (ab|1/r12|cd)



Arbitrary functions can be integrated numerically over space [1].
Apart from this, Poisson’s equation and Laplace’s equation can be solved numerically
for arbitrary charge distributions or wavefunctions [2].

The code is rather slow and only intended for debugging electron integral routines.


Requirements
Required python packages:


numpy, scipy, matplotlib
mpmath




Installation
The package is installed with
$ pip install -e .
in the top directory. To verify the proper functioning of the code
a set of tests should be run with
$ cd tests
$ python -m unittest


Getting Started
First we need to import numpy and the becke module:
import numpy as np
import becke
The multicenter grid is defined by the molecular geometry. Space is partitioned into
fuzzy Voronoi polyhedra. Each atom is the center of a spherical grid and the grids of
all atoms are superimposed. The geometry is defined as a list of tuples (Zat, (X,Y,Z))
where Zat is the atomic number, and X,Y,Z are the cartesian coordinates of the atom
in bohr:
# H2 geometry
atoms = [(1, (0.0, 0.0,-0.5)),
(1, (0.0, 0.0,+0.5))]
The resolution of the multicenter grids is controlled via:
from becke import settings
settings.radial_grid_factor = 3 # increase number of radial points by factor 3
settings.lebedev_order = 23 # angular Lebedev grid of order 23
Wavefunctions are defined as python functions, which take three numpy arrays with the
x-, y- and z-coordinates as input.
# 1s orbital on first hydrogen sA
def aoA(x,y,z):
r = np.sqrt(x**2+y**2+(z+0.5)**2)
return 1.0/np.sqrt(np.pi) * np.exp(-r)

# 1s orbital on second hydrogen sB
def aoB(x,y,z):
r = np.sqrt(x**2+y**2+(z-0.5)**2)
return 1.0/np.sqrt(np.pi) * np.exp(-r)
Now typical one- and two-electron integrals can be calculated for the atomic orbitals
by numerical integration:
# one-electron integrals
print("(a|b)= ", becke.overlap(atoms, aoA, aoB) )
print("(a|b)= ", becke.integral(atoms, lambda x,y,z: aoA(x,y,z)*aoB(x,y,z)) )
print("(a|T|b)= ", becke.kinetic(atoms, aoA, aoB) )
print("(a|V|b)= ", becke.nuclear(atoms, aoA, aoB) )
print("(a|e*r|b)= ", becke.electronic_dipole(atoms, aoA, aoB) )

# two-electron repulsion integrals
print("(aa|bb)= ", becke.electron_repulsion(atoms, aoA, aoA, aoB, aoB) )
print("(ab|ab)= ", becke.electron_repulsion(atoms, aoA, aoB, aoA, aoB) )
When computing the Laplacian or solving the Poisson equation, the return values
are functions themselves that allow to evaluate the Laplacian or electrostatic
potential on a grid (x,y,z):
# __2
# Laplacian lap(x,y,z) = \/ wfn
lap = becke.laplacian(atoms, aoA)
The Laplacian can be used to compute the kinetic energy:
print("(a|T|a)= ", -0.5 * becke.integral(atoms, lambda x,y,z: aoA(x,y,z) * lap(x,y,z) ) )
The following code solves the Poisson equation for the electron density of the
hydrogen atom and plots the electrostatic potential along the z-axis:
s = becke.overlap(atoms, aoA, aoB)
# lowest molecular orbital of hydrogen molecule
def mo(x,y,z):
return (aoA(x,y,z) + aoB(x,y,z))/np.sqrt(2*(1+s))

print("(mo|mo)= ", becke.overlap(atoms, mo, mo) )

# electrostatic potential due to electronic density
v = becke.poisson(atoms, lambda x,y,z: mo(x,y,z)**2)

import matplotlib.pyplot as plt
r = np.linspace(-2.0, 2.0, 100)
plt.plot(r, v(0*r,0*r,r), label=r"$V_{elec}$")

plt.xlabel(r"z / $a_0$")
plt.ylabel(r"electrostatic potential")
plt.show()

References


[1]
(1,2)
A.Becke, “A multicenter numerical integration scheme for polyatomic molecules”,
J.Chem.Phys. 88, 2547 (1988)


[2]
A.Becke, R.Dickson, “Numerical solution of Poisson’s equation in polyatomic molecules”,
J.Chem.Phys. 89, 2993 (1988)


Some useful information is also contained in


[3]
T.Shiozaki, S.Hirata, “Grid-based numerical Hartree-Fock solutions of polyatomic molecules”,
Phys.Rev. A 76, 040503(R) (2007)