# basic Hartee-Fock program to compute the total energies and some properties of a molecule

I am trying to write a basic Hartee-Fock program that computes the total energies and some properties of a molecule given as an input. The code uses Cython to statically compile the computationally expensive parts, and there is a basic parallelization implemented using mpi4py. However, the code is still very slow for all practical applications.
The code is currently hosted at GitHub, including a sample input for a water molecule.
Any suggestions on how this can be improved, especially execution time, are much appreciated.
Note: I am largely self-taught in Python. So any pointers to where I can learn better coding practices are also appreciated.

driver code

from numpy                      import array, float32, int32, zeros, pi, exp, dot, savetxt, diag, trace
from numpy.linalg               import norm, eigh, solve
from math                       import ceil, sqrt, log
from scipy.special              import factorial2, comb
from sys                        import argv
from itertools                  import combinations_with_replacement, combinations
from mpi4py                     import MPI
from plank.routines.oneelectron import overlapcgtos, kineticcgtos
from plank.routines.twoelectron import nuclearcgtos, ericgtos
from plank.routines.hatreefock  import computecorehamiltonian, canonicalorthogonalization, computehamiltonian
from plank.routines.hatreefock  import restrictedhatreefock, restrictedhatreefockDIIS, totalenergy
from plank.routines.hatreefock  import diisorthogonalization
from yaml                       import safe_load

with open("periodicTable.yaml") as periodic:
periodicTable   =   {y: x[y] for x in periodictable for y in x}

def getatom(atomname):
return periodicTable[atomname]

def nuclearenergy(atomobjects):
atompairs       =   list(combinations(atomobjects,2))
nuclearenergy   =   0.0
for atompair in atompairs:
center1         =   atompair[0].center
charge1         =   atompair[0].charge
center2         =   atompair[1].center
charge2         =   atompair[1].charge
distance        =   norm(center1-center2)
energy          =   charge1*charge2/distance
nuclearenergy   =   nuclearenergy+energy
return nuclearenergy

def overlap(basisobjects):
overlaps    =   []
for basisobjecttuple in basisobjects:
result  =   overlapcgtos(basisobjecttuple[0], basisobjecttuple[1])
overlaps.append((result, basisobjecttuple[0].basisindex, basisobjecttuple[1].basisindex))
return overlaps

def kinetic(basisobjects):
kineticintegrals =   []
for basisobjecttuple in basisobjects:
result  =   kineticcgtos(basisobjecttuple[0], basisobjecttuple[1])
kineticintegrals.append((result, basisobjecttuple[0].basisindex, basisobjecttuple[1].basisindex))
return kineticintegrals

def nuclear(basisobjects, atomobjects):
eniintegrals    =   []
for basisobjecttuple in basisobjects:
result      =   nuclearcgtos(basisobjecttuple[0], basisobjecttuple[1], atomobjects)
eniintegrals.append((result, basisobjecttuple[0].basisindex, basisobjecttuple[1].basisindex))
return eniintegrals

def electronic(basisobjects):
eriintegrals    =   []
for basisobjectquartet in basisobjects:
result  =   ericgtos(basisobjectquartet[0], basisobjectquartet[1], basisobjectquartet[2], basisobjectquartet[3])
eriintegrals.append((result, basisobjectquartet[0].basisindex, basisobjectquartet[1].basisindex, basisobjectquartet[2].basisindex, basisobjectquartet[3].basisindex))
return eriintegrals

def integralEngine(molecule, basisobjects, dimensions, calctype="overlap"):
results     =   []
if calctype == "overlap":
rank            =   comm.Get_rank()
nprocs          =   comm.Get_size()
shellpairs      =   list(combinations_with_replacement(basisobjects, 2))
nshellpairs     =   len(shellpairs)
nshellpair      =   ceil(nshellpairs/nprocs)
shellpairi      =   rank*nshellpair
shellpairf      =   (rank+1)*nshellpair
overlapresults  =   overlap(shellpairs[shellpairi:shellpairf])

if rank != 0:
comm.send(overlapresults, dest=0)
else:
for result in overlapresults:
results.append(result)
for proc in range(1, nprocs):
overlapresults  =   comm.recv(source=proc)
for result in overlapresults:
results.append(result)

if rank == 0:
molecule.overlapmat =   zeros((dimensions, dimensions))
for data in results:
molecule.overlapmat[data[1], data[2]]   =   data[0]
molecule.overlapmat[data[2], data[1]]   =   data[0]

molecule.overlapmat =   comm.bcast(molecule.overlapmat, root=0)

if calctype == "kinetic":
rank            =   comm.Get_rank()
nprocs          =   comm.Get_size()
shellpairs      =   list(combinations_with_replacement(basisobjects, 2))
nshellpairs     =   len(shellpairs)
nshellpair      =   ceil(nshellpairs/nprocs)
shellpairi      =   rank*nshellpair
shellpairf      =   (rank+1)*nshellpair
kineticresults  =   kinetic(shellpairs[shellpairi:shellpairf])

if rank != 0:
comm.send(kineticresults, dest=0)
else:
for result in kineticresults:
results.append(result)
for proc in range(1, nprocs):
kineticresults  =   comm.recv(source=proc)
for result in kineticresults:
results.append(result)

if rank == 0:
molecule.kineticmat =   zeros((dimensions, dimensions))
for data in results:
molecule.kineticmat[data[1], data[2]]   =   data[0]
molecule.kineticmat[data[2], data[1]]   =   data[0]

molecule.kineticmat =   comm.bcast(molecule.kineticmat, root=0)

if calctype == "eni":
rank            =   comm.Get_rank()
nprocs          =   comm.Get_size()
shellpairs      =   list(combinations_with_replacement(basisobjects, 2))
atomobjects     =   [x for x in molecule.geometry]
nshellpairs     =   len(shellpairs)
nshellpair      =   ceil(nshellpairs/nprocs)
shellpairi      =   rank*nshellpair
shellpairf      =   (rank+1)*nshellpair
eniresults      =   nuclear(shellpairs[shellpairi:shellpairf], atomobjects)

if rank != 0:
comm.send(eniresults, dest=0)
else:
for result in eniresults:
results.append(result)
for proc in range(1, nprocs):
eniresults  =   comm.recv(source=proc)
for result in eniresults:
results.append(result)

if rank == 0:
molecule.enimat =   zeros((dimensions, dimensions))
for data in results:
molecule.enimat[data[1], data[2]]   =   data[0]
molecule.enimat[data[2], data[1]]   =   data[0]

molecule.enimat =   comm.bcast(molecule.enimat, root=0)

