ifilot/hfcxx

Hartree-Fock C++ code

HFCXX

Overview

HFCXX is a Hartree-Fock code in C++. The program currently supports only single-point calculations (i.e. no geometry optimizations). The basis sets STO-3G and STO-6G are included for atoms up to N=8 (oxygen).

Note: This program is mainly made for educational purposes. There are some (obvious) things to improve as suggested by others (e.g. loading of basis set coefficients from external files). My aim is however to keep this code simple and not use too many additional features.

• Interested in learning Density Functional Theory? Have a look at DFTCXX.
• Want to know more about electronic structure calculations, have a look at my free lecture book.

Compilation

Ensure you have all the required packages installed

sudo apt install build-essential cmake libcppunit-dev

mkdir build && cd build
cmake ../src
make -j9

To test the compilation, run

Note: The testsuite contains a HF calculation of benzene using the STO-6g basis set. Since HFCXX is not optimized for speed, running the tests can take a few minutes to complete.

make test

A successful test should give an output similar to the following

Running tests...
Test project /mnt/d/PROGRAMMING/CPP/hfcxx/build
    Start 1: TestMolecules
1/1 Test #1: TestMolecules ....................   Passed   50.70 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =  50.72 sec

Usage

./hfcxx _input_file_ > _output_file_

For example, to run a HF calculation for H2, run

./hfcxx molecules/h2.in > h2.out

hfcxx's Issues

hf.cpp line 199

Hi,

I am implementing a Hartree-Fock program for my blog. I pulled your program to compare results for electron-electron integrals (I am using a different method of calculating them). Of course I still have issues in the implementation, but anyway, I was looking over your code and noticed:

Fp = trimatprod(Xp,F,X);

but then on line 199:

Eigsym es(F,Cc,molorben);

It looks to me as what you intended to have is:

Eigsym es(Fp,Cc,molorben);

I don't know exactly what Eigsym does, I didn't look into it. I suppose it only solves the ordinary eigenvalue problem (that is, not the generalized one), since you don't pass to it the overlap matrix. I think you should pass the transformed Fock operator to it.

I changed the code and executed it for H2O (not your input file, but a configuration I got from a Mathematica Journal article, which is a little bit different) and the results seem to come out a little better.

Thanks,
Adrian

Problem with diatomic molecules where both atoms have orbitals different than s

Hi,

I think I found a problem with such setup. For example for N2 molecule, with the following configuration:

basis = sto-3g
charge = 0
atoms = 2
units = bohr
N 0.0000 0.0000 0.0000
N 2.147932 0.0000 0.0000

the given result is:

Total energy [Hartree]: -144.1995252

This is with the fix of the bug I previously signaled.

According to some paper I was looking into, the Hartree-Fock limit for N2 is -108.994.

I am comparing your results with a program I've written. It has different basis sets than yours but the values are quite similar. I too have issues with diatomic molecules for this kind of configuration. It works nicely for example for H2O or any diatomic molecule I tried that has a hydrogen in it, but fails once I put two atoms with both more than s orbitals.

I'm using recurrence relations (I think they are called Hamilton & Schaefer) to solve the integrals, for this described case I could spot different terms in both the nuclear matrix and in electron-electron integrals than given by your program. Most of the values are basically identical, but some come out differently.

It doesn't mean yours are wrong since I am also getting bad results for such molecules, but it seems that I'm getting results closer to the expected value, for example for the above molecule I get -113.677. Still wrong :(

Thanks,
Adrian

Speeding up factorial and double factorial

Hi,

May I suggest an improvement?

There are ways to precompute the values and use the already computed values instead of computing them each time, using tables of values (computed even at compile time). This should speed up computations a little. Here is how I did it:

https://github.com/aromanro/HartreeFock/blob/master/HartreeFock/MathUtils.h
https://github.com/aromanro/HartreeFock/blob/master/HartreeFock/MathUtils.cpp

In my case it doesn't make much difference but you use different methods for computing the integrals where it might count more.

Processor cycles are wasted in calling Fgamma with identical parameters

In coulomb_repulsion and nuclear, there are calls like Fgamma(i + j + k, ...) within loops.

This is bad for performance. You call Fgamma Imax * Jmax * Kmax times, but only Imax + Jmax + Kmax invocations will have distinct parameters.

This can be fixed by memorizing the results & precomputing:

#define FGAMMA_CACHE_FOR(_size, expr)                                                                                                                \
    std::vector<double> fgamma_lookup_table((_size) + 1);                                                                                            \
    for (size_t i = 0; i < fgamma_lookup_table.size(); i++)                                                                                          \
    {                                                                                                                                                \
        fgamma_lookup_table[i] = Fgamma(i, expr);                                                                                                    \
    }

...
        FGAMMA_CACHE_FOR(l1 + l2 + m1 + m2 + n1 + n2, rcp2 * gamma);

        for (int i = 0; i <= l1 + l2; i++)
        {
            for (int j = 0; j <= m1 + m2; j++)
            {
                for (int k = 0; k <= n1 + n2; k++)
                {
                    sum += ax[i] * ay[j] * az[k] * fgamma_lookup_table[i + j + k];
                }
            }
        }
...
        FGAMMA_CACHE_FOR(la + lb + lc + ld + ma + mb + mc + md + na + nb + nc + nd, 0.25 * rpq2 / delta);

        for (int i = 0; i <= (la + lb + lc + ld); i++)
        {
            for (int j = 0; j <= (ma + mb + mc + md); j++)
            {
                for (int k = 0; k <= (na + nb + nc + nd); k++)
                {
                    sum += bx[i] * by[j] * bz[k] * fgamma_lookup_table[i + j + k];
                }
            }
        }

benzene gives incorrect SCF energy for sto-6g basis set

benzene gives incorrect SCF energy for sto-6g basis set

Hello,

Just to let you know. I didn't look into it. I use a different method for the integrals. From memory way back I think I may have spotted something in the kinetic integrals but that was for d orbitals anyway, which you don't use, but don't hold me to that.

Your project was a great help to get started and inspired me to write my own.

hfcxx

Number of orbitals: 36
Number of GTOs: 216
Number of iterations: 40
Total energy [Hartree]: -228.9254378

My code

Basis set = STO-6G
Basis set type = cartesian
Number of basis primitives = 144
Number of actual GTOs = 216
Number of basis functions = 36
Number of shells = 24

** more stuff ****

E(NUC) / Eh = 203.36795083
E(1E ) / Eh = -714.19179571
E(2E ) / Eh = 280.69350348
E(SCF) / Eh = -230.13034140

Which agrees with the psi4 pacakge.

Your sto-3g result checks out fine.

Cheers.