A Fortran Electronic Structure Programme (AFESP)

This project is based on the Crawford Group's excellent C++ Programming Tutorial in Chemistry.

The end goal of this project will be performing HF, MP2, CCSD, and CCSD(T), as per the original tutorial, but with additional support for multicore processors (modern CPUs, GPUs).

A report on this project as submitted to the MPhil in Scientific Computing course at the University of Cambridge, can be found here.

Progress / to-do

• Geometry read-in
• Integral read-in
• Nuclear repulsion energy calculation
• Hartree-Fock without symmetry
• Pretty printing, time reporting
• MP2 without symmetry
• CCSD without symmetry
• OpenMP parallelisation
• DIIS acceleration for Hartree-Fock SCF
• CCSD(T) without symmetry
• DIIS acceleration for CCSD iterations
• Loop optimisations - cache
• User input file: level of theory, tolerances
• BLAS acceleration of tensor contractions
• Spin-free CCSD
• Spin-free CCSD(T)
• Spin-free CCSD[T]
• Renormalised CCSD[T]
• Renormalised CCSD(T)
• Completely renormalised CCSD[T]
• Completely renormalised CCSD(T)
• GPU acceleration for CCSD(T)
• RHF tracking (write-out and read-in of SCF guesses)
• Integral generation
• All-purpose wrapper script
• Adapting Hartree-Fock with symmetry
• MPI parallelisation

Installation guide

To get started, clone this directory by

git clone [email protected]:brianz98/A-Fortran-Electronic-Structure-Programme.git

Dependencies

You need to have:

• CMake (at least 3.12)
• A Fortran compiler (gfortran 9.3.0 have been tested)
• OpenBLAS for linear algebra routines (see below for instructions)
• Fypp for pre-processing (available from pip or conda, see here)

Optional, but only if you want to source for your own integral files:

• Python 3
• Psi4 (available on conda, see here)

Installing OpenBLAS

You may already have OpenBLAS on your system. Type whereis openblas can see if the library is already installed.

If not, you can download the latest OpenBLAS version from their website, and untar, then

make USE_OPENMP=1
make install [PREFIX=/your/directory]

where the PREFIX can be omitted and OpenBLAS will be installed in /opt/OpenBLAS by default, but you'll likely need to prefix sudo to the command.

If the linker can't find the OpenBLAS library, then you need to change the last line of CMakeLists.txt from -lopenblas to where libopenblas.a is, which should be under /your/directory/lib. You may also need to add -lpthread to the last line.

Compilation

The compilation can be handled by the Bash script make.sh. It takes three arguments:

• -n: build from scratch, removes previous build/ directory and els.x binary.
• -p: production-level build, this means more optimisations and less checks. Without this flag a debug build will be carried out, where array boundary checks are enabled, together with gdb debugging options enabled.
• -g: GPU acceleration (currently only for the spin-orbital CCSD(T) module). Requires a fairly recent PGI/NVHPC Fortran compiler. This unfortunately disables OpenMP due to a compatibility issue.

Running calculations

File requirements

A directory with the following files are needed:

• els.in - input settings file, see below
• eri.dat - two-electron (electron repulsion) integrals
• geom.dat - geometry
• s.dat - overlap integrals
• t.dat - kinetic integrals
• v.dat - potential integrals

The integrals can be generated by the utility given in utils/psi4\_integrals\_nosym.py.

Input file

The input file needs to be named els.in and is a Fortran namelist file. It has the following form:

&elsinput
calc_type="CCSD(T)_spinorb",    ! See below for allowed calculation types
scf_e_tol=1e-6,                 ! HF-SCF energy convergence threshold
scf_d_tol=1e-7,                 ! HF-SCF RMS density convergence threshold
scf_diis_n_errmat=6,            ! Number of SCF DIIS error matrices
ccsd_e_tol=1e-6,                ! CCSD energy convergence threshold
ccsd_t_tol=1e-7,                ! CCSD RMS t2 convergence threshold
ccsd_diis_n_errmat=8,           ! Number of CCSD DIIS error matrices
scf_maxiter = 50,               ! Maximum number of SCF iterations
ccsd_maxiter = 50,              ! Maximum number of CCSD iterations
write_fcidump = .false.         ! (Requires at least an MP2 calculation) output ERI in MO basis
scf_read_guess = .false.,       ! Read in the AO Fock matrix from a previous geometry as initial guess 
scf_write_guess = .true.        ! Write out the AO Fock matrix
/

Allowed calculation types are:

• RHF
• MP2_spatial
• CCSD_spinorb
• CCSD_spatial
• CCSD(T)_spinorb
• CCSD(T)_spatial
• CCSD[T]_spatial
• RCCSD(T)_spatial
• RCCSD[T]_spatial
• CRCCSD(T)_spatial
• CRCCSD[T]_spatial

The final energy breakdown will contain all energies up to the highest level of calculation performed.

Using the wrapper script

A fairly capable wrapper script is provided in utils/els_wrapper. The settings are in main and are all self-explanatory. The most common use for it would be to calculate a binding curve. The outputs will be a directory named <molecule-name>-<basis-name>, containing subdirectories named <bondlength>-<angle>, and other files detailed below.

- h2o-cc-pvdz/
    - 1.20_104.45/
        - s.dat ---------------- Standard integral/geometry files
        - t.dat ---------------- Ditto
        - v.dat ---------------- Ditto
        - eri.dat -------------- Ditto
        - geom.dat ------------- Ditto
        - els.in --------------- Input file
        - els.out -------------- els.x standard output
        - guess_in.dat --------- AO Fock from previous geometry
        - guess_out.dat -------- AO Fock from this geometry
        - els_energy.dat ------- Summarised energies from this calculation
        - reference.dat -------- Summarised energies from the Psi4 calculation
    - 1.22_104.45/
        - ...
    - binding_energy_els.dat --- Numpy-readable datafile containing energies from all points
    - binding_energy_psi4.dat -- Same but from Psi4
    - h20-cc-pvdz.psi4out ------ Psi4 standard output

If only a single point is needed, you can also use utils/psi4_integrals_nosym.py, which is essentially the same code as above but without running els.x automatically.