# Calclulation of MP2 first order wave function

My question is, how can I calculate the MP2 first order wavefunction, and thus construct the MP2 density matrix? Now here is what I understand so far:

The density matrix is given as:

$$P_{ij}=\sum_{a}C_{ia}(C_{ja})^*$$

Here $C$ is the wave function coefficient matrix. To get the density matrix for the MP2 first order wavefunction, their coefficients need to be found.

The first-order wave function is given as [1]:

$$\left|\left.\varPsi_{i}^{(1)}\right> \right.=\sum_{n}c_{n}^{(1)}\left|\left.\varPsi_{n}^{(0)}\right> \right.$$

Now this gives:

$$\left< \left.\varPsi_{n}^{(0)}\right|\varPsi_{i}^{(1)}\right> =c_{n}^{(1)}$$

Which by some more steps, that can be seen in Szabo and Ostlund leads to:

$$c_{n}^{(1)}=\left< \left.\varPsi_{n}^{(0)}\right|\varPsi_{i}^{(1)}\right> =\frac{\left< \varPsi_{n}^{(0)}\left|\mathscr{V}\right|\varPsi_{i}^{(0)}\right> }{E_{i}^{(0)}-E_{n}^{(0)}}$$

In Szabo and Ostlund the energy is derived to have the expression:

$$E_{i}^{(2)}=\sum_{n}'\frac{\left< \varPsi_{i}^{(0)}\left|\mathscr{V}\right|\varPsi_{n}^{(0)}\right> \left< \varPsi_{n}^{(0)}\left|\mathscr{V}\right|\varPsi_{i}^{(0)}\right> }{E_{i}^{(0)}-E_{n}^{(0)}}$$

Though it can also be written as:

$$E_{i}^{(2)}=\sum_{n}'c_{n}^{(1)}\left< \varPsi_{i}^{(0)}\left|\mathscr{V}\right|\varPsi_{n}^{(0)}\right>$$

Now forwarding to the energy expression for the MP2 energy:

$$E_{0}^{(2)}=\sum_{i<j,a<b}\frac{\left< \varPsi_{0}^{(0)}\left|\mathscr{V}\right|\varPsi_{ij}^{ab}\right>\left< \varPsi_{ij}^{ab}\left|\mathscr{V}\right|\varPsi_{0}^{0}\right>}{E_{0}^{(0)}-E_{ij}^{ab}}=\sum_{i<j,a<b}\frac{\left< ij\left|\right|ab\right>\left< ab\left|\right|ij\right>}{\varepsilon_{i}+\varepsilon_{j}-\varepsilon_{a}-\varepsilon_{b}}$$

Now I would like to think that I could identify the wave function coefficients by comparing the two above equations to get:

$$c_{ijab}^{(1)}=\frac{\left< ab\left|\right|ij\right>}{\varepsilon_{i}+\varepsilon_{j}-\varepsilon_{a}-\varepsilon_{b}}$$

I would have expected that I would have found a vector of $c$ such that I would be able to construct the coefficent matrix as:

$$C^{MP2}=c\cdot c^T \Rightarrow P^{MP2}_{ij}=\sum_{a}C^{MP2}_{ia}(C^{MP2}_{ja})$$

But instead I got a four index coefficent out. Where did my line of thinking go wrong? And how can I proceed to construct the first order wave function, and thus get the density matrix?

References to old papers where this have been done explicit is also appreciated, because I have not been able to find any thing but recent papers about MP2 and density matrix fitting, which is not what I am looking for.

[1] Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Dover: Mineola, NY, 1989.

• You have two different kinds of coefficients here. The density matrix in your first equation is calculated from MO coefficients $C$. While the coefficients $c$ are from the development of the MP2 wave function in terms of eigenfunctions of the unpertubated Hamiltonian. Thus, those are two different density matrices you are talking about. Sep 24 '17 at 8:37
• I see I have made some bad notation. So my goal is to end up with the density matrix from the MP2 wave function. Now I thought that all (not just MO) density matricies could be made by the first equation. Assuming that the equation to find the coefficient $c_{ijab}^{(1)}$ is correct, the thing I don't get it how I go from four indecies to a single index. Sep 24 '17 at 8:44
• One might simply reshape it, but I don't know if this is helpful here. But why do you need a single index? Density matrices may have any number of indices (not sure about odd numbers, but at least even numbers). Sep 24 '17 at 8:51
• I would have thought that the $c^{(1)}$ would be a vector, thus making $C^{MP2}$ a marix and thus making $P^{MP2}$ a matrix. Now I think I would like the denisty matrix to be of two dimension, because then it would fit directly into the equation to find the molecular dipole moment. And it could be diagonalized to find the Natural Orbital occupancy numbers. Sep 24 '17 at 8:56
• Where did you get your Ansatz for $c^{(1)}$ from? I am not sure if you can do this ... Sep 24 '17 at 9:57

