Restricted Hartree-Fock with STO-3G From Scrach. How to Use

Introduction
Programs
1. Program structure
How to Use
1. Set molecular arrangement and drowing color list of atoms
2. Options
Running
Last

Introduction

In quantum chemistry, Hartree-Fock method is introduced as a method to handle a multiple particle system and calculate energys of them.

However, how to compute them in acutual programs is not explained in detail in most testbooks.

So, I was interested in how it works and wanted to understand.

That’s why I tried to code RHF/STO-3G from scrach.

Here, I introduce how to use my code.

Programs

I’ve uploaded the code of RHF/STO-3G to Github.

You can try the calculation in Python.

GitHub - k7632670-hash/quantum-chemistry: Hartree-Fock STO-3G program. Moleclar Energy, Structure Optimization, Vibration Analysis

Hartree-Fock STO-3G program. Moleclar Energy, Structure Optimization, Vibration Analysis - k7632670-hash/quantum-chemistry

Program structure

My code is composed of 6 programs.

main.py
basis.py
properties.py
integrals.py
scfdriver.py
drawing.py

main.py : You need to edit it. Main procedure.

basis.py: The list of gaussian basis set

properties.py: The list of atomic information

integrals.py: Integral methods of overlap, nuclear attraction, kinetic energy, electron repulsive integral

scfdriver.py: Self-consistent field method

drawing.py: Draw moleclar you set and molecular orbitals you computed

How to Use

A part of main.py that you need to edit is shown below:

Set molecular arrangement and drowing color list of atoms

This is H2O case.

#----------#
### Main ###
#----------#

#Set molecular arrangement. unit:Bohr
mole_ar = [
            ["H",[ 0, 1.430523, 1.430523]],
            ["H",[ 0,-1.430523, 1.430523]],
            ["O",[ 0,    0, 0 ]],
          ]

#Set colors to draw molecular
color_set = {
                "H":"blue",
                "O":"red",
                "C":"grey"
            }

Molecular arrangement is set in “mole_ar” as list type.

You type atomic symbols and coordinates.

Note that: Coordinates are written in cartesian coordinate. And Unit is Bohr.

This program is applicable to atoms from H to Ne in periodic table because gusian besis sets corrsponding to other atoms are not prepared yet.

“color_set” is not always necessary. If DRAW_MOLE is True, you need to edit it.

You also set atomic symbols and colors.

Drawing code is based on matplotlib and color is assigned by “color =” or “c=”.

You can use color names and RGBa. Please refer matplotlib web site in detail.

List of named colors — Matplotlib 3.11.1 documentation

Options

There are several options you can set.

### Set Parameters ###
DRAW_MOLE = True       #True: Enable drawing molecular
DRAW_MOLE_ORBs = True  #True: Enalbe drawing molecular orbitals

tole = 1e-8     #SCF Convergence tolerance
max_ite = 200   #SCF Max Iteration
alpha = 0.7     #SCF Mixing parameter for updating P

DRAW_MOLE: Drow molecular you set in more_ar

Draw_MOLE_ORBs: Draw molecular orbirtals calculated after SCF procedure.

tole: SCF convergence tolerance.

max_ite: Max iteration of SCF procedure

alpha: Mixing parameter to update P matrix in SCF procedure

Note that: In this code, initial C matrix is defined 0 matrix. It usually works. But, sometimes it yeilds wrong results. Then, you need to define goody initial C matric by yourself.

You find lines shown below in main.py and edit it.

#Set initial C matrix for SCF 
Cocc = np.zeros((nao,nocc))

Running

I show results of H2O case.

### Final Results ###
Molecler Arangement:
[['H', [0, 1.430523, 1.430523]], ['H', [0, -1.430523, 1.430523]], ['O', [0, 0, 0]]]
eg: [-20.272  -1.231  -0.538  -0.471  -0.395] hartree
eu: [0.521 0.585] hartree
Eelectron: -83.20874745930698 hartree
Etotal: -74.95 hartree
End:1.203709363937378 s

eg: orbital energies for occupied orbitals.

eu: orbital energies for unoccupied orbitals.

Eelectron: electronic energy.

Etotal: total energy including nuclear repulsion energy.

Here is results by Veloxchem.

               Spin-Restricted Hartree-Fock:                                                                              
               -----------------------------                                                                              
               Total Energy                       :      -74.9504116927 a.u.                                              
               Electronic Energy                  :      -83.2078170289 a.u.                                              
               Nuclear Repulsion Energy           :        8.2574053361 a.u.                                              
               ------------------------------------

Good agreement! You can see if this code works.

And other outputs are shown below:

Molecular Plot: H2O

Molecular orbitals 2D plot in XY, XZ and YZ planes

Last

This page introduced how to use RHF/STO-3G briefly.

I’m planning to expand my code to structure optimization, normal vibration analysis, IR intensity, etc. Actually, those methods are built already. But it’s not sorted well.

In addition, I’ll explain how RHF/STO-3G works, particulally gaussian integral and obara-saika recursion and so on. It was very educated to understand it.

Please forgive me if there is wrong with my code. I’ll fix it if you let me know about it.

Thank you.