Is there a simple free working code that implements Hartree–Fock or DFT?

Question

Is there somewhere in the internet a good simple free working code (with a minimal number of lines in its presentation) that implements a calculation of eigenvalues with Hartree–Fock (or DFT) algorithm for $\ce{H2}$, $\ce{He}$, $\ce{Li}$, $\ce{H2O}$, $\ce{CO2}$ (or something else, something simple) on Python, Octave, C/C++, Fortran or pseudocode (no matter), maybe an article with a description of the code? In other words, writing a code for a concrete simple molecule (e.g. $\ce{H2O}$), what is a number of lines of the code which could be spent in minimum? (Naturally, not filling one line up to the maximum, but transferring the code to a next line within the meaning.)

Answer 1

SlowQuant

There's dozens of free, working, open source codes for Hartree-Fock and DFT available online. See Susi Lehtola's review paper about open source quantum chemistry software. I'm not sure which one is the shortest or simplest, so maybe we can play a bit of code golf here.

One program does come to my mind though, which is SlowQuant because the author put an emphasis on readability rather than speed (which is unusual for quantum chemistry software, where speed is usually or at least traditionally was the priority):

SlowQuant is a molecular quantum chemistry program written in python. Even the computational demanding parts are written in python, so it lacks speed, thus the name SlowQuant.

The author of SlowQuant is also active here at MMSE which is an added bonus if you have any questions about it here. On the other hand the GitHub page suggests that the code hasn't been actively developed in the last 5 years, except for some small changes that were made slightly more than 12 months ago.

The Hartree-Fock part is in this folder. The HartreeFock.py file contains 129 lines but more than 29 lines are either whitespace or comments, so it's actually under 100 lines of "code" (the UHF.py file is only a bit longer, whereas the uhf.py file in PySCF contains more than 1000 lines!):

import numpy as np
import scipy.linalg
from slowquant.hartreefock import DIIS

def diagonlize(M):
    eigVal, eigVec = np.linalg.eigh(M)
    return eigVal, eigVec

def symm_orth(eigVal, eigVec):
    diaVal = np.diag(eigVal)
    M = np.dot(np.dot(eigVec,scipy.linalg.sqrtm(np.linalg.inv(np.diag(eigVal)))),np.matrix.transpose(eigVec))
    return M