The main problem with the last statement,

\begin{align} C^{MP2} &= c\cdot c^T \\ \Rightarrow P^{MP2}_{ij} &= \sum_{a}C^{MP2}_{ia}(C^{MP2}_{ja}), \end{align}

is that one never has "MP2 MO coefficients" in this way; you would have MP2 natural orbital (NO) coefficients that come from diagonalizing this density, so it's the other way around:

$$P^{\text{MP2}}C^{\text{MP2 NO}} = n^{\text{MP2 NO}}C^{\text{MP2 NO}},$$

where $n^{\text{MP2 NO}}$ is the diagonal matrix of natural orbital occupation numbers. This still leaves the question of how to calculate $P^{\text{MP2}}$. Conceptually, the way I view it is that the set of $\{c_{ijab}^{(1)}\}$ are the amplitudes describing how to correct the uncorrelated ground-state wavefunction; they don't specify the wavefunction on their own, so $C^{\text{SCF}}$ will probably have to appear in the total density.

All of the following assumed a restricted wavefunction. $i,j,k,\dots$ are occupied MO indices, $a,b,c,\dots$ are virtual/unoccupied MO indices, $p,q,r,\dots$ are general (any) MO indices, and $\mu,\nu,\dots$ are general AO indices. The 2nd-order correction to the density, specifically the occupied-occupied and virtual-virtual blocks are given by Eqn. (11a) in Trucks [1]:

\begin{align} P_{ij}^{\text{MP2}} &= -\frac{1}{2} \sum_{kab} c_{jk}^{ab} c_{ik}^{ab} \\ P_{ab}^{\text{MP2}} &= +\frac{1}{2} \sum_{ijc} c_{ij}^{ac} c_{ij}^{bc}, \end{align}

where I switched to using your notation. This is indeed similar to your idea from above; accumulate over both indices of the opposite space, then do a matrix multiplication.

The total density is then [Eqn. (7)]

$$\mathscr{P} = P^{\text{SCF}} + P^{\text{MP2}},$$

where $P_{pq}^{\text{SCF}} = \delta_{pq}$ because so far we're working in the MO basis. The back-transformation to the AO basis will be the same as for any other 2-index MO-basis quantity:

\begin{align} P_{\mu\nu} &= \sum_{pq}^{\text{all MOs}} C_{\mu p} P_{pq} (C_{\nu q})^{T} \\ \Rightarrow \mathbf{P}^{\text{AO}} &= \mathbf{C} \cdot \mathbf{P}^{\text{MO}} \cdot \mathbf{C}^{T}. \end{align}

Then, we can evaluate molecular properties that are represented as energy derivatives [Eqn. (6)]:

$$\left.\frac{\partial \mathscr{E}}{\partial \lambda}\right|_{\lambda=0} = \sum_{\mu\nu} \mathscr{P}_{\mu\nu} Q_{\mu\nu}^{\lambda}$$

where $\mathscr{E} = E_0^{(0)} + E_0^{(2)}$ and $Q_{\mu\nu}^{\lambda}$ is the matrix representation of the external perturbation $\lambda$ in the AO basis. This is analogous to the SCF result:

$$\left.\frac{\partial E^{\text{SCF}}}{\partial \lambda}\right|_{\lambda=0} = \sum_{\mu\nu} P_{\mu\nu}^{\text{SCF}} Q_{\mu\nu}^{\lambda}$$

## Unrelaxed vs. relaxed density

The above definition of the density is incomplete, since it describes the unrelaxed correlated density. Our usual formulation of correlated methods is to describe correlation as some to-be-determined mixture of the virtual orbitals with the occupied orbitals. There is no such mixing in the above equations. In reality, $P_{ia}^{\text{MP2}} \overset{!}{\neq} 0$, so there is something that looks like a "transition density" with non-zero occupied-virtual blocks, in order to correct for the fact that our $C^{\text{SCF}}$ has no correlation. Allowing this mixing (which corresponds to orbital response) leads to the relaxed density. However, one has to solve the coupled-perturbed Hartree-Fock (CPHF or CPSCF) equations.

See between Eqns. (8) and (9) in [1] for the CPHF equations and how to solve them for the occupied-virtual mixing coefficients and finally the $P_{ia}^{\text{MP2}}$ block. Generally, note the similarities to the CIS density equations.

