What are Pair Natural Orbitals?

Question

Recently, I came across electron correlation methods that use local pair natural orbitals (LPNOs) or domain-based local pair natural orbitals (DLPNO), such as DLPNO-MP2, DLPNO-CCSD etc. The dynamic electron correlation is constrained to local areas and it reduces the computational cost. I tried to look at some papers such as this one, but I can't figure out what the PNOs actually are (or what they mean physically).

So, what are these PNOs? How are they different to the usual orbital localization schemes such as Ruedenberg, Pipek-Mezey? And most importantly, why are they used for correlation calculation? It would help if the answer is explained in simple words, I am not much of a theoretician!

Answer 1

I think Nike has answered to all questions adequately enough. I am sharing my understanding as one of the developers of the PNO-based local coupled-cluster (CC) methods, codes for computing response properties in particular.

In coupled-cluster theory, the correlated wavefunction is described in terms of "cluster amplitudes" (which are the wavefunction parameters). In the standard formulation of CC theory, these cluster amplitudes are defined in the basis of standard canonical Hartree-Fock (HF) orbitals. For large molecules, the number of cluster amplitudes defined in terms of HF orbitals may range up to a few billions (a crude estimate!). This makes the CC calculation formidable even on the most powerful and spacious computer on the planet. The only way to accomplish such calculations is to reduce the number of cluster amplitudes, in this way achieving a compact description of the correlated wavefunction.

Localized orbitals are used in quantum chemistry for various purposes. One of them is to get a compact description of the many-electron wavefunction. While the canonical HF orbitals spread over a large molecule as the whole, localized orbitals have much smaller spatial extents. Two types of orbital bases are needed to define the cluster amplitudes in CC theory: occupied orbitals and "virtual" or unoccupied orbitals. Pipek-Mezey or Foster-Boys schemes are used for obtaining only the localized occupied orbitals. The quest for a compact local description of the virtual orbitals is still open.

Several options have been put forth for defining a compact virtual orbital basis. I will not go into the details of all those. However, it is very important to mention that the pioneering idea in this field (local MP2 in particular) came from Prof. Peter Pulay, who proposed the use of the "Projected Atomic Orbitals"(PAOs). I will skip a description of them in this context.

The use of the "Pair-natural orbitals" (PNOs) is another option to compress the virtual space. The PNOs were introduced in the context of Coupled Electron-Pair Approximations (CEPA) by Wilfried Meyer in 1970's (both PNO-CEPA and PNO-CI) and was resurrected recently by Prof. Frank Neese and co-workers (including myself) in the context of CC theory, or more specifically, the domain-based pair-natural orbital CC (DLPNO-CC) approach. Several applications involving very large molecular systems (even a whole crambin protein!) has demonstrated that the PNOs provide the most compact description of the virtual orbital space, and the DLPNO-CC [both DLPNO-CCSD and DLPNO-CCSD(T)] method can truly achieve a linear scaling of the memory cost and the computational cost (in terms of wall times) with respect to system size.

The way PNOs are derived in the DLPNO-CC approach is rather involved. Here is a simple recipe to obtain the PNOs (this is absolutely not how it works for DLPNO-CC and the description below should only be taken as a conceptually simplified guideline):

1. For a large molecule, get the optimized HF-SCF MOs.
2. Localize the occupied orbitals using Pipek-Mezey or the Foster-Boys scheme.
3. Obtain MP2 guess for the cluster amplitudes using localized occupied orbitals and the canonical HF virtual orbitals.
4. Define "pair density" for each pair of localized occupied orbitals (i,j). These pair densities are defined solely in terms of virtual orbitals. (I will skip the formula, as quite frankly, I do not know how to write equations on this platform).
5. Diagonalize the pair density matrix. This gives the "pair-natural orbital" occupation numbers and the PNO coefficient vectors. (Note: The term "natural orbitals" in general is used for designating a set of orbitals that diagonalizes the one-particle density matrix, e.g., MCSCF natural orbitals. The name "pair-natural orbitals" also derives from the same concept. The PNOs diagonalize the pair density matrix for every pair of localized occupied orbitals.)
6. The final step is to expand the PNOs in terms of certain basis functions. The DLPNO-CC approach expands them in terms of the PAOs mentioned above.

