Python package for evaluating integrals of Gaussian type orbitals in electronic structure calculations
Project description
PyQInt
Table of Contents
Purpose
PyQInt is a Python package for calculating one- and two-electron integrals as encountered in electronic structure calculations. Since integral evaluation can be quite computationally intensive, they are programmed in C++ and connected to Python using Cython.
Installation
Anaconda
Open Anaconda prompt and type
conda install -c ifilot pyqint
PyPi
Open a terminal and type
pip install pyqint
Usage
Overlap integrals
from pyqint import PyQInt, cgf, gto
import numpy as np
from copy import deepcopy
# construct integrator object
integrator = PyQInt()
# build cgf for hydrogen separated by 1.4 a.u.
cgf1 = cgf([0.0, 0.0, 0.0])
cgf1.add_gto(0.154329, 3.425251, 0, 0, 0)
cgf1.add_gto(0.535328, 0.623914, 0, 0, 0)
cgf1.add_gto(0.444635, 0.168855, 0, 0, 0)
# create a copy of the CGF
cgf2 = deepcopy(cgf1)
cgf2.p[2] = 1.4
# construct empty matrix
S = np.zeros((2,2))
S[0,0] = integrator.overlap(cgf1, cgf1)
S[0,1] = S[1,0] = integrator.overlap(cgf1, cgf2)
S[1,1] = integrator.overlap(cgf2, cgf2)
# output result
print(S)
Dipole integrals
from pyqint import PyQInt, cgf, gto, Molecule
import numpy as np
# construct integrator object
integrator = PyQInt()
# build a water molecule
mol = Molecule("H2O")
mol.add_atom('O', 0.00000, -0.07579, 0.0000, unit='angstrom')
mol.add_atom('H', 0.86681, 0.60144, 0.0000, unit='angstrom')
mol.add_atom('H', -0.86681, 0.60144, 0.0000, unit='angstrom')
cgfs, nuclei = mol.build_basis('sto3g')
# construct empty tensor
N = len(cgfs)
# loop over CGFS and cartesian directions
D = np.zeros((3,N,N))
for k in range(0,3):
for j in range(0,N):
for i in range(0,N):
D[k,j,i] = integrator.dipole(cgfs[i], cgfs[j], k)
# output result
print(D)
Kinetic integrals
from pyqint import PyQInt, cgf, gto
import numpy as np
from copy import deepcopy
# construct integrator object
integrator = PyQInt()
# build cgf for hydrogen separated by 1.4 a.u.
cgf1 = cgf([0.0, 0.0, 0.0])
cgf1.add_gto(0.154329, 3.425251, 0, 0, 0)
cgf1.add_gto(0.535328, 0.623914, 0, 0, 0)
cgf1.add_gto(0.444635, 0.168855, 0, 0, 0)
# create a copy of the CGF
cgf2 = deepcopy(cgf1)
cgf2.p[2] = 1.4
# construct empty matrix
T = np.zeros((2,2))
T[0,0] = integrator.kinetic(cgf1, cgf1)
T[0,1] = T[1,0] = integrator.kinetic(cgf1, cgf2)
T[1,1] = integrator.kinetic(cgf2, cgf2)
# output result
print(T)
Nuclear attraction integrals
from pyqint import PyQInt, cgf, gto
import numpy as np
from copy import deepcopy
# construct integrator object
integrator = PyQInt()
# build cgf for hydrogen separated by 1.4 a.u.
cgf1 = cgf([0.0, 0.0, 0.0])
cgf1.add_gto(0.154329, 3.425251, 0, 0, 0)
cgf1.add_gto(0.535328, 0.623914, 0, 0, 0)
cgf1.add_gto(0.444635, 0.168855, 0, 0, 0)
# create a copy of the CGF
cgf2 = deepcopy(cgf1)
cgf2.p[2] = 1.4
# Build nuclear attraction integrals
V1 = np.zeros((2,2))
V1[0,0] = integrator.nuclear(cgf1, cgf1, cgf1.p, 1)
V1[0,1] = V1[1,0] = integrator.nuclear(cgf1, cgf2, cgf1.p, 1)
V1[1,1] = integrator.nuclear(cgf2, cgf2, cgf1.p, 1)
V2 = np.zeros((2,2))
V2[0,0] = integrator.nuclear(cgf1, cgf1, cgf2.p, 1)
V2[0,1] = V2[1,0] = integrator.nuclear(cgf1, cgf2, cgf2.p, 1)
V2[1,1] = integrator.nuclear(cgf2, cgf2, cgf2.p, 1)
# print result
print(V1,V2)
Two-electron integrals
from pyqint import PyQInt, cgf, gto
import numpy as np
from copy import deepcopy
# construct integrator object
integrator = PyQInt()
# build cgf for hydrogen separated by 1.4 a.u.
cgf1 = cgf([0.0, 0.0, 0.0])
cgf1.add_gto(0.154329, 3.425251, 0, 0, 0)
cgf1.add_gto(0.535328, 0.623914, 0, 0, 0)
cgf1.add_gto(0.444635, 0.168855, 0, 0, 0)
# create a copy of the CGF
cgf2 = deepcopy(cgf1)
cgf2.p[2] = 1.4
T1111 = integrator.repulsion(cgf1, cgf1, cgf1, cgf1)
T1122 = integrator.repulsion(cgf1, cgf1, cgf2, cgf2)
T1112 = integrator.repulsion(cgf1, cgf1, cgf1, cgf2)
T2121 = integrator.repulsion(cgf2, cgf1, cgf2, cgf1)
T1222 = integrator.repulsion(cgf1, cgf2, cgf2, cgf2)
T2211 = integrator.repulsion(cgf2, cgf2, cgf1, cgf1)
print(T1111)
print(T1122)
print(T1112)
print(T2121)
print(T1222)
print(T2211)
Construction of Contracted Gaussian Functions
from pyqint import PyQInt, Molecule
import numpy as np
# construct integrator object
integrator = PyQInt()
# build hydrogen molecule
mol = Molecule()
mol.add_atom('H', 0.0, 0.0, 0.0)
mol.add_atom('H', 0.0, 0.0, 1.4)
cgfs, nuclei = mol.build_basis('sto3g')
print(cgfs, nuclei)
Parallel evaluation of integrals
From a collection of Contracted Gaussian Functions, the complete set of overlap, kinetic, nuclear attraction and two-electron integrals
can be quickly evaluated using the build_integrals
function. Using the npar
argument, the number of threads to be spawned can be set.
from pyqint import PyQInt, Molecule
import numpy as np
import multiprocessing
# construct integrator object
integrator = PyQInt()
# build hydrogen molecule
mol = Molecule()
mol.add_atom('H', 0.0, 0.0, 0.0)
mol.add_atom('H', 0.0, 0.0, 1.4)
cgfs, nuclei = mol.build_basis('sto3g')
# evaluate all integrals
ncpu = multiprocessing.cpu_count()
S, T, V, teint = integrator.build_integrals(cgfs, nuclei, npar=ncpu, verbose=False)
print(S, T, V, teint)
Plotting molecular orbitals
from pyqint import PyQInt, Molecule
import matplotlib.pyplot as plt
import numpy as np
# coefficients (calculated by Hartree-Fock using a sto-3g basis set)
coeff = [8.37612e-17, -2.73592e-16, -0.713011, -1.8627e-17, 9.53496e-17, -0.379323, 0.379323]
# construct integrator object
integrator = PyQInt()
# build water molecule
mol = Molecule('H2O')
mol.add_atom('O', 0.0, 0.0, 0.0)
mol.add_atom('H', 0.7570, 0.5860, 0.0)
mol.add_atom('H', -0.7570, 0.5860, 0.0)
cgfs, nuclei = mol.build_basis('sto3g')
# build grid
x = np.linspace(-2, 2, 50)
y = np.linspace(-2, 2, 50)
xx, yy = np.meshgrid(x,y)
zz = np.zeros(len(x) * len(y))
grid = np.vstack([xx.flatten(), yy.flatten(), zz]).reshape(3,-1).T
res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x)))
# plot wave function
plt.imshow(res, origin='lower', extent=[-2,2,-2,2], cmap='PiYG')
plt.colorbar()
plt.title('1b$_{2}$ Molecular orbital of H$_{2}$O')
Performing Hartree-Fock calculations
from pyqint import PyQInt, Molecule, HF
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
def main():
# calculate sto-3g coefficients for h2o
cgfs, coeff = calculate_co()
# visualize orbitals
fig, ax = plt.subplots(2,3, figsize=(18,10))
for i in range(0,2):
for j in range(0,3):
dens = plot_wavefunction(cgfs, coeff[:,i*3+j])
limit = max(abs(np.min(dens)), abs(np.max(dens)) )
im = ax[i,j].imshow(dens, origin='lower', interpolation='bilinear',
extent=[-2,2,-2,2], cmap='PiYG', vmin=-limit, vmax=limit)
ax[i,j].set_xlabel('Distance a.u.')
ax[i,j].set_ylabel('Distance a.u.')
divider = make_axes_locatable(ax[i,j])
cax = divider.append_axes('right', size='5%', pad=0.05)
fig.colorbar(im, cax=cax, orientation='vertical')
def calculate_co():
mol = Molecule()
mol.add_atom('C', 0.0, -0.5, 0.0)
mol.add_atom('O', 0.0, 0.5, 0.0)
result = HF().rhf(mol, 'sto3g')
return result['cgfs'], result['orbc']
def plot_wavefunction(cgfs, coeff):
# build integrator
integrator = PyQInt()
# build grid
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
xx, yy = np.meshgrid(x,y)
zz = np.zeros(len(x) * len(y))
grid = np.vstack([xx.flatten(), yy.flatten(), zz]).reshape(3,-1).T
res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(y), len(x)))
return res
if __name__ == '__main__':
main()
Constructing isosurfaces
By installing the PyTessel package, you can construct an isosurface from a wave function scalar field.
from pyqint import PyQInt, Molecule, HF
import numpy as np
from pytessel import PyTessel
def main():
# calculate sto-3g coefficients for h2o
cgfs, coeff = calculate_co()
# build isosurface of the fifth MO
# isovalue = 0.1
# store result as .ply file
build_isosurface('co_04.ply', cgfs, coeff[:,4], 0.1)
def build_isosurface(filename, cgfs, coeff, isovalue):
# generate some data
sz = 100
integrator = PyQInt()
grid = integrator.build_rectgrid3d(-5, 5, sz)
scalarfield = np.reshape(integrator.plot_wavefunction(grid, coeff, cgfs), (sz, sz, sz))
unitcell = np.diag(np.ones(3) * 10.0)
pytessel = PyTessel()
vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), isovalue)
pytessel.write_ply(filename, vertices, normals, indices)
def calculate_co():
mol = Molecule()
mol.add_atom('C', 0.0, -0.5, 0.0)
mol.add_atom('O', 0.0, 0.5, 0.0)
result = HF().rhf(mol, 'sto3g')
return result['cgfs'], result['orbc']
if __name__ == '__main__':
main()
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distributions
Built Distributions
Hashes for pyqint-0.10.0.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | d97b72913c98ee0ef3aa0096107784ff6dcc6db0b44a9d2c1c1e162fc7d4fefd |
|
MD5 | 607f3a1058f0f875573feadeb5180fb0 |
|
BLAKE2b-256 | f6e539d31bf0c53eba48b27e3b99cdfbcf355d62caf3b893e81f4aa958fd8060 |
Hashes for pyqint-0.10.0.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 0eb56976a70350c228e9cc3242900b8454f78136e7e4e86d14449b13f7e59102 |
|
MD5 | 847a878b51ffb80c58c2978392e05d37 |
|
BLAKE2b-256 | 2dfb794f69172c911f6c0b2662dcf5df672239c1db8511c9bb613f945474f212 |
Hashes for pyqint-0.10.0.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 3e1fcc682f200f7feb7a037ef55b3c316f265d4ecfb85577a76b11c2d62b2936 |
|
MD5 | ce5143a6f919d032a32a42f55bd9af2e |
|
BLAKE2b-256 | cd50741b6a3cc245da29993d7d927b64e140de394c0e6b483dfd50762d5d7fa1 |
Hashes for pyqint-0.10.0.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 928a87581b3ab0c01181108c08c72de9a2ba7fdcfb1dcd330f57d88b0a18017b |
|
MD5 | a9ac82f6d7ab9a8a69f39305735c026b |
|
BLAKE2b-256 | 98b7ba14504fb8a9e1706c5b514680504f3c7cc3a0c34c740cd1dcab224b95c3 |
Hashes for pyqint-0.10.0.0-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 3f101b246774657b2d84bc7d3965128a454f0cf9591f9d176c653666030a6698 |
|
MD5 | 3221b08b6d39f57e16698173a2c7cd4c |
|
BLAKE2b-256 | 65d50bb0fc3536e150db5ab51ccfa6c1da601f4d31685b8b00a57a4cc10c5ee4 |