## CPHF equations and first-order properties

As stated, when differentiating the energy with respect to some perturbation $\lambda$ to calculate first-order properties, one needs to contract the density with the operator property matrix. If you transform Eqn. (6) from the AO into the MO basis,

\begin{align} \left.\frac{\partial \mathscr{E}}{\partial \lambda}\right|_{\lambda=0} &= \sum_{pq} \mathscr{P}_{pq} Q_{pq}^{\lambda} \\ &= \sum_{pq} P_{pq}^{\text{HF}} Q_{pq}^{\lambda} + \sum_{pq} P_{pq}^{\text{MP2}} Q_{pq}^{\lambda}, \end{align}

the SCF contribution simplifies due to the nature of the SCF density matrix in the MO basis,

$$P_{pq}^{\text{HF}} = \left\{ \begin{array}{lr} 2 & \text{if }p = q\text{, p is an occupied MO} \\ 0 & \text{otherwise} \end{array} \right.,$$

\begin{align} \left<\hat{\lambda}\right> &= \sum_{pq} P_{pq}^{\text{HF}} Q_{pq}^{\lambda} \\ &= 2 \sum_{i} Q_{ii}^{\lambda} \\ &= \mathrm{tr}\left(\mathbf{P}^{\text{AO}}\mathbf{Q}^{\text{AO}}\right) \\ &= \mathrm{tr}\left(\mathbf{P}^{\text{MO}}\mathbf{Q}^{\text{MO}}\right) \end{align}

for one-electron operators and wavefunctions that obey the Hellmann-Feynman theorem. I absorbed the 2 into $\mathbf{Q}^{\text{MO}}$, but this depends on the program implementations and their definition of the density. The point is that this contribution is relatively simple, and illustrates how for one-electron operators (dipole, quadrupole, angular momentum, ...), evaluating them for a converged Hellmann-Feynman wavefunction is very fast: one the density and property integrals are in the same basis, either do a matrix multiply and then trace, or do an element-wise (Schur) product and accumulate all the elements, so it scales as $O(N^3)$. No solution of CPHF equations is required. The same goes for the gradient with respect to nuclear coordinates $\frac{\partial E^{\text{SCF}}}{\partial x}$; it isn't a one-electron operator, and requires a mathematical trick to make $\frac{\partial P^{\text{SCF}}}{\partial x}$ disappear, but the SCF gradient doesn't require solving the CPHF equations either. This isn't true for $\frac{\partial E^{\text{MP2}}}{\partial x}$, which requires $\frac{\partial P^{\text{MP2}}}{\partial x}$, where you need to solve the CPHF equations to even form the full $P^{\text{MP2}}$.

Consider the remaining MP2 part of the property evaluation; the total density contains contributions from the occupied-occupied (oo) block, the virtual-virtual (vv) block, and the occupied-virtual (ov) block. I'm going to switch from $P^{\text{MP2}}$ and $c$ to $D$ and $t$ to make typing easier.

$$\sum_{pq} D_{pq} Q_{pq}^{\lambda} = \sum_{ij} D_{ij} Q_{ij}^{\lambda} + \sum_{ab} D_{ab} Q_{ab}^{\lambda} + 2\sum_{ia} D_{ia} Q_{ia}^{\lambda}$$

The first two terms can already be calculated from the unrelaxed density; $D_{ia}$ is the term that must come from solving the CPHF equations.

In general, the form of solving CPHF equations looks like $\mathbf{AU} = -\mathbf{V}$, where $\mathbf{A}$ is the orbital Hessian (2nd derivative of the energy with respect to orbital rotations), $\mathbf{V}$ is the gradient of some perturbation $\hat{V}$. In [1], it's $\lambda$, but Neese [3,4] and the Scandinavians usually use $V$ or $Q$ when the perturbation is associated with an observable. This perturbation gradient, termed the right-hand side (or RHS), describes how the perturbation might cause occupation of virtual orbitals, and the orbital Hessian describes how the response of one MO to some perturbation is coupled to the other MOs [5,6].

Start by rewriting $\sum_{ia} D_{ia} Q_{ia}^{\lambda}$ as the combination of a perturbation-dependent response $\mathbf{U}$ and the perturbation gradient $\mathbf{X}$:

$$2\sum_{ia} D_{ia} Q_{ia}^{\lambda} = 2\sum_{ia} X_{ai} U_{ai}^{\lambda}$$

$\mathbf{X}$ serves the same purpose as the $\mathbf{V}$ mentioned above, but rather than coming from an external perturbation, it is an intermediate quantity that describes the MP2 wavefunction.

First, the perturbation-dependent response. This part is not specific to MP2, and is also important for higher-order SCF-based molecular response [3].

\begin{align} U_{ai}^{\lambda} &= \sum_{bj} (\mathbf{A}^{-1})_{ai,bj} Q_{bj}^{\lambda} (\varepsilon_{j} - \varepsilon_{b})^{-1} \\ A_{ai,bj} &= 1 + \left( \left<ab||ij\right> - \left<aj||ib\right> \right) (\varepsilon_{a} - \varepsilon_{b})^{-1} \end{align}

The $\varepsilon$ are MO energies. If we don't assume the use of canonical MOs, this energy denominator is actually $F_{ab}S_{ij} - F_{ij}S_{ab}$.

Jumping ahead, the relationship between the ov density, the orbital Hessian, and the MP2 RHS is [Eqn. (9)]

$$\sum_{bj} D_{bj} (\varepsilon_{j} - \varepsilon_{b}) A_{bj,ai} = X_{ai}$$

See how this is structurally similar to solving for $U$, but the RHS is different. To show the equivalence between all quantities,

\begin{align} 2\sum_{ai} X_{ai} U_{ai}^{\lambda} &= 2\sum_{ai} X_{ai} \sum_{bj} (\mathbf{A}^{-1})_{ai,bj} Q_{bj}^{\lambda} (\varepsilon_{j} - \varepsilon_{b})^{-1} \\ &= 2\sum_{bj} D_{bj} Q_{bj}^{\lambda}. \end{align}

The explicit equation for $X_{ai}$ is [Eqn. (11b)]

$$X_{ai} = \sum_{jk} D_{kj} \left<ki||ja\right> + \sum_{bc} D_{bc} \left<bi||ca\right> - \frac{1}{2} \left[ t_{kj}^{ab} \left<ib||kj\right> + t_{ij}^{cb} \left<cb||aj\right> \right],$$

which is in fact the MP2 Lagrangian, $L_{ai}$ [2, Eqn. (11)]. Again, notice the similarity to the CIS Lagrangian. The Lagrangian is required because the wavefunction doesn't need to be relaxed in order to calculate the correlation energy, so you pay the price here. These equations must be evaluated for MP2 geometry optimizations (during the gradient part) or for fully-relaxed MP2 molecular properties. To be clear, if one is only interested in the density and not molecular properties, only $\mathbf{DA} = \mathbf{X}$ needs to be solved for, and not the equations involving $\mathbf{U}$. In terms of scaling, SCF-based CPHF is about $O(N^4)$ (after back-transforming the $\mathbf{A}^{-1}\mathbf{Q}$ terms to be AO-direct) but might require fewer iterations than the SCF itself. According to [2], the MP2 gradient is something like $O(ON^4) + O(O^{2}V^{3})$, but they can do similar tricks during iterations to cut down on the number of fifth-order steps. In my experience, MP2 gradient calculations are about as expensive as an SCF polarizability calculation.

### References

1. Trucks, G. W.; Salter, E. A.; Sosa, C.; Bartlett, R. J. Theory and implementation of the MBPT density matrix. An application to one-electron properties. Chem. Phys. Lett. 1988, 147 (4), 359-366. DOI: 10.1016/0009-2614(88)80249-5

2. Frisch, M. J.; Head-Gordon, M.; Pople, J. A. A Direct MP2 gradient method. Chem. Phys. Lett. 1990, 166 (3), 275-280. DOI: 10.1016/0009-2614(90)80029-D

3. Neese, F. Prediction of molecular properties and molecular spectroscopy with density functional theory: From fundamental theory to exchange-coupling. Coordination Chemistry Reviews 2009, 253 (5-6), 526-563. DOI: 10.1016/j.ccr.2008.05.014

4. I cannot recommend the above paper enough; it contains excellent descriptions of how electric and magnetic response are calculated with explicit equations that are actually usable. There, you can see how CPHF is used for molecular response properties and the connection with excitation energies, found by solving $\mathbf{AU} = \mathbf{0}$ (diagonalization rather than solving a set of linear equations).

5. Toulouse, J. Introduction to the calculation of molecular properties by response theory.

6. Parrish, R. M. Psi4 Development: Theory Supplement, CPHF for Dipole Polarizabilities.

• If you would like an elaboration on how the CPHF equations are solved and what each term means, let me know and that part can be expanded. Unlike correlated methods, I actually know how CPHF works. Sep 24 '17 at 18:49
• An elaboration of how the CPHF equations are solved would be very appricaited! Sep 24 '17 at 19:37