becke-multicenter-integration 0.0.2

Molecular integrals on multicenter grids (Becke's integration scheme)

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)