Programming Projects: Unrestricted MP2
Overview
In this project, you will develop code to evaluate the electronic energy using unrestricted second-order Møller-Pleset perturbation theory (MP2) within the Psi4 electronic structure package. MP2 is the simplest wave-function-based post-Hartree-Fock method for evaluating the dynamical correlation energy.
The MP2 correlation energy
In spin-orbital notation, the MP2 correlation energy is given by
$$E_\text{MP2} = E_\text{HF} + \frac{1}{4} \sum_{ijab} \frac{|\langle ij || ab \rangle|^2}{\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b}$$ Here, EHF represents the unrestricted Hartree-Fock reference energy, which includes the nuclear repulsion energy. The four-index sum reflects the idea that the MP2 energy includes contributions to the correlation energy that arise from all doubly substituted electronic configurations; the factor of 1/4 ensures that we do not double count excitations out of the reference. The indices i and j refer to spin orbitals that are occupied in the reference function, and the labels a and b refer to spin orbitals that are unoccupied in the reference function (virtual orbitals). Again, these sumations are over spin orbitals, which means that the index i spans all occupied orbitals of both α and β spin.
The denominator in the MP2 energy expression consists of molecular orbital energies; these quantities are the eigenvalues of the α- and β-spin Fock matrices. The symbol <ij||ab> denotes an antisymmetrized two-electron repulsion integral in physicists' notation:
$$\langle ij || ab \rangle = \langle ij | ab \rangle - \langle ij | ba \rangle $$ The two-electron integrals in physicists' notation are related to those in chemists' notation used in the previous tutorials by
$$\langle ij | ab \rangle = (ia|jb) $$ Note that the orbital labels refer to molecular spin orbitals. The molecular orbital (MO) basis is the basis of orthonormal orbitals generated by the Hartree-Fock procedure that diagonalizes the α- and β-spin Fock matrices. So, given a set of two-electron integrals in the atomic orbital basis, (μν|λσ), the MO-basis integrals can be obtained through the orbital transformation
$$(ia|jb) = \sum_{\mu\nu\lambda\sigma} (\mu\nu|\lambda \sigma) C_{\mu i} C_{\nu a} C_{\lambda j} C_{\sigma b}$$ Here, the transformation matrix, C, is the same orbital transformation matrix determined through the Hartree-Fock procedure. Note, however, that this transformation is given here in a very naive form. As written, this transformation would scale with the eigth power of system size! Instead, one should perform a series of partial transformations that each scale only as the fifth power of system size:
$$(i\nu|\lambda\sigma) = \sum_{\mu} (\mu\nu|\lambda \sigma) C_{\mu i} $$ $$(ia|\lambda\sigma) = \sum_{\nu} (i\nu|\lambda \sigma) C_{\nu a} $$ $$(ia|j\sigma) = \sum_{\lambda} (ia|\lambda \sigma) C_{\lambda j} $$ $$(ia|jb) = \sum_{\sigma} (ia|j \sigma) C_{\sigma b} $$ This integral transformation step is the source of the overall fifth-power scaling of the MP2 approach.
Density fitting in MP2
As in Hartree-Fock theory, the density fitting (DF) approximation to the two-electron repulsion integrals can reduce the storage requirements and computational effort required to evaluate the MP2 energy. The primary utility of the DF approximation is in reducing the scaling of the integral transformation step for MP2 from the fifth power with system size to the fourth power with system size. The three-index integral transformation can be carried out as two partial transformations that carry fourth-power scaling:
$$(i\nu|Q) = \sum_\mu (\mu\nu|Q) C_{\mu i}$$ $$(ia|Q) = \sum_\nu (i\nu|Q) C_{\nu a}$$ The MO-basis four-index electron repulsion integral (ERI) tensor can then be constructed at fifth-power cost as
$$(ia|jb) = \sum_Q (ia|Q)(Q|jb)$$ So, the DF approximation does not reduce the formal scaling of the MP2 energy evaluation, but it does reduce the number of fifth-power-scaling steps. Conventional MP2 requires four fifth-power-scaling partial integral transformations, while DF-MP2 requires only one fifth-power-scaling contraction to construct the four-index ERI tensor.
It is important to note here that the auxiliary basis used within a DF-MP2 procedure is usually different than that used within a density-fitted Hartree-Fock procedure. The reason is that the auxiliary basis sets are separately optimized to give accurate results for either Hartree-Fock or MP2, but not both. Keep this in mind when comparing your results to those from Psi4's built-in DF-MP2 function.
Procedure
If you have completed the previous Hartree-Fock tutorials, you can add the MP2 energy evaluation to the end of your plugin, after the Hartree-Fock procedure converges. Otherwise, create a new Psi4 plugin, as described here. Before moving on, open the input file in your plugin directory (input.dat) and look at the molecule block
molecule {
O
H 1 R
H 1 R 2 A
R = .9
A = 104.5
# add this line!
symmetry c1
}
This plugin will not make use of point group symmetry, so you should add a line specifying that the computation be performed using C1 symmetry. Now, add a flag indicating that the computation will use density fitting techniques in the Hartree-Fock procedure:
set {
#add this line!
scf_type df
basis sto-3g
}
For the remainder of this project, I've assumed that the rest of the input is unchanged. Keep this in mind when comparing to my results.
Now, we must make sure that Psi4 also recognizes that we will be using the DF approximation within our MP2 procedure. Close input.dat and open pymodule.py. Locate the line that contains "Ensure IWL has been written". The line following this one is meant to ensure that a plugin has access to the four-index integrals in cases where the Hartree-Fock procedure used density fitting. Since we will use the DF approximation in our MP2 procedure, we do not need four-index integrals. Comment out that line and add a few additional lines that ensure that the reference wave function passed into the plugin knows what density fitting basis set is appropriate for our MP2 procedure:
# Ensure IWL files have been written when not using DF/CD
#proc_util.check_iwl_file_from_scf_type(psi4.core.get_option('SCF', 'SCF_TYPE'), ref_wfn)
aux_basis = psi4.core.BasisSet.build(ref_wfn.molecule(), "DF_BASIS_MP2",
psi4.core.get_option("DFMP2", "DF_BASIS_MP2"),
"RIFIT", psi4.core.get_global_option('BASIS'))
ref_wfn.set_basisset("DF_BASIS_MP2", aux_basis)
In order to use Psi4's density fitting capabilities and the Matrix, Vector, and Wavefunction classes, we must include the corresponding headers in plugin.cc. If you are making this plugin from scratch, then you will need to add the following:
#include "psi4/libmints/wavefunction.h" #include "psi4/libmints/matrix.h" #include "psi4/libmints/vector.h" #include "psi4/libmints/basisset.h" #include "psi4/lib3index/dftensor.h" #include "psi4/libqt/qt.h"
As in the previous tutorials, we can ask the reference wave function passed into the plugin for some information about the system (the number of α electrons, the orbital transformation matrices, etc.):
// total number of alpha electrons
int na = ref_wfn->nalpha();
// total number of beta electrons
int nb = ref_wfn->nbeta();
// total number of basis functions
int nso = ref_wfn->nso();
// total number of molecular orbitals
int nmo = ref_wfn->nmo();
std::shared_ptr<Matrix> Ca = ref_wfn->Ca();
std::shared_ptr<Matrix> Cb = ref_wfn->Cb();
If you are writing your MP2 code as an extension to your previous Hartree-Fock plugin, use the orbital transformation matrices optimized by your procedure. Now, use the DFTensor class to generate three-index integrals in the atomic orbital basis. Be sure to use the auxiliary basis set that is appropriate for MP2, as opposed to that for Hartree-Fock theory (Psi4 refers to this basis as "DF_BASIS_MP2").
// grab the molecule from the wavefunction that was passed into the plugin
std::shared_ptr<Molecule> mol = ref_wfn->molecule();
// get primary basis:
std::shared_ptr<BasisSet> primary = ref_wfn->get_basisset("ORBITAL");
// get auxiliary basis:
std::shared_ptr<BasisSet> auxiliary = ref_wfn->get_basisset("DF_BASIS_MP2");
// total number of auxiliary basis functions
int nQ = auxiliary->nbf();
// construct the three-index integrals
std::shared_ptr<DFTensor> DF (new DFTensor(primary,auxiliary,Ca,na,nso-na,nb,nso-nb,options));
std::shared_ptr<Matrix> Qso = DF->Qso();
You can now use the α- and β-spin MO transformatrices to transform the three-index integrals from the atomic orbital basis to the molecular orbital basis. For unrestricted MP2, you will need to separately transform α- and β-spin tensors so that you can eventually construct (αα|αα), (αα|ββ), and (ββ|ββ) type integrals.
The last necessary ingredient for the MP2 energy is the list of molecular orbital energies. You can grab these from the reference wavefunction:
std::shared_ptr<Vector> epsilon_a = std::shared_ptr<Vector>(new Vector(nmo));
epsilon_a->copy(reference_wavefunction_->epsilon_a().get());
std::shared_ptr<Vector> epsilon_b = std::shared_ptr<Vector>(new Vector(nmo));
epsilon_b->copy(reference_wavefunction_->epsilon_b().get());
If you are writing this plugin as an extension of a previously completed Hartree-Fock plugin, then you need not ask the reference wave function for the orbital energies. Instead, you should use the eigenvalues of the Fock matrix constructed using your converged Hartree-Fock solution.
For the water molecule in a minimal basis set (using the default input file with the small modifications given above), the DF-MP2 correlation energy should be -0.031081575913 Eh. Once your code can correctly recover the correlation energy for this neutral species, verify that you recover the correct correlation energy for an open-shell molecule:
molecule {
1 2
O
H 1 R
H 1 R 2 A
R = .9
A = 104.5
# add this line!
symmetry c1
}
set {
basis sto-3g
# add these lines!
scf_type df
reference uhf
}
energy('mymp2')
In this case, the DF-MP2 correlation energy should be -0.024767575359 Eh.