Molecular Geometry Optimization using Direct Restricted Hartree - Fock<center>

Larry Tesler<center>

CHM6470<center>

Fall 2014<center>

LTPsi - Basic Quantum Chemistry Suite in Python

Current Features:

Restricted Hartree Fock (RHF)
- DIIS Convergence Acceleration

Several Basis Sets
- STO-3G, STO-6G, 3-21G, 6-31G, 6-31G, 6-31++G

Fast Two-Electron Integrals using Rys Polynomials
- Implemented in Cython for speed!!!
- Numerical Rys Quadrature
- Chebyshev-Gauss Quadrature for f orbitals

Dipole Moment Calculation

Gradient Descent Optimization

Mulliken Population Analysis

Møller - Plesset Pertubation Theory (MP2) Energy Correction

Molecule Rendering with OpenGL

Molecular orbital visulization with Mayavi

Can clone from GitHub: https://github.com/LT12/LTPsi

References

Restricted Hartree - Fock / MP2 / DIIS / Molecular Properties:

Crawdad Chemistry Programming Tutorial
- http://sirius.chem.vt.edu/wiki/doku.php?id=crawdad:programming

Modern Quantum Chemistry by Szabo and Ostlund

One - Electron Integrals:

Overlap Integral: http://www.mathematica-journal.com/2012/02/evaluation-of-gaussian-molecular-integrals/

Kinetic Energy Integral: H. Taketa, S. Huzinaga, K. O-ohata, Gaussian-Expansion Methods for Molecular Integrals, J. Phys. Soc. Japan, 21 (1966), p. 2312

Nuclear - Electron Attraction Integral: http://www.mathematica-journal.com/2014/09/evaluation-of-gaussian-molecular-integrals-3/

Dipole Integral: http://www.chem.helsinki.fi/~manninen/aqc2012/Session180412.pdf

Two - Electron Integrals:

Flocke N. , Lotrich V.,Efficient electronic integrals and their generalized derivatives for object oriented implementations of electronic structure calculations.,J. Comput Chem. 2008 Dec;29(16):2722-36. doi: 10.1002/jcc.21018.

Step 1 - Calculate Nuclear - Nuclear Repulsion Energy:

$$ V_{NN} = \sum_{i}^{N}\sum_{j}^{i-1}\frac{Z_{i}Z_{j}}{|\vec{r_{i}}-\vec{r_{j}}|}$$



In [ ]:

    
def calcNuclRepulsion(self):

    self.NuclRepEnergy = 0
    Z = getAtomicCharge

    if self.numAtom > 1:
        for i in xrange(self.numAtom):
            for j in xrange(i):

                R_ij =  self.cartMatrix[i,:] - self.cartMatrix[j,:]
                nR_ij = np.sqrt(np.dot(R_ij,R_ij))

                Z_iZ_j = Z(self.atomType[i]) * Z(self.atomType[j])

                self.NuclRepEnergy += Z_iZ_j/nR_ij

Step 2 - Calculate Molecular Integrals:

$$S_{ij}, \,\,\, H^{\text{Core}}_{ij},\,\,\, (ij|kl)$$



In [ ]:

    
Ints = MolecularIntegrals(self.orbList, self.atomType, self.cartMatrix, 1)

self.CoreHam, self.OverlapMatrix, self.ERT = Ints.CoreHam, Ints.OLMatrix, Ints.ERT

Step 3 - Calculate Orthogonalization Matrix:

We need to solve Roothaan-Hall Equations:

$$\mathbf{FC}=\mathbf{SC\epsilon}$$

In an orthonormal AO basis, the overlap matrix is the identity matrix, hence the Roothan Equations become:

