MP2 density matrix
Why
This issue is motivated by the need for efficient communication between Octopus and Qiskit regarding the ongoing collaboration between MPSD and IBM. Indeed, currently, the only way to compute one-body and two-body matrix elements is starting from a Hartree-Fock or a DFT solution and then calculating the matrix elements using those states. However, using those states is inefficient.
As an example, let's consider the Hydrogen molecule. In order to retrieve 75% of the total correlation energy (compared with other codes), it is necessary to use 200 unoccupied states. This leads to a file dimension of roughly 4GB (for the two-body matrix elements). This makes the transfer and the reading/writing of such data very slow.
One solution to be able to retrieve significant portions of the correlation energy with a reduced number of states is to use natural orbitals, which essentially are a linear combination of HF states [1].
How
First, we start with a Hartree-Fock calculation. Octopus supports the computation of the one-body and two-body matrix elements with the theory levels DFT and HF. Let's focus on HF.
- We need to introduce a new Output option for computing the MP2 one-body density matrix
- After the SCF cycle, an HF calculation will give us the MO on the grid. At this point, if the user requested the MP2 one-body density matrix, call the routines that build it. See doc, [1]
- With reference to [1], Section Theoretical Background, the MP2 one-body density matrix is defined as
\gamma_{ab} = \frac{1}{2} \sum_{ijc} t_{ij}^{ac} t_{ij}^{bc}
, wherei, j
denote occupied orbitals,a, b
denote virtual (unoccupied) orbitals andt_{ij}^{ab} = \frac{\bra{ij}\ket{ab}}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b}
, where\epsilon_i
refer to the Hartree-Fock molecular orbital energies, and\bra{ij}\ket{ab}
is an antisymmetrized two-electron integral:- Physicists’ notation
\bra{ij}\ket{ab} = \int \phi_i^* \phi_j^* \frac{1}{r_{12}} \phi_a \phi_b
- Chemists’ notation
\left[ij|ab\right] = \int \phi_i^* \phi_j \frac{1}{r_{12}} \phi_a^* \phi_b
- SHALL WE SUPPORT BOTH NOTATIONS?
- Physicists’ notation
- The dimensions of
\gamma_{ab}
correspond to the number of unoccupied states in the HF calculation - As one can see from the formulas above, in order to compute the MP2 matrix, we need the eigenvalues of the HF solution and the two-body electron integrals. Since Octopus already has a routine to compute such elements, we will keep the MO states on the grid, and just use the existing routines. Therefore, the steps for this task are:
- Compute the two body electron integrals
- Build the MP2 matrix
- With reference to [1], Section Theoretical Background, the MP2 one-body density matrix is defined as
- Diagonalize the density matrix
-
\gamma V = nV
, where the eigenvectorsV
are the virtual NO and the eigenvaluesn
represent the occupation of each natural orbital. - Note that we are interested in only a few natural orbitals, so we can use an iterative eigensolver to reduce computation time. On the other hand,
\gamma
will not be big (at least for finite systems)
-
- At this stage, we need to build the new state object that will handle the natural orbitals states. In order to preserve the parallelization scheme and to use the Octopus routines, we duplicate the electronic states object (use
states_elec_copy
) - Populate the new states object: loop over the eigenvectors from the MP2 matrix and use them as coefficients to weight the HF orbitals in the grid. In this way, we have a new states state object which contains the natural orbitals on the grid
-
\Psi_j^{NO} = \sum_{i=0}^{n} V_{i,j}\phi_i
wheren
is the number of unoccupied HF states,V
is the eigenvector from the diagonalization of the MP2 matrix and\phi
is the unoccupied HF state
-
- Finally, call the routines for one-body and two-body matrix elements passing the new states object and output the matrix elements in natural orbital basis
Test
A hydrogen atom with one occupied and one unoccupied state -> no correlation, but everything is analytical so that we can test the routines
Validate one-body and two-body matrix element for k-points different than \Gamma
with PySCF
References
[1] Ashutosh Kumar and T. Daniel Crawford, Frozen Virtual Natural Orbitals for Coupled-Cluster Linear-Response Theory, J.Phys.Chem.A2017, 121, 708−716