if calctype == "eri":
rank            =   comm.Get_rank()
nprocs          =   comm.Get_size()
shellquartets  =   []
for basisobjectA in basisObjects:
for basisobjectB in basisObjects:
for basisobjectC in basisObjects:
for basisobjectD in basisObjects:
index1  =   basisobjectA.basisindex
index2  =   basisobjectB.basisindex
index3  =   basisobjectC.basisindex
index4  =   basisobjectD.basisindex
index12 =   (index1*(index1+1)//2)+index2
index34 =   (index3*(index3+1)//2)+index4
if index12 >= index34:
shellquartets.append((basisobjectA, basisobjectB, basisobjectC, basisobjectD))
nshellquartets  =   len(shellquartets)
nshellquartet   =   ceil(nshellquartets/nprocs)
shellpairi      =   rank*nshellquartet
shellpairf      =   (rank+1)*nshellquartet
eriresults      =   electronic(shellquartets[shellpairi:shellpairf])

if rank != 0:
comm.send(eriresults, dest=0)
else:
for result in eriresults:
results.append(result)
for proc in range(1, nprocs):
eriresults  =   comm.recv(source=proc)
for result in eriresults:
results.append(result)

if rank == 0:
molecule.erimat =   zeros((dimensions, dimensions, dimensions, dimensions))
for result in results:
molecule.erimat[result[1], result[2], result[3], result[4]]  =   result[0]
molecule.erimat[result[3], result[4], result[1], result[2]]  =   result[0]
molecule.erimat[result[2], result[1], result[4], result[3]]  =   result[0]
molecule.erimat[result[4], result[3], result[2], result[1]]  =   result[0]
molecule.erimat[result[2], result[1], result[3], result[4]]  =   result[0]
molecule.erimat[result[4], result[3], result[1], result[2]]  =   result[0]
molecule.erimat[result[1], result[2], result[4], result[3]]  =   result[0]
molecule.erimat[result[3], result[4], result[2], result[1]]  =   result[0]

molecule.erimat =   comm.bcast(molecule.erimat, root=0)

def scfEngine(molecule, dimensions):
density     =   zeros((dimensions, dimensions))
orthomat    =   canonicalorthogonalization(molecule.overlapmat)
diisortho   =   diisorthogonalization(molecule.overlapmat)
core        =   computecorehamiltonian(molecule.kineticmat, molecule.enimat)
for iteration in range(1, molecule.maxiterations):
if iteration <= 2 or molecule.DIIS == False:
orthofock, fock, hamiltonian, orbitals, coeffs, density, rmsdensity, maxdensity =   restrictedhatreefock(molecule.nelectrons, core, diisortho, density, molecule.erimat)
molenergy   =   totalenergy(molecule.nuclearenergy, fock, density, core)
print("{:10.0f}{:20.5f}{:20.5f}{:20.5f}".format(iteration, molenergy, log(rmsdensity), log(maxdensity)))
if iteration > 2 and molecule.DIIS == True:
orthofock, fock, hamiltonian, orbitals, coeffs, density, rmsdensity, maxdensity =   restrictedhatreefockDIIS(molecule, molecule.nelectrons, core, diisortho, density, molecule.erimat)
molenergy   =   totalenergy(molecule.nuclearenergy, fock, density, core)
print("{:10.0f}{:20.5f}{:20.5f}{:20.5f}".format(iteration, molenergy, log(rmsdensity), log(maxdensity)))
if log(rmsdensity) <= -20.0 and log(maxdensity) <= -20.0:
molecule.scfstatus  =   True
molecule.density    =   density
molecule.fock       =   fock
molecule.orbitals   =   orbitals
molecule.coeffs     =   coeffs
print("SCF CONVERGENCE ACHIEVED IN {:10.0f} ITERATIONS".format(iteration))
break

class Basis(object):
""" docstring for Basis

Basis class constructs one basis object of type (s,p,d or f) depending on the basis set
passed to the code. The generation of basis object happens during the construction of the
atom object in the class Atom, with the list of basis objects tied to the atom object.
"""
def __init__(self, shell, exponents, coefficients, center):
super(Basis, self).__init__()
self.shell          =   array(shell)
self.exponents      =   array(exponents)
self.coefficients   =   array(coefficients)
self.normcoeffs     =   zeros(self.coefficients.size)
self.center         =   array(center)
self.basisindex     =   0

self.normalizepGTO()

def normalizepGTO(self):
ll, mm, nn          =   self.shell
totalangmomentum    =   sum(self.shell)
prefactorpGTO       =   pow(2, 2*totalangmomentum)*pow(2, 1.5)/factorial2(2*ll-1)/factorial2(2*mm-1)/factorial2(2*nn-1)/pow(pi, 1.5)

for index, exponent in enumerate(self.exponents):
self.normcoeffs[index]  =   sqrt(pow(exponent, totalangmomentum)*pow(exponent, 1.5)*prefactorpGTO)

prefactorcGTO   =   pow(pi, 1.5)*factorial2(2*ll-1)*factorial2(2*mm-1)*factorial2(2*nn-1)/pow(2.0, totalangmomentum)
normalfactor    =   0.0

for index1, coefficient1 in enumerate(self.coefficients):
for index2, coefficient2 in enumerate(self.coefficients):
t_normalfactor  =   (self.normcoeffs[index1]*self.normcoeffs[index2]*self.coefficients[index1]*self.coefficients[index2])/(pow(self.exponents[index1]+self.exponents[index2], totalangmomentum+1.5))
normalfactor    =   normalfactor + t_normalfactor

normalfactor    =   prefactorcGTO*normalfactor
normalfactor    =   pow(normalfactor, -0.5)

for index, coefficient in enumerate(self.coefficients):
self.coefficients[index]  =   self.coefficients[index]*normalfactor

class Atom(object):
"""docstring for Atom."""

def __init__(self, atomname, center, charge, basisset='sto-3g'):
super(Atom, self).__init__()
self.atomname   =   atomname
self.center     =   array(center)*1.8897259886
self.basisset   =   basisset
self.shells     =   []
self.charge     =   charge

self.basisfile  =   "./basis/"+self.atomname+"."+self.basisset+".gbs"
with open(self.basisfile) as target:
assert (_atomname == self.atomname), "Wrong basis set"
for lnumber, line in enumerate(basisdata[1:]):
if "S" in line and "P" not in line:
nprims          =   int(line.split()[1])
pgtodata        =   [x.replace('D', 'E').split() for x in basisdata[lnumber+2:lnumber+2+nprims]]
exponents       =   [float(x[0]) for x in pgtodata]
coefficients    =   [float(x[1]) for x in pgtodata]
shell00         =   [0, 0, 0]
self.shells.append(Basis(shell00, exponents, coefficients, self.center))

if "P" in line:
nprims          =   int(line.split()[1])
pgtodata        =   [x.replace('D', 'E').split() for x in basisdata[lnumber+2:lnumber+2+nprims]]
coefficient1    =   [float(x[1]) for x in pgtodata]
coefficient2    =   [float(x[2]) for x in pgtodata]
exponents       =   [float(x[0]) for x in pgtodata]
shell00         =   [0, 0, 0]
self.shells.append(Basis(shell00, exponents, coefficient1, self.center))
shell11         =   [1, 0, 0]
self.shells.append(Basis(shell11, exponents, coefficient2, self.center))
shell12         =   [0, 1, 0]
self.shells.append(Basis(shell12, exponents, coefficient2, self.center))
shell13         =   [0, 0, 1]
self.shells.append(Basis(shell13, exponents, coefficient2, self.center))

if "D" in line and "+" not in line:
nprims          =   int(line.split()[1])
pgtodata        =   [x.replace('D', 'E').split() for x in basisdata[lnumber+2:lnumber+2+nprims]]
exponents       =   [float(x[0]) for x in pgtodata]
coefficients    =   [float(x[1]) for x in pgtodata]
shell20         =   [2, 0, 0]
self.shells.append(Basis(shell20, exponents, coefficients, self.center))
shell21         =   [1, 1, 0]
self.shells.append(Basis(shell21, exponents, coefficients, self.center))
shell22         =   [1, 0, 1]
self.shells.append(Basis(shell22, exponents, coefficients, self.center))
shell23         =   [0, 2, 0]
self.shells.append(Basis(shell23, exponents, coefficients, self.center))
shell24         =   [0, 1, 1]
self.shells.append(Basis(shell24, exponents, coefficients, self.center))
shell25         =   [0, 0, 2]
self.shells.append(Basis(shell25, exponents, coefficients, self.center))

class Geometry(object):
"""docstring for Geometry."""

def __init__(self, inputfile):
super(Geometry, self).__init__()
self.inputfile      =   inputfile
self.charge         =   0
self.natoms         =   0
self.nelectrons     =   0
self.calctype       =   'energy'
self.basisset       =   'sto-3g'
self.geometry       =   []
self.atomobjects    =   []
self.density        =   []
self.scfstatus      =   False
self.maxiterations  =   100
self.overlapmat     =   []
self.kineticmat     =   []
self.enimat         =   []
self.erimat         =   []
self.nuclearenergy  =   0.0
self.totalenergy    =   0.0
self.diisfock       =   []
self.diiserror      =   []
self.DIIS

with open(self.inputfile) as file:
calcbasis       =   inputdata[0].split()
self.calctype   =   calcbasis[0]
self.basisset   =   calcbasis[1]
chargemul       =   inputdata[1].split()
self.charge     =   int(chargemul[0])
multiplicity    =   int(chargemul[1])
try:
assert multiplicity == 1, "[ERROR] : Only closed-shell singlets are currently supported."
except AssertionError as Message:
exit(Message)
natoms          =   int(inputdata[2].split()[0])
coorddata       =   [x for x in inputdata[3:3+natoms]]
for atom in coorddata:
atomdata        =   atom.split()
atomname        =   atomdata[0]
charge          =   getatom(atomname)
coords          =   [float(x) for x in atomdata[1:]]
self.nelectrons +=  charge
self.geometry.append(Atom(atomname, coords, charge, self.basisset))
iterdiisline     =   inputdata[-1].split()
if iterdiisline  != []:
self.maxiterations  =   int(iterdiisline[0])
self.DIIS           =   True if (iterdiisline[1] == 'T') else False

inputfilename       =   argv[1]
logfilename         =   inputfilename[:-4]+".log"
molecule            =   Geometry(inputfilename)
comm                =   MPI.COMM_WORLD
atomObjects         =   [x for x in molecule.geometry]
basisObjects        =   [y for x in atomObjects for y in x.shells]
nbasisObjects       =   len(basisObjects)
for counter, x in enumerate(basisObjects):
x.basisindex    =   counter

integralEngine(molecule, basisObjects, nbasisObjects, "overlap")
integralEngine(molecule, basisObjects, nbasisObjects, "kinetic")
integralEngine(molecule, basisObjects, nbasisObjects, "eni")
integralEngine(molecule, basisObjects, nbasisObjects, "eri")

rank    =   comm.Get_rank()
if rank == 0:
with open("overlap."+inputfilename[:-4]+".txt", "w") as t:
for row in molecule.overlapmat:
for column in row:
t.write("{:10.3f}".format(column))
t.write("\n")

with open("kinetic."+inputfilename[:-4]+".txt", "w") as t:
for row in molecule.kineticmat:
for column in row:
t.write("{:10.3f}".format(column))
t.write("\n")

with open("nuclear."+inputfilename[:-4]+".txt", "w") as t:
for row in molecule.enimat:
for column in row:
t.write("{:10.3f}".format(column))
t.write("\n")

with open("electronic."+inputfilename[:-4]+".txt", "w") as t:
for axis1 in molecule.erimat:
for axis2 in axis1:
for axis3 in axis2:
for axis4 in axis3:
t.write("{:10.3f}".format(axis4))
t.write("\n")
t.write("\n")
t.write("\n")

molecule.nuclearenergy  =   nuclearenergy(atomObjects)
scfEngine(molecule, nbasisObjects)

with open("density."+inputfilename[:-4]+".txt", "w") as d:
savetxt(d, molecule.density, "%10.3f")

MPI.Finalize()


cython sections

oneelectron.pyx

import cython   as  cython
import numpy    as  np
cimport numpy   as  np

from scipy.special  import  comb, factorial2
from libc.math      cimport exp,  pow
from numpy          import  dot,  pi

@cython.boundscheck(False)
@cython.wraparound(False)
def gaussiantheorem(np.ndarray[np.float64_t] center1, double exponent1, np.ndarray[np.float64_t] center2, double exponent2):
cdef np.ndarray gaussiancenter  =   np.zeros(3)
cdef double gaussianexponent    =   0.0
cdef double gaussianintegral    =   0.0
gaussiancenter      =   ((exponent1*center1)+(exponent2*center2))/(exponent1+exponent2)
gaussianexponent    =   (exponent1*exponent2)/(exponent1+exponent2)
gaussianintegral    =   exp(-1*gaussianexponent*dot(center1-center2, center1-center2))
return gaussiancenter, gaussianintegral

@cython.boundscheck(False)
@cython.wraparound(False)
cdef double overlappgtos(double center1, double exponent1, int shell1, double center2, double exponent2, int shell2, double gaussiancenter):
cdef double t_overlap   =   0.0
cdef double auxiliary   =   0.0
cdef int counter1, counter2
for counter1 in range(0, shell1+1):
for counter2 in range(0, shell2+1):
if (counter1+counter2)%2 == 0:
auxiliary   =   comb(shell1, counter1)
auxiliary   =   auxiliary*comb(shell2, counter2)
auxiliary   =   auxiliary*factorial2(counter1+counter2-1)
auxiliary   =   auxiliary*pow(gaussiancenter-center1, shell1-counter1)
auxiliary   =   auxiliary*pow(gaussiancenter-center2, shell2-counter2)
auxiliary   =   auxiliary/pow(2*(exponent1+exponent2), 0.5*(counter1+counter2))
t_overlap   =   t_overlap+auxiliary
return t_overlap

@cython.boundscheck(False)
@cython.wraparound(False)
def overlapcgtos(basisobject1, basisobject2):
cdef np.ndarray center1         =   basisobject1.center
cdef np.ndarray exponents1      =   basisobject1.exponents
cdef np.ndarray coefficients1   =   basisobject1.coefficients
cdef np.ndarray shell1          =   basisobject1.shell
cdef np.ndarray normcoeffs1     =   basisobject1.normcoeffs

cdef np.ndarray center2         =   basisobject2.center
cdef np.ndarray exponents2      =   basisobject2.exponents
cdef np.ndarray coefficients2   =   basisobject2.coefficients
cdef np.ndarray shell2          =   basisobject2.shell
cdef np.ndarray normcoeffs2     =   basisobject2.normcoeffs

cdef double overlaptotal        =   0.0
cdef index1                     =   0
cdef index2                     =   0
cdef exponent1                  =   0
cdef exponent2                  =   0
cdef np.ndarray gaussiancenter  =   np.zeros(3)
cdef double gaussianintegral    =   0.0
for index1, exponent1 in enumerate(exponents1):
for index2, exponent2 in enumerate(exponents2):
gaussiancenter, gaussianintegral    =   gaussiantheorem(center1, exponent1, center2, exponent2)
overlapx        =   overlappgtos(center1[0], exponent1, shell1[0], center2[0], exponent2, shell2[0], gaussiancenter[0])
overlapy        =   overlappgtos(center1[1], exponent1, shell1[1], center2[1], exponent2, shell2[1], gaussiancenter[1])
overlapz        =   overlappgtos(center1[2], exponent1, shell1[2], center2[2], exponent2, shell2[2], gaussiancenter[2])
overlap         =   overlapx*overlapy*overlapz*gaussianintegral*pow(pi/(exponent1+exponent2), 1.5)
overlaptotal    =   overlaptotal+(normcoeffs1[index1]*normcoeffs2[index2]*coefficients1[index1]*coefficients2[index2])*overlap
return overlaptotal

@cython.boundscheck(False)
@cython.wraparound(False)
def kineticcgtos(basisobject1, basisobject2):
cdef np.ndarray center1         =   basisobject1.center
cdef np.ndarray exponents1      =   basisobject1.exponents
cdef np.ndarray coefficients1   =   basisobject1.coefficients
cdef np.ndarray shell1          =   basisobject1.shell
cdef np.ndarray normcoeffs1     =   basisobject1.normcoeffs

cdef np.ndarray center2         =   basisobject2.center
cdef np.ndarray exponents2      =   basisobject2.exponents
cdef np.ndarray coefficients2   =   basisobject2.coefficients
cdef np.ndarray shell2          =   basisobject2.shell
cdef np.ndarray normcoeffs2     =   basisobject2.normcoeffs

cdef double Kx                  =   0.0
cdef double Ky                  =   0.0
cdef double Kz                  =   0.0
cdef np.ndarray gaussiancenter  =   np.zeros(3)
cdef double gaussianintegral    =   0.0
for index1, exponent1 in enumerate(exponents1):
for index2, exponent2 in enumerate(exponents2):
gaussiancenter, gaussianintegral    =   gaussiantheorem(center1, exponent1, center2, exponent2)
overlapx    =   overlappgtos(center1[0], exponent1, shell1[0], center2[0], exponent2, shell2[0], gaussiancenter[0])
overlapy    =   overlappgtos(center1[1], exponent1, shell1[1], center2[1], exponent2, shell2[1], gaussiancenter[1])
overlapz    =   overlappgtos(center1[2], exponent1, shell1[2], center2[2], exponent2, shell2[2], gaussiancenter[2])

overlapx11  =   overlappgtos(center1[0], exponent1, shell1[0]-1, center2[0], exponent2, shell2[0]-1, gaussiancenter[0])
overlapx12  =   overlappgtos(center1[0], exponent1, shell1[0]+1, center2[0], exponent2, shell2[0]-1, gaussiancenter[0])
overlapx13  =   overlappgtos(center1[0], exponent1, shell1[0]-1, center2[0], exponent2, shell2[0]+1, gaussiancenter[0])
overlapx14  =   overlappgtos(center1[0], exponent1, shell1[0]+1, center2[0], exponent2, shell2[0]+1, gaussiancenter[0])

overlapy11  =   overlappgtos(center1[1], exponent1, shell1[1]-1, center2[1], exponent2, shell2[1]-1, gaussiancenter[1])
overlapy12  =   overlappgtos(center1[1], exponent1, shell1[1]+1, center2[1], exponent2, shell2[1]-1, gaussiancenter[1])
overlapy13  =   overlappgtos(center1[1], exponent1, shell1[1]-1, center2[1], exponent2, shell2[1]+1, gaussiancenter[1])
overlapy14  =   overlappgtos(center1[1], exponent1, shell1[1]+1, center2[1], exponent2, shell2[1]+1, gaussiancenter[1])

overlapz11  =   overlappgtos(center1[2], exponent1, shell1[2]-1, center2[2], exponent2, shell2[2]-1, gaussiancenter[2])
overlapz12  =   overlappgtos(center1[2], exponent1, shell1[2]+1, center2[2], exponent2, shell2[2]-1, gaussiancenter[2])
overlapz13  =   overlappgtos(center1[2], exponent1, shell1[2]-1, center2[2], exponent2, shell2[2]+1, gaussiancenter[2])
overlapz14  =   overlappgtos(center1[2], exponent1, shell1[2]+1, center2[2], exponent2, shell2[2]+1, gaussiancenter[2])

kx          =   shell1[0]*shell2[0]*overlapx11
kx          +=  -2*exponent1*shell2[0]*overlapx12
kx          +=  -2*exponent2*shell1[0]*overlapx13
kx          +=  4*exponent1*exponent2*overlapx14

ky          =   shell1[1]*shell2[1]*overlapy11
ky          +=  -2*exponent1*shell2[1]*overlapy12
ky          +=  -2*exponent2*shell1[1]*overlapy13
ky          +=  4*exponent1*exponent2*overlapy14

kz          =   shell1[2]*shell2[2]*overlapz11
kz          +=  -2*exponent1*shell2[2]*overlapz12
kz          +=  -2*exponent2*shell1[2]*overlapz13
kz          +=  4*exponent1*exponent2*overlapz14

Kx          +=  0.5*kx*overlapy*overlapz*gaussianintegral*pow(pi/(exponent1+exponent2), 1.5)*normcoeffs1[index1]*coefficients1[index1]*normcoeffs2[index2]*coefficients2[index2]
Ky          +=  0.5*ky*overlapx*overlapz*gaussianintegral*pow(pi/(exponent1+exponent2), 1.5)*normcoeffs1[index1]*coefficients1[index1]*normcoeffs2[index2]*coefficients2[index2]
Kz          +=  0.5*kz*overlapx*overlapy*gaussianintegral*pow(pi/(exponent1+exponent2), 1.5)*normcoeffs1[index1]*coefficients1[index1]*normcoeffs2[index2]*coefficients2[index2]
return Kx+Ky+Kz


twoelectron.pyx

import cython   as  cython
import numpy    as  np
cimport numpy   as  np

from scipy.special  import  comb, factorial2, factorial, hyp1f1
from libc.math      cimport exp,  pow, sqrt
from numpy          import  dot,  pi
from numpy.linalg   import norm

@cython.boundscheck(False)
@cython.wraparound(False)
def gaussiantheorem(np.ndarray[np.float64_t] center1, double exponent1, np.ndarray[np.float64_t] center2, double exponent2):
cdef np.ndarray gaussiancenter  =   np.zeros(3)
cdef double gaussianexponent    =   0.0
cdef double gaussianintegral    =   0.0
gaussiancenter      =   ((exponent1*center1)+(exponent2*center2))/(exponent1+exponent2)
gaussianexponent    =   (exponent1*exponent2)/(exponent1+exponent2)
gaussianintegral    =   exp(-1*gaussianexponent*dot(center1-center2, center1-center2))
return gaussiancenter, gaussianintegral

@cython.boundscheck(False)
@cython.wraparound(False)
cdef double expansioncoeff1(int expansionindex, int shell1, float center1, int shell2, float center2, float gaussiancenter):
cdef double t_expansioncoeff    =   0.0
cdef int counter                =   0
cdef double auxiliary           =   0.0
for counter in range(max(0, expansionindex-shell2), min(expansionindex, shell1)+1):
auxiliary           =   comb(shell1, counter)
auxiliary           =   auxiliary*comb(shell2, expansionindex-counter)
auxiliary           =   auxiliary*pow(gaussiancenter-center1, shell1-counter)
auxiliary           =   auxiliary*pow(gaussiancenter-center2, shell2+counter-expansionindex)
t_expansioncoeff    =   t_expansioncoeff+auxiliary
return t_expansioncoeff

@cython.boundscheck(False)
@cython.wraparound(False)
cdef double expansioncoeff2(int expansionindex1, int expansionindex2, int expansionindex3, int shell1, float center1, int shell2, float center2, float atomcenter, float gaussiancenter, float gamma):
cdef double t_expansioncoeff    =   0.0
cdef double epsilon             =   1.0/(4*gamma)
t_expansioncoeff                =   expansioncoeff1(expansionindex1, shell1, center1, shell2, center2, gaussiancenter)
t_expansioncoeff                =   t_expansioncoeff*factorial(expansionindex1)
t_expansioncoeff                =   t_expansioncoeff*pow(-1, expansionindex1+expansionindex3)
t_expansioncoeff                =   t_expansioncoeff*pow(gaussiancenter-atomcenter, expansionindex1-(2*expansionindex2)-(2*expansionindex3))
t_expansioncoeff                =   t_expansioncoeff*pow(epsilon, expansionindex2+expansionindex3)
t_expansioncoeff                =   t_expansioncoeff/factorial(expansionindex3, exact=True)
t_expansioncoeff                =   t_expansioncoeff/factorial(expansionindex2, exact=True)
t_expansioncoeff                =   t_expansioncoeff/factorial(expansionindex1-(2*expansionindex2)-(2*expansionindex3), exact=True)
return t_expansioncoeff

@cython.boundscheck(False)
@cython.wraparound(False)
cdef double boysfunction(int boysindex, double boysparameter):
return hyp1f1(0.5+boysindex, 1.5+boysindex, -1*boysparameter)/((2*boysindex)+1)

@cython.boundscheck(False)
@cython.wraparound(False)
def nuclearcgtos(basisobject1, basisobject2, atomobjects):
cdef double t_eniintegral   =   0.0
cdef np.ndarray center1         =   basisobject1.center
cdef np.ndarray exponents1      =   basisobject1.exponents
cdef np.ndarray coefficients1   =   basisobject1.coefficients
cdef np.ndarray shell1          =   basisobject1.shell
cdef np.ndarray normcoeffs1     =   basisobject1.normcoeffs

cdef np.ndarray center2         =   basisobject2.center
cdef np.ndarray exponents2      =   basisobject2.exponents
cdef np.ndarray coefficients2   =   basisobject2.coefficients
cdef np.ndarray shell2          =   basisobject2.shell
cdef np.ndarray normcoeffs2     =   basisobject2.normcoeffs
for index1, exponent1 in enumerate(exponents1):
for index2, exponent2 in enumerate(exponents2):
t_nuclearintegral   =   0.0
gamma               =   exponent1+exponent2
gaussiancenter, gaussianintegral    =   gaussiantheorem(center1, exponent1, center2, exponent2)
for atom in atomobjects:
for expansionindex1A in range(0, shell1[0]+shell2[0]+1):
for expansionindex1B in range(0, int(expansionindex1A//2)+1):
for expansionindex1C in range(0, int((expansionindex1A-(2*expansionindex1B))//2)+1):
auxiliary1A =   expansioncoeff2(expansionindex1A, expansionindex1B, expansionindex1C, shell1[0], center1[0], shell2[0], center2[0], atom.center[0], gaussiancenter[0], gamma)

for expansionindex2A in range(0, shell1[1]+shell2[1]+1):
for expansionindex2B in range(0, int(expansionindex2A//2)+1):
for expansionindex2C in range(0, int((expansionindex2A-(2*expansionindex2B))//2)+1):
auxiliary2A =   expansioncoeff2(expansionindex2A, expansionindex2B, expansionindex2C, shell1[1], center1[1], shell2[1], center2[1], atom.center[1], gaussiancenter[1], gamma)

for expansionindex3A in range(0, shell1[2]+shell2[2]+1):
for expansionindex3B in range(0, int(expansionindex3A//2)+1):
for expansionindex3C in range(0, int((expansionindex3A-(2*expansionindex3B))//2)+1):
auxiliary3A =   expansioncoeff2(expansionindex3A, expansionindex3B, expansionindex3C, shell1[2], center1[2], shell2[2], center2[2], atom.center[2], gaussiancenter[2], gamma)
boysindex   =   (expansionindex1A+expansionindex2A+expansionindex3A)-2*(expansionindex1B+expansionindex2B+expansionindex3B)-(expansionindex1C+expansionindex2C+expansionindex3C)
boysparam   =   dot(atom.center-gaussiancenter, atom.center-gaussiancenter)*gamma
auxiliary4A =   boysfunction(boysindex, boysparam)
t_nuclearintegral   =   t_nuclearintegral+(auxiliary1A*auxiliary2A*auxiliary3A*auxiliary4A)*atom.charge*-1.0
t_eniintegral   =   t_eniintegral+(t_nuclearintegral*coefficients1[index1]*coefficients2[index2]*normcoeffs1[index1]*normcoeffs2[index2]*gaussianintegral*(2*pi/gamma))
return t_eniintegral

@cython.boundscheck(False)
@cython.wraparound(False)
cdef double enigaussian(int shell1, int shell2, int nodes, float center1, float center2, float exponent1, float exponent2):
cdef double gamma       =   exponent1+exponent2
cdef double exponent    =   exponent1*(exponent2/gamma)

if (nodes < 0) or (nodes > (shell1+shell2)):
return 0
elif nodes == shell1 == shell2 == 0:
return exp(-1*exponent*pow(center1-center2, 2))
elif shell2 == 0:
return (1/(2*gamma))*enigaussian(shell1-1, shell2, nodes-1, center1, center2, exponent1, exponent2) - (exponent*(center1-center2)/exponent1)*enigaussian(shell1-1, shell2, nodes, center1, center2, exponent1, exponent2) + (nodes+1)*enigaussian(shell1-1, shell2, nodes+1, center1, center2, exponent1, exponent2)
else:
return (1/(2*gamma))*enigaussian(shell1, shell2-1, nodes-1, center1, center2, exponent1, exponent2) + (exponent*(center1-center2)/exponent2)*enigaussian(shell1, shell2-1, nodes, center1, center2, exponent1, exponent2) + (nodes+1)*enigaussian(shell1, shell2-1, nodes+1, center1, center2, exponent1, exponent2)

@cython.boundscheck(False)
@cython.wraparound(False)
cdef double auxiliaryhermiteeri(int nodes1, int nodes2, int nodes3, int boysorder, float exponent, float xdist, float ydist, float zdist, float dist):
cdef double boysparam   =   exponent*pow(dist, 2)
cdef double auxhermite  =   0.0

if nodes1 == nodes2 == nodes3 == 0:
auxhermite  +=  pow(-2*exponent, boysorder)*boysfunction(boysorder, boysparam)
elif nodes1 == nodes2 == 0:
if nodes3 > 1:
auxhermite  +=  (nodes3-1)*auxiliaryhermiteeri(nodes1, nodes2, nodes3-2, boysorder+1, exponent, xdist, ydist, zdist, dist)
auxhermite  +=  zdist*auxiliaryhermiteeri(nodes1, nodes2, nodes3-1, boysorder+1, exponent, xdist, ydist, zdist, dist)
elif nodes1 == 0:
if nodes2 > 1:
auxhermite  +=  (nodes2-1)*auxiliaryhermiteeri(nodes1, nodes2-2, nodes3, boysorder+1, exponent, xdist, ydist, zdist, dist)
auxhermite  +=  ydist*auxiliaryhermiteeri(nodes1, nodes2-1, nodes3, boysorder+1, exponent, xdist, ydist, zdist, dist)
else:
if nodes1 > 1:
auxhermite  +=  (nodes1-1)*auxiliaryhermiteeri(nodes1-2, nodes2, nodes3, boysorder+1, exponent, xdist, ydist, zdist, dist)
auxhermite  +=  xdist*auxiliaryhermiteeri(nodes1-1, nodes2, nodes3, boysorder+1, exponent, xdist, ydist, zdist, dist)
return auxhermite

@cython.boundscheck(False)
@cython.wraparound(False)
cdef double eripgtos(exponent1, shell1, center1, exponent2, shell2, center2, exponent3, shell3, center3, exponent4, shell4, center4):
cdef double     eripgto             =   0.0
cdef np.ndarray gaussiancenterA     =   np.zeros(3)
cdef double     gaussianintegralA   =   0.0
cdef np.ndarray gaussiancenterB     =   np.zeros(3)
cdef double     gaussianintegralB   =   0.0
gaussiancenterA, gaussianintegralA  =   gaussiantheorem(center1, exponent1, center2, exponent2)
gaussiancenterB, gaussianintegralB  =   gaussiantheorem(center3, exponent3, center4, exponent4)
cdef double xdist                   =   gaussiancenterA[0]-gaussiancenterB[0]
cdef double ydist                   =   gaussiancenterA[1]-gaussiancenterB[1]
cdef double zdist                   =   gaussiancenterA[2]-gaussiancenterB[2]
cdef double distance                =   norm(gaussiancenterA-gaussiancenterB)
cdef double combexponentA           =   exponent1+exponent2
cdef double combexponentB           =   exponent3+exponent4
cdef double combexponent            =   combexponentA*(combexponentB/(combexponentA+combexponentB))

cdef int index1                     =   0
cdef int index2                     =   0
cdef int index3                     =   0
cdef int index4                     =   0
cdef int index5                     =   0
cdef int index6                     =   0
for index1 in range(0, shell1[0]+shell2[0]+1):
for index2 in range(0, shell1[1]+shell2[1]+1):
for index3 in range(0, shell1[2]+shell2[2]+1):
for index4 in range(0, shell3[0]+shell4[0]+1):
for index5 in range(0, shell3[1]+shell4[1]+1):
for index6 in range(0, shell3[2]+shell4[2]+1):
erigaussianA    =   enigaussian(shell1[0], shell2[0], index1, center1[0], center2[0], exponent1, exponent2)
erigaussianB    =   enigaussian(shell1[1], shell2[1], index2, center1[1], center2[1], exponent1, exponent2)
erigaussianC    =   enigaussian(shell1[2], shell2[2], index3, center1[2], center2[2], exponent1, exponent2)
erigaussianD    =   enigaussian(shell3[0], shell4[0], index4, center3[0], center4[0], exponent3, exponent4)
erigaussianE    =   enigaussian(shell3[1], shell4[1], index5, center3[1], center4[1], exponent3, exponent4)
erigaussianF    =   enigaussian(shell3[2], shell4[2], index6, center3[2], center4[2], exponent3, exponent4)
erigaussianG    =   auxiliaryhermiteeri(index1+index4, index2+index5, index3+index6, 0, combexponent, xdist, ydist, zdist, distance)
result          =   erigaussianA*erigaussianB*erigaussianC*erigaussianD*erigaussianE*erigaussianF*erigaussianG*pow(-1, index4+index5+index6)
eripgto         =   eripgto+result
eripgto *=  2*pow(pi, 2.5)/(combexponentA*combexponentB*sqrt(combexponentA+combexponentB))
return eripgto

@cython.boundscheck(False)
@cython.wraparound(False)
def ericgtos(basisobject1, basisobject2, basisobject3, basisobject4):
cdef np.ndarray center1         =   basisobject1.center
cdef np.ndarray exponents1      =   basisobject1.exponents
cdef np.ndarray coefficients1   =   basisobject1.coefficients
cdef np.ndarray shell1          =   basisobject1.shell
cdef np.ndarray normcoeffs1     =   basisobject1.normcoeffs

cdef np.ndarray center2         =   basisobject2.center
cdef np.ndarray exponents2      =   basisobject2.exponents
cdef np.ndarray coefficients2   =   basisobject2.coefficients
cdef np.ndarray shell2          =   basisobject2.shell
cdef np.ndarray normcoeffs2     =   basisobject2.normcoeffs

cdef np.ndarray center3         =   basisobject3.center
cdef np.ndarray exponents3      =   basisobject3.exponents
cdef np.ndarray coefficients3   =   basisobject3.coefficients
cdef np.ndarray shell3          =   basisobject3.shell
cdef np.ndarray normcoeffs3     =   basisobject3.normcoeffs

cdef np.ndarray center4         =   basisobject4.center
cdef np.ndarray exponents4      =   basisobject4.exponents
cdef np.ndarray coefficients4   =   basisobject4.coefficients
cdef np.ndarray shell4          =   basisobject4.shell
cdef np.ndarray normcoeffs4     =   basisobject4.normcoeffs

cdef double ericgto =   0.0
for index1, exponent1 in enumerate(exponents1):
for index2, exponent2 in enumerate(exponents2):
for index3, exponent3 in enumerate(exponents3):
for index4, exponent4 in enumerate(exponents4):
result  =   eripgtos(exponent1, shell1, center1, exponent2, shell2, center2, exponent3, shell3, center3, exponent4, shell4, center4)
result  *=  coefficients1[index1]*normcoeffs1[index1]*coefficients2[index2]*normcoeffs2[index2]*coefficients3[index3]*normcoeffs3[index3]*coefficients4[index4]*normcoeffs4[index4]
ericgto =   ericgto+result
return ericgto


hartreefock.pyx

import cython   as  cython
import numpy    as  np
cimport numpy   as  np

from scipy.special  import  comb, factorial2
from libc.math      cimport exp,  pow
from numpy          import  dot,  pi, diag, amax, trace, isclose
from numpy.linalg   import  eigh, solve

@cython.boundscheck(False)
@cython.wraparound(False)
def computecorehamiltonian(np.ndarray[np.float64_t, ndim=2]kineticmat, np.ndarray[np.float64_t, ndim=2]enimat):
cdef int nbasis =   kineticmat[0].size
cdef np.ndarray[np.float64_t, ndim=2] corehamiltonian   =   np.zeros((nbasis, nbasis))
corehamiltonian =   kineticmat+enimat
return corehamiltonian

@cython.boundscheck(False)
@cython.wraparound(False)
def canonicalorthogonalization(np.ndarray[np.float64_t, ndim=2] overlapmat):
cdef int nbasis =   overlapmat[0].size
cdef np.ndarray[np.float64_t, ndim=1]   eigenvalues     =   np.zeros(nbasis)
cdef np.ndarray[np.float64_t, ndim=2]   t_eigenvalues   =   np.zeros((nbasis, nbasis))
cdef np.ndarray[np.float64_t, ndim=2]   eigenvectors    =   np.zeros((nbasis, nbasis))
cdef np.ndarray[np.float64_t, ndim=2]   orthogonalmat   =   np.zeros((nbasis, nbasis))
eigenvalues, eigenvectors   =   eigh(overlapmat)
t_eigenvalues               =   diag((eigenvalues)**-0.5)
orthogonalmat               =   dot(eigenvectors, t_eigenvalues)
return orthogonalmat

@cython.boundscheck(False)
@cython.wraparound(False)
def computehamiltonian(np.ndarray[np.float64_t, ndim=2] densitymat, np.ndarray[np.float64_t, ndim=4] erimat):
cdef int nbasis =   densitymat[0].size
cdef int index1
cdef int index2
cdef int index3
cdef int index4
cdef np.ndarray[np.float64_t, ndim=2]   hamiltonian =   np.zeros((nbasis, nbasis))
for index1 in range(0, nbasis):
for index2 in range(0, nbasis):
for index3 in range(0, nbasis):
for index4 in range(0, nbasis):
hamiltonian[index1, index2] =   hamiltonian[index1, index2]+(densitymat[index3, index4]*(erimat[index1, index2, index3, index4]-0.5*(erimat[index1, index4, index3, index2])))
return hamiltonian

@cython.boundscheck(False)
@cython.wraparound(False)
def computedensity(int nelectrons, np.ndarray[np.float64_t, ndim=2] overlapmat):
cdef int nbasis =   overlapmat[0].shape[0]
cdef np.ndarray[np.float64_t, ndim=2]   densitymat  =   np.zeros((nbasis, nbasis))
cdef int index1
cdef int index2
cdef int index3
for index1 in range(0, nbasis):
for index2 in range(0, nbasis):
for index3 in range(0, int(nelectrons//2)):
densitymat[index1, index2]  +=  2.0*overlapmat[index1, index3]*overlapmat[index2, index3]
return densitymat

@cython.boundscheck(False)
@cython.wraparound(False)
def diisorthogonalization(np.ndarray[np.float64_t, ndim=2] overlapmat):
eigenvalues, eigenvectors   =   eigh(overlapmat)
halfeigmat                  =   diag((eigenvalues)**-0.5)
rightmat                    =   dot(eigenvectors, halfeigmat)
orthomat                    =   dot(rightmat, eigenvectors.T)
return orthomat

@cython.boundscheck(False)
def diisEngine(molecule, np.ndarray[np.float64_t, ndim=2] fock, np.ndarray[np.float64_t, ndim=2] density):
cdef int dimensions         =   fock[0].size
eigenvalues, eigenvectors   =   eigh(molecule.overlapmat)
halfeigmat                  =   diag((eigenvalues)**-0.5)
rightmat                    =   dot(eigenvectors, halfeigmat)
orthomat                    =   dot(rightmat, eigenvectors.T)
t_errorvector               =   dot(fock, dot(density, molecule.overlapmat))-dot(molecule.overlapmat, dot(density, fock))
leftmat                     =   dot(orthomat.T, t_errorvector)
errorvector                 =   dot(leftmat, orthomat)
molecule.diisfock.append(fock)
molecule.diiserror.append(errorvector)
cdef int fockdimension               =   len(molecule.diisfock)
if fockdimension > 6:
del molecule.diisfock[0]
del molecule.diiserror[0]
fockdimension   =   fockdimension-1
cdef np.ndarray[np.float64_t, ndim=2] bmatrix   =   np.zeros((fockdimension+1, fockdimension+1))
bmatrix[-1,:]               =   -1.0
bmatrix[:,-1]               =   -1.0
bmatrix[-1,-1]              =   0.0
for index1 in range(0, fockdimension):
for index2 in range(0, index1+1):
bmatrix[index1, index2] =   trace(dot(molecule.diiserror[index1].T, molecule.diiserror[index2]))
bmatrix[index2, index1] =   trace(dot(molecule.diiserror[index1].T, molecule.diiserror[index2]))
cdef np.ndarray[np.float64_t, ndim=1] residue    =   np.zeros(fockdimension+1)
residue[-1]                 =   -1.0
cdef np.ndarray[np.float64_t, ndim=1] weights    =   solve(bmatrix, residue)
assert isclose(sum(weights[:-1]),1.0)
cdef np.ndarray[np.float64_t, ndim=2] newfock    =   np.zeros((dimensions, dimensions))
for index, fockelement in enumerate(molecule.diisfock):
newfock =   newfock+(weights[index]*fockelement)
return newfock

@cython.boundscheck(False)
@cython.wraparound(False)
def restrictedhatreefock(int nelectrons, np.ndarray[np.float64_t, ndim=2] corehamiltonian, np.ndarray[np.float64_t, ndim=2] orthogonalmat, np.ndarray[np.float64_t, ndim=2] densitymat, np.ndarray[np.float64_t, ndim=4] erimat):
cdef int nbasis =   corehamiltonian[0].size
cdef np.ndarray[np.float64_t, ndim=2]   hamiltonianmat  =   computehamiltonian(densitymat, erimat)
cdef np.ndarray[np.float64_t, ndim=2]   fockmatrix      =   hamiltonianmat+corehamiltonian
cdef np.ndarray[np.float64_t, ndim=2]   orthofockmatrix =   dot(orthogonalmat.conj().T, dot(fockmatrix, orthogonalmat))
cdef np.ndarray[np.float64_t, ndim=1]   orbitalenergies =   np.zeros(nbasis)
cdef np.ndarray[np.float64_t, ndim=2]   orthocoeffs     =   np.zeros((nbasis, nbasis))
orbitalenergies, orthocoeffs                            =   eigh(orthofockmatrix)
cdef np.ndarray[np.float64_t, ndim=2]   canonicalcoeffs =   dot(orthogonalmat, orthocoeffs)
cdef np.ndarray[np.float64_t, ndim=2]   t_densitymat    =   computedensity(nelectrons, canonicalcoeffs)
cdef double maxdensitydiff
cdef double rmsdensitydiff
maxdensitydiff, rmsdensitydiff                          =   checkconvergence(densitymat, t_densitymat)
return orthofockmatrix, fockmatrix, hamiltonianmat, orbitalenergies, canonicalcoeffs, t_densitymat, maxdensitydiff, rmsdensitydiff

@cython.boundscheck(False)
@cython.wraparound(False)
def restrictedhatreefockDIIS(molecule, int nelectrons, np.ndarray[np.float64_t, ndim=2] corehamiltonian, np.ndarray[np.float64_t, ndim=2] orthogonalmat, np.ndarray[np.float64_t, ndim=2] densitymat, np.ndarray[np.float64_t, ndim=4] erimat):
cdef int nbasis =   corehamiltonian[0].size
cdef np.ndarray[np.float64_t, ndim=2]   hamiltonianmat  =   computehamiltonian(densitymat, erimat)
cdef np.ndarray[np.float64_t, ndim=2]   fockmatrix      =   hamiltonianmat+corehamiltonian
cdef np.ndarray[np.float64_t, ndim=2]   diisfock        =   diisEngine(molecule, fockmatrix, densitymat)
cdef np.ndarray[np.float64_t, ndim=2]   orthofockmatrix =   dot(orthogonalmat.conj().T, dot(diisfock, orthogonalmat))
cdef np.ndarray[np.float64_t, ndim=1]   orbitalenergies =   np.zeros(nbasis)
cdef np.ndarray[np.float64_t, ndim=2]   orthocoeffs     =   np.zeros((nbasis, nbasis))
orbitalenergies, orthocoeffs                            =   eigh(orthofockmatrix)
cdef np.ndarray[np.float64_t, ndim=2]   canonicalcoeffs =   dot(orthogonalmat, orthocoeffs)
cdef np.ndarray[np.float64_t, ndim=2]   t_densitymat    =   computedensity(nelectrons, canonicalcoeffs)
cdef double maxdensitydiff
cdef double rmsdensitydiff
maxdensitydiff, rmsdensitydiff                          =   checkconvergence(densitymat, t_densitymat)
return orthofockmatrix, fockmatrix, hamiltonianmat, orbitalenergies, canonicalcoeffs, t_densitymat, maxdensitydiff, rmsdensitydiff

@cython.boundscheck(False)
@cython.wraparound(False)
def checkconvergence(np.ndarray[np.float64_t, ndim=2] densitymat, np.ndarray[np.float64_t, ndim=2] t_densitymat):
cdef np.ndarray[np.float64_t, ndim=2]   densitydiff     =   t_densitymat-densitymat
cdef double maxdensitydiff  =   amax(densitydiff)
cdef double rmsdensitydiff  =   pow(sum(sum((densitydiff)**2))/4.0, 0.5)
return maxdensitydiff, rmsdensitydiff

@cython.boundscheck(False)
@cython.wraparound(False)
def totalenergy(float nucenergy, np.ndarray[np.float64_t, ndim=2] fockmatrix, np.ndarray[np.float64_t, ndim=2] densitymatrix, np.ndarray[np.float64_t, ndim=2] corehamiltonian):
cdef double t_energy    =   0.0
cdef int nbasis         =   fockmatrix[0].size
for index1 in range(0, nbasis):
for index2 in range(0, nbasis):
t_energy    =   t_energy+(0.5*densitymatrix[index1,index2]*(fockmatrix[index1,index2]+corehamiltonian[index1,index2]))
return t_energy+nucenergy
$$$$

• Did you consider using Numba? Commented Jun 21, 2023 at 13:59
• I did not. I was under the impression that static compilation was a better way to go. Because many programs have already implemented routines using Cython; PyQuante for example. Commented Jun 21, 2023 at 14:37

Use snake_case for naming, it's a standard practice in Python. nuclear_energy instead of nuclearenergy, get_atom instead of getatom, etc.. Much easier to read.

Don't use camelCase and, more importantly, don't mix naming conventions together:

periodictable   =   safe_load(periodic)
periodicTable   =   {y: x[y] for x in periodictable for y in x}


At first I thought you're assigning periodictable twice.

Don't pollute the global scope by creating names here. The piece above creates 2 names periodictable and periodicTable that will be accessible from anywhere in the code. Instead create a PeriodicTable class with a get_atom member method.

Let's look at this function:

def nuclearenergy(atomobjects):
atompairs       =   list(combinations(atomobjects,2))
nuclearenergy   =   0.0
for atompair in atompairs:
center1         =   atompair[0].center
charge1         =   atompair[0].charge
center2         =   atompair[1].center
charge2         =   atompair[1].charge
distance        =   norm(center1-center2)
energy          =   charge1*charge2/distance
nuclearenergy   =   nuclearenergy+energy
return nuclearenergy


It can be rewritten as:

def get_nuclear_energy(atoms):
return sum([pair[0].charge * pair[1].charge / norm(pair[0].center - pair[1].center) for pair in combinations(atoms, 2)])


This is shorter and more readable while also being about 6% faster (according to a simple test via PyCharm profiling). The only information we lost is that norm(center1-center2) is distance, but I assume people who are going to contribute to this code are already familiar with the formulas.

Now this function:

def overlap(basisobjects):
overlaps    =   []
for basisobjecttuple in basisobjects:
result  =   overlapcgtos(basisobjecttuple[0], basisobjecttuple[1])
overlaps.append((result, basisobjecttuple[0].basisindex, basisobjecttuple[1].basisindex))
return overlaps


We see a similar structure: create an empty container-accumulator and fill it with stuff in a cycle. This is verbose and slow, since you have to call append() every time you're adding to the container. It can be rewritten like this:

def get_overlaps(bases):
return [(overlapcgtos(basis[0], basis[1]), basis[0].basis_index, basis[1].basis_index) for basis in bases]


Can't check the performance without compiling pyx, but should be about 44% faster (tested append in a cycle vs list comprehension via PyCharm profiler).

The same goes for the following 3 functions.

def integralEngine ...

def scfEngine ...


Can't tell what those functions are for by reading their name. Functions should be named with verbs that represent the action the function is doing, e.g. start_scf_engine.

if calctype == "overlap":
...
if calctype == "kinetic":
...
if calctype == "eni":
...
if calctype == "eri":
...


Consider defining an enum for that, so things don't break when you decide to change the string in one place and forget to change it in another. Also use elif to get rid of unnecessary checks.

pow(2, 2*totalangmomentum)*pow(2, 1.5)


Python has a pow operator: **

2 ** (2 * totalangmomentum) * 2 ** 2.5


[x for x in molecule.geometry]


This statement does nothing with molecule.geometry. If you wanted to copy it, there is a function for that in the copy package: copy.copy(molecule.geometry)

Consider using if __name__ == '__main__'` construction to run the script.

In general, try running your code with a profiler and see what parts are the slowest. Try to speed them up first and don't bother with stuff that is already fast. And consider trying Numba, as @greybeard suggested in the comments.

• @greybeard, no. Is there a point? Commented Jun 24, 2023 at 16:10
• @greybeard haven't noticed, ma bad Commented Jun 24, 2023 at 19:43