def HartreeFock(input, VNN, Te, S, VeN, Vee, deTHR=10**-6,rmsTHR=10**-6,Maxiter=100, DO_DIIS='Yes', DIIS_steps=6, print_SCF='Yes'):
    # ###############################
    #
    # VNN = nuclear repulsion
    # Te  = kinetic energy
    # S   = overlap integrals
    # VeN = nuclear attraction
    # Vee = two electron integrals
    # 
    # ###############################
    
    #Core Hamiltonian
    Hcore = VeN+Te
    
    #Diagonalizing overlap matrix
    Lambda_S, L_S = diagonlize(S)
    #Symmetric orthogonal inverse overlap matrix
    S_sqrt = symm_orth(Lambda_S, L_S)
    
    #Initial Density
    F0prime = np.dot(np.dot(np.matrix.transpose(S_sqrt),Hcore),np.matrix.transpose(S_sqrt))
    eps0, C0prime = diagonlize(F0prime)
    
    C0 = np.matrix.transpose(np.dot(S_sqrt, C0prime))
    
    #Only using occupied MOs
    C0 = C0[0:int(input[0,0]/2)]
    C0T = np.matrix.transpose(C0)
    D0 = np.dot(C0T, C0)        
    
    # Initial Energy
    E0el = 0
    for i in range(0, len(D0)):
        for j in range(0, len(D0[0])):
            E0el += D0[i,j]*(Hcore[i,j]+Hcore[i,j])
    
    #SCF iterations
    
    if print_SCF == 'Yes':
        print("Iter\t Eel\t")
        """
        output = open('out.txt', 'a')
        output.write('Iter')
        output.write("\t")
        output.write('Eel')
        output.write("\t \t \t \t \t")
        output.write('Etot')
        output.write("\t \t \t \t")
        output.write('dE')
        output.write("\t \t \t \t \t")
        output.write('rmsD')
        if DO_DIIS == 'Yes':
            output.write("\t \t \t \t \t")
            output.write('DIIS')
        output.write("\n")
        output.write('0')
        output.write("\t \t")
        output.write("{:14.10f}".format(E0el))
        output.write("\t \t")
        output.write("{:14.10f}".format(E0el+VNN[0]))
        """
    
    for iter in range(1, Maxiter):
        #New Fock Matrix
        J = np.einsum('pqrs,sr->pq', Vee,D0)
        K = np.einsum('psqr,sr->pq', Vee,D0)
        F = Hcore + 2.0*J-K
        
        if DO_DIIS == 'Yes':
            #Estimate F by DIIS
            if iter == 1:
                F, errorFock, errorDens, errorDIIS = DIIS.runDIIS(F,D0,S,iter,DIIS_steps,0,0)
            else:
                F, errorFock, errorDens, errorDIIS  = DIIS.runDIIS(F,D0,S,iter,DIIS_steps,errorFock, errorDens)
            
        Fprime = np.dot(np.dot(np.transpose(S_sqrt),F),S_sqrt)
        eps, Cprime = diagonlize(Fprime)
        
        C = np.dot(S_sqrt, Cprime)
        
        CT = np.matrix.transpose(C)
        CTocc = CT[0:int(input[0,0]/2)]
        Cocc = np.matrix.transpose(CTocc)
        
        D = np.dot(Cocc, CTocc)
        
        #New SCF Energy
        Eel = 0
        for i in range(0, len(D)):
            for j in range(0, len(D[0])):
                Eel += D[i,j]*(Hcore[i,j]+F[i,j])

        #Convergance
        dE = Eel - E0el
        rmsD = 0
        for i in range(0, len(D0)):
            for j in range(0, len(D0[0])):
                rmsD += (D[i,j] - D0[i,j])**2
        rmsD = (rmsD)**0.5
        
        if print_SCF == 'Yes':
            if DO_DIIS == 'Yes':
                if errorDIIS != 'None':
                    print(iter,"\t","{:14.10f}".format(Eel),"\t","{:14.10f}".format(Eel+VNN[0]),"\t","{: 12.8e}".format(dE),"\t","{: 12.8e}".format(rmsD),"\t","{: 12.8e}".format(errorDIIS))
                else:
                    print(iter,"\t","{:14.10f}".format(Eel),"\t","{:14.10f}".format(Eel+VNN[0]),"\t","{: 12.8e}".format(dE),"\t","{: 12.8e}".format(rmsD))
            else:
                print(iter,"\t","{:14.10f}".format(Eel),"\t","{:14.10f}".format(Eel+VNN[0]),"\t","{: 12.8e}".format(dE),"\t","{: 12.8e}".format(rmsD))
    
        D0 = D
        E0el = Eel
        if np.abs(dE) < deTHR and rmsD < rmsTHR:
            break

    
    return Eel+VNN[0], C, F, 2*D

It depends on a DIIS.py subroutine which is only a few dozen lines long:

import numpy as np
from time import sleep

def DIIS(F,D,S,Efock,Edens,numbF):
    
    Edens[1:] = Edens[:len(Edens)-1]
    Efock[1:] = Efock[:len(Efock)-1]
    Edens[0] = D
    Efock[0] = F
    
    Err = np.zeros((numbF,len(F),len(F)))
    for i in range(0, numbF):
        Err[i] =  np.dot(np.dot(Efock[i],Edens[i]),S)
        Err[i] -= np.dot(np.dot(S, Edens[i]),Efock[i])
    Emax = np.max(np.abs(Err))
    
    #Construct B and b0
    B = -1*np.ones((numbF+1,numbF+1))
    B[numbF,numbF] = 0
    b0 = np.zeros(numbF+1)
    b0[numbF] = -1
    
    for i in range(0, numbF):
        for j in range(0, numbF):
            B[i,j] = np.trace(np.dot(Err[i], Err[j]))
    
    c = np.linalg.solve(B, b0)

    Fprime = np.zeros((len(F),len(F)))
    
    for i in range(0, numbF):
        Fprime += c[i]*Efock[i]
    
    return Fprime, Efock, Edens, Emax
    