How do the PNOs achieve a compact description of the virtual space? Once the pair density matrices are diagonalized, the PNO occupation numbers are compared to a user-defined threshold, TcutPNO, in the context of DLPNO-CC. All PNOs with occupation numbers less than TCutPNO are discarded. For each pair of localized occupied orbitals (i,j), one thus obtains a much fewer number of PNOs to describe the virtual space than the billions of canonical virtual HF MOs. The cluster amplitudes are defined for each pair (i,j) only in terms of the corresponding PNOs. This gives a highly compressed description of the correlated wavefunction.

I hope the descriptions above give some conceptual background and also addresses the questions. I did skip a lot of minute details, however.

Answer 2

I'll answer each of your three questions separately, but the one you say is "most important" will go first 😊

And most importantly, why are they used for correlation calculation?

They can significantly reduce the cost of a calculation on a big system, especially when there is a large number of "virtual" orbitals (unoccupied orbitals) in the basis set. Effectively they allow the size of the virtual space to be reduced. You mentioned MP2 and CCSD, which scale rapidly with the number $N$ of orbitals: $\mathcal{O}(N^5)$ for MP2 and $\mathcal{O}(N^6)$ for CCSD, so when $N$ is large (for example 4000 orbitals for a 40 atom system) it can become absolutely crucial to effectively reduce $N$ from a cost perspective. Without PNO based methods, it can be extremely difficult to do MP2 or CCSD on such a large number of atoms, even with a TZ basis set, but with LNO (similar to PNO) it was possible to do CCSD(T) on a molecule with 1023 atoms in a QZ basis set (44712 orbitals). For a small number of atoms (for example 10), in not too large a basis set (for example QZ) then PNO based methods are probably not worth the trouble and slight loss of accuracy which occurs in implementations of PNO-MP2 and PNO-CCSD.

What are these PNOs?

The term was first proposed in 1966 by Edmiston and Krauss as "pseudonatural orbitals" because, as Mayer described it in some context they can be considered an approximation to natural orbitals ("natural orbitals" are eigenvectors of the 1-electron density matrix), even though they can be very different from natural orbitals. Later people started to refer to them as "pair natural orbitals" instead of "pseudonatural orbitals" but even people that call them pair natural orbitals mean the same thing as Edmiston and Krauss did. Pair natural orbitals are eigenvectors of the "pair density matrix".

Since you said:

It would help if the answer is explained in simple words, I am not much of a theoretician!

I might be getting a bit too zealous by getting into more detail, but perhaps others will appreciate it. PNOs are eigenvectors of the density matrix for the "independent pair wavefunctions" (I'll use the notation in the aforementioned paper by Mayer):

$$ \tag{1} \Psi_0 + \sum_i \tilde{C}_P^{ai} \Phi_P^{ai} + \sum_{ij}\tilde{C}_P^{ij} \Phi_P^{ij}, $$

where $\Phi_P^{mn}$ is a Slater determinant (configuration) obtained by coupling two electrons with orbitals $m$ and $n$ with a double hole state $P$ (which is defined in the bottom-left corner of the 2nd page of Mayer's paper), and the coefficients $\tilde{C}$ variationally minimize the energy of $\Psi_P$.

In perhaps Frank Neese's earliest work on the subject (circa 2009) he and co-authors say:

"each electron pair is treated by the most rapidly converging expansion of external orbitals, that, by definition, is provided by the natural orbitals that are specific for this pair [76]",

where [76] is this 1955 paper by Lowdin.

How are they different to the usual orbital localization schemes such as Ruedenberg, Pipek-Mezey?

In the aforementioned paper by Neese et al. they say this right in the abstract:

"The internal space is spanned by localized internal orbitals. The external space is greatly compressed through the method of pair natural orbitals PNOs".

By "internal space" they mean occupied orbitals, and by "external space" they mean unoccupied orbitals. Basically: they localize the occupied orbitals by schemes such as the one by Pipek-Mezey, and they use PNO for the unoccupied orbitals.