$$\mathbf{F'C'}=\mathbf{C'\epsilon}$$

In order to convert the Fock matrix into an orthonormal basis set, we must calculate the orthogonalization matrix as follows:

$$ S^{-\frac{1}{2}}=\mathbf{V\Lambda ^ {-\frac{1}{2}}V^{T}}$$



In [ ]:

    
def symOrthogMatrix(self):

    eig, eigVM = eigh(self.OverlapMatrix,turbo = True)

    eigT = np.diag(np.reciprocal(np.sqrt(eig)))

    self.symS = np.dot(eigVM, np.dot(eigT, eigVM.T))

Step 4 - Transform Fock Matrix into Orthonormal Basis and Calculate AO Coefficients

Transform Fock matrix into orthonormal AO basis

$$ F' = S^{-\frac{1}{2}^{T}} F S^{-\frac{1}{2}} $$

Solve Roothaan Equation:

$$ \mathbf{F'C'}=\mathbf{C'}\epsilon$$

Convert AO coefficents back to non-orthonormal AO basis:

$$ C = S^{-\frac{1}{2}}C'$$



In [ ]:

    
def transFock(self, Fock):

    self.transFockM = np.dot(self.symS.T,np.dot(Fock,self.symS))

    self.eig,eigVM = eigh(self.transFockM,turbo=True)

    self.AOcoefs = np.dot(self.symS,eigVM)

Step 5 - Calculate Density Matrix

$$D_{ij} = \sum_{\phi}^{\text{Occ. }\phi}C_{i}^{\phi}C_{j}^{\phi}$$



In [ ]:

    
def compDensityMatrix(self):
    self.densityM = np.einsum('im ,jm ->ij', self.AOcoefs[:,:self.nOcc],
                                             self.AOcoefs[:,:self.nOcc])

Step 6 - Calculate SCF Electronic Energy

$$E_{e^{-}} = \sum_{i}^{N}\sum_{j}^{N}D_{ij}(H^{\text{core}}_{ij}+F_{ij})$$

-In order to utilize Numpy's einsum function, sum will be distributed:

$$E_{e^{-}} = \sum_{i}^{N}\sum_{j}^{N}D_{ij}H^{\text{core}}_{ij}+\sum_{i}^{N}\sum_{j}^{N}D_{ij}F_{ij}$$



In [ ]:

    
def compElectronicEnergy(self, Fock):
    self.ElectrEnergy = np.einsum('ij,ij', self.densityM,
                                           self.CoreHam + Fock)

Step 7 - Calculate New Fock Matrix

$$ F = H^{\text{core}} + \sum_{k}^N\sum_{l}^N D_{kl}(\,2(ij|kl)-(ik|jl)\,)$$



In [ ]:

    
def compFock(self):

    self.FockM = self.CoreHam +\
                 2 * np.einsum('kl,ijkl->ij',self.densityM,self.ERT)\
                  - np.einsum('kl,ikjl->ij',self.densityM,self.ERT)

    #Storing Fock Matricies for DIIS    
    self._fockList[self._FockCount % 6] = np.copy(self.FockM)
    self._FockCount += 1
    if self.itercount>2:
        self._DIIS_switch = True
        self._Fock_switch = True

Iterate Until Convergence is reached:



In [ ]:

    
def SCF(self):

    self.symOrthogMatrix()

    #initial density guess
    if self.densityM is None:
        self.transFock(self.CoreHam)
        self.compDensityMatrix()

    self.compElectronicEnergy(self.CoreHam)

    self.itercount, prevElecEnergy, self._DIIS_switch, self._FockCount = 0, 0, False, 0
    self._fockList,self._errorList = np.ndarray(6,dtype=np.ndarray),np.ndarray(6,dtype=np.ndarray)
    self._Fock_switch = True

    while True:

        if abs(self.ElectrEnergy-prevElecEnergy) < 1E-12: break
        prevElecEnergy = self.ElectrEnergy

        if self.itercount == 2:
            self._DIIS_switch = True
            self._Fock_switch = False

        if self._DIIS_switch == True:
            self.DIIS()
        else:
            self.compFock()

        if self._Fock_switch == True:
            self.compErrorMatrix()

        self.transFock(self.FockM)
        self.compDensityMatrix()
        self.compElectronicEnergy(self.FockM)
        self.itercount += 1

    self.SCFEnergy = self.ElectrEnergy + self.NuclRepEnergy
    self.TotalEnergy = self.SCFEnergy

MP2

First, we need to transform two-electron integrals from AO basis into MO basis
- This can be done naively ($O(N^{8})$) as follows:
  
  $$ (pq|rs) = \sum_{i}\sum_{j}\sum_{k}\sum_{l}C_{i}^{p}C_{j}^{q}C_{k}^{r}C_{l}^{s}(ij|kl)$$
- A smarter way ($O(N^{5})$) would be to use partial quarter transformations:
  
  $$ (pq|rs) = \sum_{i}C_{i}^{p}\left[\sum_{j}C_{j}^{q}\left[\sum_{k}C_{k}^{r}\left[\sum_{l}C_{l}^{s}(ij|kl)\right]\right]\right]$$



In [ ]:

    
#Convert Electron Repulsion Tensor from AO basis to MO basis
MOERT = np.einsum('ijkl, ls -> ijks',MOERT,C)
MOERT = np.einsum('ijks, kr -> ijrs',MOERT,C)
MOERT = np.einsum('ijrs, jq -> iqrs',MOERT,C)
MOERT = np.einsum('iqrs, ip -> pqrs',MOERT,C)

Now the MP2 energy correction can be calculated:

$$E_{\text{MP2}} = \sum_{i}^{\text{Occ. }\phi}\sum_{j}^{\text{Occ. }\phi}\sum_{a}^{\text{UnOcc. }\phi}\sum_{b}^{\text{UnOcc. }\phi}\frac{(ij|ab)(2(ij|ab)-(ib|ja))}{(\epsilon_{i}+\epsilon_{i})-(\epsilon_{a}+\epsilon_{b})}$$



In [ ]:

    
for i,j,a,b in itertools.product(xrange(nOcc),xrange(nOcc),
                                 xrange(nOcc,N),xrange(nOcc,N)):

    self.MP2Energy += MOERT[i,a,j,b] *\
                      (2 * MOERT[i,a,j,b] - MOERT[i,b,j,a])/\
                      (OrbEnergy[i] + OrbEnergy[j]
                       - OrbEnergy[a] - OrbEnergy[b])



In [ ]:

    
def MP2(self):

    #define N-number of orbitals, nOcc-number of occupied orbitals,
    #MOERT-ERT in MO basis, Orbenergy- energy of each orbital in
    #ordered manner, C-Matrix of AO coefficients
    N,nOcc,MOERT,OrbEnergy,C = len(self.orbList),self.nOcc,\
                               np.copy(self.ERT),self.eig,\
                               self.AOcoefs

    #Convert Electron Repulsion Tensor from AO basis to MO basis
    MOERT = np.einsum('ijkl, ls -> ijks',MOERT,C)
    MOERT = np.einsum('ijks, kr -> ijrs',MOERT,C)
    MOERT = np.einsum('ijrs, jq -> iqrs',MOERT,C)
    MOERT = np.einsum('iqrs, ip -> pqrs',MOERT,C)

    self.MP2Energy = 0

    for i,j,a,b in itertools.product(xrange(nOcc),xrange(nOcc),
                                     xrange(nOcc,N),xrange(nOcc,N)):

        self.MP2Energy += MOERT[i,a,j,b] *\
                          (2 * MOERT[i,a,j,b] - MOERT[i,b,j,a])/\
                          (OrbEnergy[i] + OrbEnergy[j]
                           - OrbEnergy[a] - OrbEnergy[b])

    self.TotalEnergy += self.MP2Energy