def runDIIS(Fnew,Dnew,Snew,iter,Steps,errF,errD):
    
    if iter == 1:
        Efock = np.zeros((Steps, len(Fnew), len(Fnew)))
        Edens = np.zeros((Steps, len(Fnew), len(Fnew)))
        Efock[0] = Fnew
        Edens[0] = Dnew
        F = Fnew
        Emax = 'None'
    else:
        if iter < Steps:
            F, Efock, Edens, Emax = DIIS(Fnew,Dnew,Snew,errF,errD,iter)
        else:
            F, Efock, Edens, Emax = DIIS(Fnew,Dnew,Snew,errF,errD,Steps)
    
    return F, Efock, Edens, Emax

While it appears above that the Hartree-Fock algorithm was implemented in SlowQuant with only about 100 lines of actual code, there's also about 150 more lines of code for the integrals. This could all be made more compact at the expense of readability, and certainly MATLAB (Octave) and/or Julia would allow for even more compact code.

Answer 2

As mentioned above, there is practically a myriad of free working quantum chemistry codes with implementations of Hartree-Fock and Density-Functional Theory.

As a counter-example to the code mentioned above, I want to highlight a recent quantum chemistry code:

The QUantum Interaction Computational Kernel (QUICK)

• It's FOSS (Free and Open Source Software), meaning you can become a developer!
• It's easy to use, and made for speed. It's a good code for large-scale and multi-scale quantum chemistry, with a niche but growing community.
• The code is young but is made "the right way", meaning that at it's core is the acceleration of electron repulsion integrals (ERIs).
• It has MPI parallelization for CPU platforms
• Massive parallel GPU acceleration through CUDA for NVIDIA GPUs
• Multi-GPU support via MPI + CUDA, also across multiple compute nodes
• HF, UHF, DFT and QM/MM (with some new developments soon... stay tuned!)
• Makes use of LIBXC, so you can use your favorite functional!
• Supports DFT-D (e.g. can perform B3LYP-D3(BJ) calculations, wB97X-D)

Here is a sample DFT input file for pSB3 from the QUICK test suite:

 DFT LIBXC=HYB_GGA_XC_B3LYP BASIS=6-31G DIPOLE CUTOFF=1.0d-10 DENSERMS=1.0d-6 CHARGE=+1 
 C -2.74724163 -0.83655480  0.85891890
 C -1.45690243 -0.47166414  0.99917288
 C -0.62772841 -0.22145348 -0.15324144
 C  0.68944541  0.15260156 -0.20171919
 C  1.48823343  0.36448923  0.95078019
 N  2.73140279  0.71794292  0.90370531
 H  1.15299007  0.29662735 -1.16123028
 H -1.11028708 -0.34629454 -1.10667741
 H  1.08266370  0.23672479  1.93583665
 H -1.04825008 -0.36804581  1.98774361
 H  3.26866760  0.85959058  1.73675486
 H  3.20838915  0.86445595  0.03379823
 H -3.36885475 -1.02411938  1.71303219
 H -3.20113376 -0.95331433 -0.10812450
   #TOTAL_ENERGY=  -249.766901607

• This input file calls the LIBXC version of B3LYP; there is a native B3LYP functional as well in QUICK.

Easy to git ; easy to compile; easy to run. If you're interested in looking "under the hood", the modules (e.g. one-electron integrals module) are well-documented and referenced. The DIIS algorithm is the default, and you can choose from several basis sets. It's pretty solid!

You can get QUICK here: QUICK GitHub

Learn QUICK here with tutorials: QUICK user manual

1. DFT implementation paper in JCTC here
2. QM/MM implementation paper in JCIM here
3. Massive multi-GPU acceleration of HF and DFT paper in JCTC here

I think QUICK is a nice code: easy to use, free and fast! Definitely check it out.

Answer 3

I will complement the above answers by highlighting the work that culminated in Psi4NumPy distributed on github and described in a paper. Here the authors focused on creating a workflow to quickly and efficiently develop new algorithms and implementations for useful quantum chemical code.

The main psi4 is distributed separately, but the freedom to use all of python's libraries and numpy's array programming with a fast and FOSS package is notable.

In my experience, the relative performance of the few free (as in beer) computational chemistry software packages I have actually used for production purposes is as follows:

• Scripting: psi4 > ORCA > NWChem
• Stability: ORCA > NWChem > psi4
• Speed: ORCA > psi4 > NWChem
• Accuracy: ORCA > psi4 > NWChem
• Freedom: psi4 > NWChem > ORCA