aromanro/hartreefock

A program implementing the Hartree–Fock (also post-HF: MP2, CCSD(T), CIS and TDHF/RPA)/self-consistent field method (also DIIS) with Gaussian orbitals

HartreeFock

A program implementing the Hartree–Fock/self-consistent field method (and more, see the README at the end)

Description is available here: http://compphys.go.ro/the-hartree-fock-program/

Some Hartree-Fock theory here: http://compphys.go.ro/the-hartree-fock-method/

Some things more general (Schrödinger equation, Born-Oppenheimer approximation, variational principle), here: http://compphys.go.ro/how-to-solve-a-quantum-many-body-problem/

PROGRAM IN ACTION

HOW TO COMPUTE

Here are some things about usage:

Using the classes should be easy. Here is how to grab some atoms from the 'basis':

	Systems::AtomWithShells H1, H2, O, N, C, He, Li, Ne, Ar;

	for (auto &atom : basis.atoms)
	{
		if (atom.Z == 1) H1 = H2 = atom;
		else if (atom.Z == 2) He = atom;
		else if (atom.Z == 3) Li = atom;
		else if (atom.Z == 8) O = atom;
		else if (atom.Z == 6) C = atom;
		else if (atom.Z == 7) N = atom;
		else if (atom.Z == 10) Ne = atom;
		else if (atom.Z == 18) Ar = atom;
	}

Here is how to set the H2O molecule with the coordinates from the 'Mathematica Journal' (referenced in the code):

	H1.position.X = H2.position.X = O.position.X = 0;

	H1.position.Y = 1.43233673;
	H1.position.Z = -0.96104039;
	H2.position.Y = -1.43233673;
	H2.position.Z = -0.96104039;

	O.position.Y = 0;
	O.position.Z = 0.24026010;

	Systems::Molecule H2O;

	H2O.atoms.push_back(H1);
	H2O.atoms.push_back(H2);
	H2O.atoms.push_back(O);

	H2O.Init();

And here is how you calculate:

HartreeFock::RestrictedHartreeFock HartreeFockAlgorithm;
HartreeFockAlgorithm.alpha = 0.5;
HartreeFockAlgorithm.initGuess = 0;

HartreeFockAlgorithm.Init(&H2O);
double result = HartreeFockAlgorithm.Calculate();

You can do computation for a single atom, too, for now by putting it into a dummy molecule with a single atom in it. For example for He:

  Systems::Molecule Heatom;
  Heatom.atoms.push_back(He);
  Heatom.Init();

LATEST ADDITIONS

I added more basis sets, it's not limited to STO-nG. While doing that, I found and fixed a bug in the integrals repository that manifested for orbitals having L > 1. Now it seems to work fine.

Basis sets added:

• Split valence orbitals: 3-21G, 6-21G, 6-31G
• Besides 'split valence', polarization on heavy atoms: 6-31G*
• 'Split valence', polarization on heavy atoms and hydrogen and diffusion functions on heavy atoms: 6-31+G**

You may add more of them, the parsed format is nwchem.

Tutorials from Crawford Group

Most of those are already implemented:

https://github.com/CrawfordGroup/ProgrammingProjects

Projects for MP2, DIIS, CCSD(T), CIS and TDHF/RPA are implemented. DIIS and MP2 are implemented for both the restricted and unrestricted methods.

Projects 1 and 2 are not implemented here, maybe I'll add a python workbook for those in the python repository, 3 was already implemented obviously before finding this. Project 7 is also not implemented, also 9 and 13.

hartreefock's Issues

Add post Hartree-Fock method(s) to improve results

For example:
https://en.wikipedia.org/wiki/M%C3%B8ller%E2%80%93Plesset_perturbation_theory
https://en.wikipedia.org/wiki/Coupled_cluster
https://en.wikipedia.org/wiki/Configuration_interaction

3D charting of orbitals

A surface or volumetric visualization with VTK

Vibrational spectra

Having this #16 implemented, vibrational modes and frequencies could be computed.

It could even go further, computing line intensities.

Results are wrong if orbitals with L>1 are involved

First I noticed the wrong results by testing some * basis sets, but it turned out that the results are wrong even for STOnG if a D orbital is involved. As results come lower than the Hartree Fock limit, this is very probably a bug in integrals computation.

Compute gradient, Hessian matrix

The derivatives could be computed numerically (central difference would be my choice).
In fact, they could be computed 'analytically', but I don't think I would use that method for this project, I suppose the post-HF things would complicate things even more.

They could be used for optimizations computations, that is, computation of molecular structure.
Also vibrational modes & frequencies could be computed.
Ideally, also the lines intensity...

Add more basis sets

In principle the program should be able to use a different number of Gaussians in the inner shell than in the valence shell, so it could work with orbital sets like 3-21G or 6-21G. Maybe the parsing code should be changed for that, but it should be possible to use them.

Add DIIS

https://en.wikipedia.org/wiki/DIIS

Implement Davidson-Liu algorithms for CIS

As described in this tutorial: https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2313

Add Car-Parrinello?

Similar, but more complex, with what it's already implemented in python: https://github.com/aromanro/PythonCompphys/blob/master/Car-Parrinello.ipynb

The complexity stems from those derivatives, that's what stopped me to implement it first here. It's way simpler when using only 's orbitals'.

Implement CISD or CID

For now CIS for restricted method is implemented, to be used in TDHF/RPA.

Maybe I'll also implement at least CID.

Add electric field term to the Hamiltonian

Only the 'core' part needs to be changed, it's not so hard.

The advantage would be to be able to calculate the electric dipole moment using the derivative of energy with respect to the electric field.

For details, see 6.101 and the explanation before it from Szabo & Ostlund, Modern Quantum Chemistry. Introduction to Advanced Electronic Structure Theory.

Add DFT

Another thing that probably won't happen soon - but it's something that is possible - is to add DFT to the project.

Add dynamic dipole polarizability and photoabsorption spectrum computation

Together with some other things (not only time dependent but also time independent), they are for example briefly described in
'Introduction to the calculation of molecular properties by response theory' by Joulien Toulouse.

Add * orbitals

Add unoccupied/polarization orbitals to the basis set (as in STO-3G* or 3-21G*).

Add option to chart Hartrees instead of eV

The energy is usually given in Hartrees in the literature, so this would help checking the results

Structural computations

Having this #16 implemented, it's just an optimization problem.

The code could compute the molecule structure.

Implement CC for the unrestricted method

Currently I'm implementing https://en.wikipedia.org/wiki/Coupled_cluster for the restricted method, hopefully that will succeed, but probably I'll postpone indefinitely the implementation for the unrestricted method.

I'll leave the issue opened for the future...

Add population analysis

As in the last paragraph here: https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2303

Add possibility to specify occupation of orbitals

This would allow for example to compute excited states.

Implement CIS and TDHF/RPA for unrestricted method?

I'll see if I'll do it...

Add computation of dipole-moment integrals

Would be nice to be able to compute more, like polarizability or transition dipole moment.
Probably won't happen soon.

Add GW?

Another one possible: https://en.wikipedia.org/wiki/GW_approximation
There are plenty of articles on this (for example: https://aip.scitation.org/doi/10.1063/1.4718428) and a good book, "Interacting Electrons: Theory and Computational Approaches" by Richard Martin et al.