Natural orbital functional for multiplets

TL;DR

This paper introduces a new natural orbital functional tailored for electronic multiplet states, accurately capturing electron correlation effects and conserving total spin, with promising results for transition-metal ionization potentials.

Contribution

A novel natural orbital functional for multiplets based on a new 2RDM reconstruction that preserves spin and captures electron correlation.

Findings

Accurately predicts ionization potentials of transition metals.

Fulfills (2,2)-positivity and spin conservation.

Matches high-level theoretical and experimental data.

Abstract

A natural orbital functional for electronic systems with any value of the spin is proposed. This energy functional is based on a new reconstruction for the two-particle reduced density matrix (2RDM) of the multiplet, that is, of the mixed quantum state that allows all possible spin projections. The mixed states of maximum spin multiplicity are considered. This approach differs from the methods routinely used in electronic structure calculations that focus on the high-spin component or break the spin symmetry. In the ensemble, there are no interactions between electrons with opposite spins in singly occupied orbitals, as well as new inter-pair alpha-beta-contributions are proposed in the 2RDM. The proposed 2RDM fulfills (2,2)-positivity necessary N-representability conditions and guarantees the conservation of the total spin. The NOF for multiplets is able to recover the non-dynamic and the intrapair dynamic electron correlation. The missing dynamic correlation is recovered by the NOF-MP2 method. Calculation of ionization potentials of the first-row transition-metal atoms is presented as test case. The values obtained agree with those reported at the coupled cluster singles and doubles level of theory with perturbative triples and experimental data.

Tables

Table 1. Table 1: Ionizations potentials (kcal/mol) using the cc-pVTZ basis set. CCSD(T) and experimental values taken from Table IV of Ref. (Balabanov and Peterson, 2005 ) .

Atom		X		X+		PNOF7	CCSD(T)	PNOF7-MP2	EXP
Sc		2D		3D		143.8	146.8	148.4	151.3
Ti		3F		4F		151.7	152.8	156.7	157.5
V		4F		5D		151.1	149.4	149.2	155.2
Cr		7S		6S		141.8	147.9	160.7	156.0
Mn		6S		7S		167.1	166.0	174.4	171.4
Fe		5D		6D		183.1	176.2	177.5	182.3
Co		4F		3F		182.8	173.3	165.7	181.5
Ni		3F		2D		187.8	164.5	166.7	175.1
Cu		2S		1S		178.9	167.0	167.6	178.2
Zn		1S		2S		189.0	208.0	208.3	216.6
MAE						7.9	7.3	6.5

Equations

$$E={\tilde{E}}_{hf}+E^{corr}={\tilde{E}}_{hf}+E^{dyn}+E^{sta}$$

$$\hat{H}=\sum _{ik}{\mathcal{H}}_{ki}{\hat{a}}_k^{†}{\hat{a}}_i+\frac{1}{2}\sum _{ijkl}⟨kl|ij⟩{\hat{a}}_k^{†}{\hat{a}}_l^{†}{\hat{a}}_j{\hat{a}}_i$$

$$\hat{\mathfrak{D}}=\sum _{M=-S}^S{\omega }_M|SM⟩⟨SM|$$

$$E=\sum _{ik}{\mathcal{H}}_{ki}{\Gamma }_{ki}+\sum _{ijkl}⟨kl|ij⟩D_{kl,ij}$$

$$\begin{array}[]{c}\Gamma_{ki}={\displaystyle{\textstyle{\displaystyle\sum_{M=-S}^{S}}}}\omega_{M}\left\langle SM\right|\hat{a}_{k}^{\dagger}\hat{a}_{i}\left|SM\right\rangle\\ D_{kl,ij}={\displaystyle{\textstyle{\displaystyle\frac{1}{2}\sum_{M=-S}^{S}}}}\omega_{M}\left\langle SM\right|\hat{a}_{k}^{\dagger}\hat{a}_{l}^{\dagger}\hat{a}_{j}\hat{a}_{i}\left|SM\right\rangle\end{array}$$

$$E=\sum _i{n}_i{\mathcal{H}}_{ii}+\sum _{ijkl}D[n_i,n_j,n_k,n_l]⟨kl|ij⟩$$

$$\theta_{k}=\sigma\left(1\right)\sigma\left(2\right)...\sigma\left(\mathrm{N_{I}}\right),\quad\sigma\left(i\right)=\left\{\begin{array}[]{c}\alpha\left(i\right)\\ \beta\left(i\right)\end{array}\right.$$

$${\hat{S}}_z{\theta }_k=\frac{1}{2}\left(\mu -\nu \right){\theta }_k$$

$$\alpha \left(1\right)\alpha \left(2\right),\alpha \left(1\right)\beta \left(2\right),\beta \left(1\right)\alpha \left(2\right),\beta \left(1\right)\beta \left(2\right)$$

$$\left[{\alpha }^2{\beta }^2\right]=\alpha \alpha \beta \beta +\alpha \beta \alpha \beta +\alpha \beta \beta \alpha +\beta \alpha \alpha \beta +\beta \alpha \beta \alpha +\beta \beta \alpha \alpha$$

$${\hat{S}}^2\left[{\alpha }^{\mu }{\beta }^{\nu }\right]=\frac{{\mathrm{N}}_{\mathrm{I}}}{2}\left(\frac{{\mathrm{N}}_{\mathrm{I}}}{2}+1\right)\left[{\alpha }^{\mu }{\beta }^{\nu }\right]$$

$${\hat{S}}_z\left[{\alpha }^{\mu }{\beta }^{\nu }\right]=\frac{1}{2}\left(\mu -\nu \right)\left[{\alpha }^{\mu }{\beta }^{\nu }\right]$$

$$\left\{\begin{array}[]{c}\qquad\qquad\alpha\left(1\right)\alpha\left(2\right),\qquad\qquad M=1\\ \frac{1}{\sqrt{2}}\left[\alpha\left(1\right)\beta\left(2\right)+\beta\left(1\right)\alpha\left(2\right)\right],\>M=0\\ \qquad\qquad\beta\left(1\right)\beta\left(2\right),\qquad\qquad M=-1\end{array}\right.$$

$$<{\hat{S}}_z>=\frac{1}{{\mathrm{N}}_{\mathrm{I}}+1}\sum _{M=-{\mathrm{N}}_{\mathrm{I}}/2}^{{\mathrm{N}}_{\mathrm{I}}/2}M=0$$

$$\sum _{p\in {\Omega }_g}{n}_p=1,\Omega {}_g\in \Omega {}_{\mathrm{II}}$$

$$2\sum _{p\in {\Omega }_{\mathrm{II}}}{n}_p=2\sum _{g=1}^{{\mathrm{N}}_{\mathrm{II}}/2}\sum _{p\in {\Omega }_g}{n}_p={\mathrm{N}}_{\mathrm{II}}$$

$$2\sum _{p\in {\Omega }_{\mathrm{I}}}{n}_p=2\sum _{g={\mathrm{N}}_{\mathrm{II}}/2+1}^{{\mathrm{N}}_{\mathrm{\Omega }}}{n}_g={\mathrm{N}}_{\mathrm{I}},$$

$$2\sum _{p\in \Omega }{n}_p=2\sum _{p\in {\Omega }_{\mathrm{II}}}{n}_p+2\sum _{p\in {\Omega }_{\mathrm{I}}}{n}_p={\mathrm{N}}_{\mathrm{II}}+{\mathrm{N}}_{\mathrm{I}}=\mathrm{N}$$

$$\begin{array}[]{c}D_{pq,rt}^{\alpha\beta}={\displaystyle\frac{\Pi_{pr}}{2}}\delta_{pq}\delta_{rt}\delta_{p\Omega_{g}}\delta_{r\Omega_{g}}\;(g\leq\frac{N_{\mathrm{II}}}{2})\\ \\ \Pi_{pr}=\left\{\begin{array}[]{c}\sqrt{n_{p}n_{r}}\qquad p=r\textrm{ or }p,r>\frac{N_{\mathrm{II}}}{2}\\ -\sqrt{n_{p}n_{r}}\qquad p=g\textrm{ or }r=g\end{array}\right.\end{array}$$

$$D_{pq,rt}^{\sigma \sigma }=\frac{n_p{n}_q}{2}\left({\delta }_{pr}{\delta }_{qt}-{\delta }_{pt}{\delta }_{qr}\right){\delta }_{p{\Omega }_f}{\delta }_{q{\Omega }_g},{}_{\sigma =\alpha ,\beta }$$

$$\begin{array}[]{c}D_{pq,rt}^{\alpha\beta}={\displaystyle{\displaystyle\frac{n_{p}n_{q}}{2}}}\delta_{pr}\delta_{qt}\delta_{p\Omega_{f}}\delta_{q\Omega_{g}}-{\displaystyle\frac{\Phi_{p}\Phi_{r}}{2}}\delta_{p\Omega_{f}}\delta_{r\Omega_{g}}\cdot\\ \\ \cdot\left\{\begin{array}[]{c}\delta_{pq}\delta_{rt}\quad f\leq\frac{\mathrm{N_{II}}}{2}\textrm{ or }g\leq\frac{\mathrm{N_{II}}}{2}\\ \delta_{pt}\delta_{qr}\qquad\frac{\mathrm{N_{II}}}{2}<f,g\leq\mathrm{N}_{\Omega}\end{array}\right.\end{array}$$

$$\begin{array}[]{c}\mathrm{<}\hat{S}^{2}\mathrm{>}={\displaystyle\frac{\mathrm{N}\left(4-\mathrm{N}\right)}{4}}{\displaystyle+\sum\limits_{pq}}\left\{D_{pq,pq}^{\alpha\alpha}+D_{pq,pq}^{\beta\beta}-2D_{pq,qp}^{\alpha\beta}\right\}\end{array}$$

$$\sum _{pq}\to \sum _{g=1}^{{\mathrm{N}}_{\Omega }}\sum _{p,q\in {\Omega }_g}+\sum _{f\ne g=1}^{{\mathrm{N}}_{\Omega }}\sum _{p\in {\Omega }_f}\sum _{q\in {\Omega }_g}$$

$$\begin{array}[]{c}\sum\limits_{pq}D_{pq,pq}^{\sigma\sigma}={\displaystyle\frac{1}{2}}{\displaystyle\sum_{f\neq g=1}^{\mathrm{N}_{\Omega}}\;\sum_{p\in\Omega_{f}}\sum_{q\in\Omega_{g}}}{\displaystyle n_{p}n_{q}}={\displaystyle\frac{1}{2}\left\{{\displaystyle\sum_{f\neq g=1}^{\frac{\mathrm{N_{II}}}{2}}}+\right.}\\ \\ {\displaystyle\left.\sum_{f=1}^{\frac{\mathrm{N_{II}}}{2}}\;\sum_{g=\frac{\mathrm{N_{II}}}{2}+1}^{\mathrm{N}_{\Omega}}+\sum_{f=\frac{\mathrm{N_{II}}}{2}+1}^{\mathrm{\mathrm{N}_{\Omega}}}\;\sum_{g=1}^{\frac{\mathrm{N_{II}}}{2}}+\sum_{f\neq g=\frac{\mathrm{N_{II}}}{2}+1}^{\mathrm{N}_{\Omega}}\right\}}{\displaystyle\sum_{p\in\Omega_{f}}\sum_{q\in\Omega_{g}}}{\displaystyle n_{p}n_{q}}\\ \\ ={\displaystyle\frac{1}{2}\left\{{\displaystyle\sum_{f\neq g=1}^{\frac{\mathrm{N_{II}}}{2}}}1+\sum_{f=1}^{\frac{\mathrm{N_{II}}}{2}}\;\sum_{g=\frac{\mathrm{N_{II}}}{2}+1}^{\mathrm{N}_{\Omega}}\frac{1}{2}+\sum_{f=\frac{\mathrm{N_{II}}}{2}+1}^{\mathrm{\mathrm{N}_{\Omega}}}\;\sum_{g=1}^{\frac{\mathrm{N_{II}}}{2}}\frac{1}{2}\;+\right.}\\ \\ \left.{\displaystyle\sum_{f\neq g=\frac{\mathrm{N_{II}}}{2}+1}^{\mathrm{N}_{\Omega}}\frac{1}{4}}\right\}={\displaystyle\frac{\mathrm{N_{II}\left(N_{II}-2\right)}}{8}+\frac{\mathrm{N_{II}N_{I}}}{4}+\frac{\mathrm{N_{I}\left(N_{I}-1\right)}}{8}}\end{array}$$

$$\sum _{p\in {\Omega }_g}{n}_p=n_g=1/2,\Omega {}_g\in \Omega {}_{\mathrm{I}}$$

$$\begin{array}[]{c}\sum\limits_{pq}D_{pq,qp}^{\alpha\beta}={\displaystyle\frac{1}{2}}{\displaystyle\sum_{g=1}^{\frac{\mathrm{N_{II}}}{2}}\;}{\displaystyle\sum_{p\in\Omega_{g}}}{\displaystyle n_{p}}-{\displaystyle\frac{1}{2}}{\displaystyle\sum_{f\neq g=\frac{\mathrm{N_{II}}}{2}+1}^{\mathrm{N}_{\Omega}}\;\sum_{p\in\Omega_{f}}\sum_{q\in\Omega_{g}}}{\displaystyle n_{p}n_{q}}\\ ={\displaystyle\frac{1}{2}}\left\{{\displaystyle{\displaystyle\sum_{g=1}^{\frac{\mathrm{N_{II}}}{2}}}\,1-\sum_{f\neq g=\frac{\mathrm{N_{II}}}{2}+1}^{\mathrm{N}_{\Omega}}\frac{1}{4}}\right\}={\displaystyle\frac{\mathrm{N_{II}}}{4}}-{\displaystyle\frac{\mathrm{N_{I}\left(N_{I}-1\right)}}{8}}\end{array}$$

$$\begin{array}[]{c}\mathrm{<}\hat{S}^{2}\mathrm{>}={\displaystyle\frac{\mathrm{\left(N_{I}+N_{II}\right)}\left(4-\mathrm{N_{I}-N_{II}}\right)}{4}}\\ \\ +2\left[{\displaystyle\frac{\mathrm{N_{II}\left(N_{II}-2\right)}}{8}+\frac{\mathrm{N_{II}N_{I}}}{4}+\frac{\mathrm{N_{I}\left(N_{I}-1\right)}}{8}}\right]\\ \\ {\displaystyle-\frac{\mathrm{N_{II}}}{2}}+{\displaystyle\frac{\mathrm{N_{I}\left(N_{I}-1\right)}}{4}}={\displaystyle\frac{\mathrm{N_{I}}}{2}\left(\frac{\mathrm{N_{I}}}{2}+1\right)}\end{array}$$

$$E=\sum _{g=1}^{{\mathrm{N}}_{\mathrm{II}}/2}{E}_g+\sum _{g={\mathrm{N}}_{\mathrm{II}}/2+1}^{{\mathrm{N}}_{\Omega }}{\mathcal{H}}_{gg}+\sum _{f\ne g}^{{\mathrm{N}}_{\Omega }}{E}_{fg}$$

$$E_g=2\sum _{p\in {\Omega }_g}{n}_p{\mathcal{H}}_{pp}+\sum _{q,p\in {\Omega }_g}{\Pi }_{qp}{\mathcal{K}}_{pq},\Omega {}_g\in {\Omega }_{\mathrm{II}}$$

$$E_{fg}=\sum _{p\in {\Omega }_f}\sum _{q\in {\Omega }_g}\left[n_q{n}_p\left(2{\mathcal{J}}_{pq}-{\mathcal{K}}_{pq}\right)-{\Phi }_q{\Phi }_p{\mathcal{K}}_{pq}\right]$$

Full text

I Introduction

In the absence of an external field at absolute zero temperature, the non-relativistic Hamiltonian used in electronic calculations does not contain spin coordinates. The eigenvectors of such a Hamiltonian are eigenvectors of the square of the total spin ($\hat{S}^{2}$) and of one of its components, usually chosen to be $\hat{S}_{z}$. Consequently, the ground state of a many-electron system with total spin $S$ is a multiplet, i.e., a mixed quantum state that allows all possible $S_{z}$ values.

Despite its success, the Kohn-Sham density functional theory (Kohn and Sham, 1965), including the spin-dependent-density formalism, can not describe the degeneracy of the spin-multiplet components. The proper description of the ensemble of pure spin states is provided by the one-particle reduced density matrix (1RDM) functional theory (Gilbert, 1975). In 1980 (Valone, 1980), Valone extended Levy’s functional (Levy, 1979) to include all ensemble N-representable 1RDMs in its domain. A series of conditions fulfilled by this exact functional have been recently investigated (Rohr and Pernal, 2011). Unfortunately, computational schemes based on the constrained search formulation are several times more expensive than solving directly the Schrödinger equation; therefore, the construction of the functional requires a practical approach.

Actually, we must only construct the 1RDM functional for the electron-electron potential energy $V_{ee}$, since the non-interacting part of the electronic Hamiltonian is a one-particle operator. $V_{ee}$ is a well-known functional of the two-particle reduced density matrix (2RDM), hence a convenient approach is to reconstruct the 2RDM in terms of the 1RDM. However, this does not fully reconstruct the ground-state energy, only the 2RDM. Approximating the energy implies that theorems established for the exact 1RDM functional may not be fulfilled since the 2RDM dependence remains (Donnelly, 1979). To guarantee the existence of an N-electron system compatible with a proposed functional, we must observe the N-representability conditions (Mazziotti, 2012) on the reconstructed 2RDM. Therefore, we are dealing with approximate 1RDM methods, where the 2RDM continues to play a dominant, albeit hidden, role (Piris, 2018a).

In most applications, the spectral decomposition of the 1RDM is used to express it in terms of the naturals orbitals (NOs) and their occupation numbers (ONs). In this representation, the energy expression is referred to as natural orbital functional (NOF). Within NOF approximations, spin-dependent formalisms have been also proposed (Goedecker and Umrigar, 1998; Lathiotakis et al., 2005), which do not conserve spin. NOFs that correctly reproduce the expectation values of spin operators have been reported (Leiva and Piris, 2007; Piris et al., 2009, 2010a; Quintero-Monsebaiz et al., 2019), however, only for the high-spin component of the multiplet.

The aim of this work is to propose a NOF that allows to describe an electronic system with any value of $S$, that is, a new reconstruction of the 2RDM for spin-multiplets. For the latter, the mixed states of maximum multiplicity are considered.

The chemical versatility of transition metal (TM) atoms is derived from their ability to exist in multiple spin states, so there is considerable interest in having the ability to predict reliable properties for these systems. One of these properties is the ionization potential (IP). The correct description of TM IP is a great challenge since not only a significant rearrangement of the electronic configurations takes place, but the resulting contraction of the electron cloud implies a considerable contribution of the dynamic correlation. Hence, multi-reference methods are needed for an accurate description (Thomas et al., 2015). Unfortunately, multireference methods are more complicated in implementation and much less easy to use than single reference methods.

Recently (Piris, 2017, 2018b), a single-reference method for the electron correlation was introduced taking as reference the Slater determinant formed with the NOs of an approximate NOF. In this approach, called natural orbital functional - second-order Møller–Plesset (NOF-MP2) method, the total energy of an N-electron system can be attained by the expression

[TABLE]

where $\tilde{E}_{hf}$ is the Hartree-Fock (HF) energy obtained with the NOs, the dynamic energy ($E^{dyn}$) is derived from a modified MP2 perturbation theory, while the non-dynamic energy ($E^{sta}$) is obtained from the static component of the employed NOF. An important feature of the method is that double counting is avoided by taking the amount of static and dynamic correlation in each orbital as a function of its occupancy.

The second aim of the present work is to extend the NOF-MP2 method to the spin multiplets and apply it to the description of the first-row TM IPs.

This article is organized as follows. We start in section II with the basic concepts and notations related to the NOF theory. Section III is devoted to the extension of the NOF-MP2 method to spin multiplets. Finally, I present the results obtained for the IPs of the TM Sc-Zn elements in Section IV.

II Natural Orbital Functional Theory

Consider the N-electron Hamiltonian

[TABLE]

where $\mathcal{H}_{ki}$ denote the matrix elements of the one-particle part of the Hamiltonian involving the kinetic energy and the potential energy operators, $\left\langle kl|ij\right\rangle$ are the two-particle interaction matrix elements, $\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$ are the familiar fermion creation and annihilation operators associated with the complete orthonormal spin-orbital set $\left\{\left|i\right\rangle\right\}$. For a given spin $S$, there are $\left(2S+1\right)$ energy degenerate eigenvectors $\left|SM\right\rangle$, so a mixed state is defined by the N-particle density matrix statistical operator

[TABLE]

In Eq. (3), $\omega_{M}$ are positive real numbers that sum one, so that $\mathfrak{\hat{D}}$ corresponds to a weighted sum of all accessible pure states. For equiprobable pure states, we take $\omega_{M}=(2S+1)^{-1}$.

The expectation value of (2) reads as

[TABLE]

where the 1RDM and 2RDM are

[TABLE]

Normalization of Löwdin is used, in which the traces of the 1RDM and 2RDM are equal to the number of electrons and the number of electron pairs, respectively. The last term in Eq. (4) is $V_{ee}$, an explicit functional of the 2RDM. To construct the functional $V_{ee}\left[\Gamma\right]$, we employ the representation where the 1RDM is diagonal ($\hat{\Gamma}={\displaystyle{\textstyle\sum}_{i}n_{i}}\left|i\right\rangle\left\langle i\right|$). Restriction on the ONs to the range $0\leq n_{i}\leq 1$ represents a necessary and sufficient condition for ensemble N-representability of the 1RDM Coleman (1963). This leads to a NOF, namely:

[TABLE]

where $D[n_{i},n_{j},n_{k},n_{l}]$ represents the reconstructed ensemble 2RDM from the ONs. For $\hat{S}_{z}$ eigenvectors, density matrix blocks that conserve the number of each spin type are non-vanishing, however, only three of them are independent, namely $D^{\alpha\alpha}$, $D^{\alpha\beta}$, and $D^{\beta\beta}$ (Piris, 2007).

Consider $\mathrm{N_{I}}$ single electrons and $\mathrm{N_{II}}$ paired electrons, so that $\mathrm{N_{I}}+\mathrm{N_{II}}=\mathrm{N}$. Assume further that all spins corresponding to $\mathrm{N_{II}}$ electrons are coupled as a singlet, thence the $\mathrm{N_{I}}$ electrons determine the spin of the system. In fact, excluding the singlet part, the $\hat{S}_{z}$ eigenfunctions $\theta_{k}$ are simply products of $\mathrm{N_{I}}$ one-electron functions ($\alpha$’s or $\beta$’s):

[TABLE]

[TABLE]

where $\mu$ is the number of $\alpha$’s, and $\nu$ is the number of $\beta$’s, respectively. For a given $\mu$, we have \left(\begin{array}[]{c}\mathrm{N_{I}}\\ \mu\end{array}\right)=\left(\begin{array}[]{c}\mathrm{N_{I}}\\ \nu\end{array}\right) primitive spin functions (7) which belong to the same eigenvalue of $\hat{S}_{z}$. For instance, we have four primitive spin functions in the two-electron case, namely,

[TABLE]

Let us define the square brackets of spins functions, $\left[\alpha^{\mu}\beta^{\nu}\right]$, as the sum of all primitive spin functions with $\mu$ $\alpha$-spin functions and $\nu$ $\beta$-spin functions. Example:

[TABLE]

It is not difficult to demonstrate (Paunez, 2000) that the square brackets are simultaneous eigenfunctions of $\hat{S}^{2}$ and $\hat{S}_{z}$ operators:

[TABLE]

[TABLE]

Note that the eigenvalue $M=\left(\mu-\nu\right)/2$ can take values in the range $-\mathrm{N_{I}}/2\leq M\leq\mathrm{N_{I}}/2$. All the square brackets are simultaneous eigenfuntions of $\hat{S}^{2}$ and $\hat{S}_{z}$ with the quantum numbers $\mathrm{N_{I}}/2$ and $M$, respectively. This is the mixed state of highest multiplicity: $2S+1=\mathrm{N_{I}}+1$, and it is the only state that belongs to the quantum number $S=\mathrm{N_{I}}/2$ (Paunez, 2000). For example, for $\mathrm{N_{I}}=2$ the triplet state ($S=1$) is

[TABLE]

Interestingly, the expected value of $\hat{S}_{z}$ for the whole ensemble is zero,

[TABLE]

In the absence of single electrons ($\mathrm{N_{I}}=0$), the energy must obviously be reduced to a NOF that well describes singlet states. Let us employ the electron-pairing approach in NOF theory (Piris, 2018c) for this purpose. First, the spin-restricted theory is adopted, in which a single set of spatial orbitals $\left\{\left|p\right\rangle\right\}$ is used for $\alpha$ and $\beta$ spins. Accordingly, all spatial orbitals are double occupied in the ensemble, so that occupancies for particles with $\alpha$ and $\beta$ spins are equal: $n_{p}^{\alpha}=n_{p}^{\beta}=n_{p}.$

Let us divide in turn the orbital space $\Omega$ into two subspaces: $\Omega=\Omega_{\mathrm{I}}\oplus\Omega_{\mathrm{II}}$. The orbital space $\Omega_{\mathrm{II}}$ is composed of $\mathrm{N_{II}}/2$ mutually disjoint subspaces: $\Omega_{\mathrm{II}}=\Omega_{1}\oplus\Omega_{2}\oplus...\oplus\Omega_{\mathrm{N_{II}}/2}$, $\Omega_{f}\cap\Omega_{g}=\textrm{\O}$ (Piris, 2018c). Each subspace $\Omega{}_{g}\in\Omega_{\mathrm{II}}$ contains one orbital $g$ below the level $\mathrm{N_{II}}/2$, and $\mathrm{N}_{g}$ orbitals above it, so $dim\left\{\Omega_{g}\right\}=\mathrm{N}_{g}+1$. In this work, $\mathrm{N}_{g}$ is equal to a fixed number corresponding to the maximum value allowed by the basis set used in calculations.

Taking into account the spin, the total occupancy for a given subspace $\Omega{}_{g}$ is 2, which is reflected in additional sum rule, namely,

[TABLE]

It follows that

[TABLE]

Similarly, $\Omega_{\mathrm{I}}$ is composed of $\mathrm{N_{I}}$ mutually disjoint subspaces: $\Omega_{\mathrm{I}}=\Omega_{\mathrm{N_{II}}/2+1}\oplus...\oplus\Omega_{\mathrm{N_{\Omega}}}$. Here, $\mathrm{\mathrm{N}_{\Omega}=}\mathrm{N_{II}}/2+\mathrm{N_{I}}$ denotes the total number of suspaces in $\Omega$. In contrast to $\Omega_{\mathrm{II}}$, each subspace $\Omega{}_{g}\in\Omega_{\mathrm{I}}$ contains only one orbital $g$ with $\mathrm{N_{II}}/2<g\leq N_{\Omega}$, i.e., $dim\left\{\Omega_{g}\right\}=1$. Besides, we take $2n_{g}=1$, so each orbital is fully occupied individually, but we do not know whether the electron has $\alpha$ or $\beta$ spin: $n_{g}^{\alpha}=n_{g}^{\beta}=n_{g}=1/2$. It follows that

[TABLE]

so taking into account Eq. (16), the trace of the 1RDM is verified equal to the number of electrons:

[TABLE]

It is worth noting that the splitting of the orbital space into the subspaces $\Omega_{\mathrm{I}}$ and $\Omega_{\mathrm{II}}$ is determined by the value of the total spin $S$. The latter determines the number of orbitals in the subspace $\Omega_{\mathrm{I}}$, which is equal to the number of unpaired electrons. If a system has $S=1/2$, only one orbital composes the subspace $\Omega_{\mathrm{I}}$, analogously, if a system has $S=1$, we have two orbitals in $\Omega_{\mathrm{I}}$, for $S=3/2$, we have three orbitals in $\Omega_{\mathrm{I}}$, and so on. The rest of the orbitals can be distributed in the subspace $\Omega_{\mathrm{II}}$ that constitutes a singlet, and to which this division is reduced in the case $S=0$ as expected.

The electron-pairing approach is used to divide the subspace $\Omega_{\mathrm{II}}$ because of the efficiency it has demonstrated in the description of singlet states (Piris, 2018c), but the method proposed here can be implemented to extent other approximate NOFs used in singlets. It is important to remember that the pairing scheme of the orbitals is allowed to vary along the optimization process until the most favorable orbital interactions are found. Therefore, the orbitals do not remain fixed in the optimization process, they adapt to the problem. In Fig. 1, an illustrative example is shown. In this example, $S=1$ and $\mathrm{N_{I}}=2$, so two orbitals make up the subspace $\Omega_{\mathrm{I}}$, whereas six electrons ($\mathrm{N_{II}}=6$) distributed in three subspaces $\left\{\Omega_{1},\Omega_{2},\Omega_{3}\right\}$ make up the subspace $\Omega_{\mathrm{II}}$. The maximum value allowed by the basis set is $\mathrm{N}_{g}=4$.

The reconstruction of the 2RDM must be subject to N-representability conditions. In general, $D$ has a dependence on four indexes, computationally expensive. We shall employ the two-index reconstruction proposed in reference (Piris, 2006). This particular reconstruction is based on the introduction of two auxiliary matrices $\Delta$ and $\Pi$ expressed in terms of the ONs to reconstruct the cumulant part of the 2RDM. The necessary N-representability D, Q, and G conditions of the 2RDM, also known as (2,2)-positivity conditions (Mazziotti, 2012), impose strict inequalities on the off-diagonal elements of $\Delta$ and $\Pi$ matrices (Piris et al., 2010b). Appropriate forms of $\Delta$ and $\Pi$ have led to several $\mathcal{JKL}$-only functionals. $\mathcal{J}$ and $\mathcal{K}$ refer to the usual Coulomb and exchange integrals, $\left\langle pq|pq\right\rangle$ and $\left\langle pq|qp\right\rangle$ respectively, while $\mathcal{L}$ denotes the exchange-time-inversion integral (Piris, 1999) $\left\langle pp|qq\right\rangle$.

Let us extent the last member of this functional family, PNOF7 (Piris, 2017; Mitxelena et al., 2018; Piris, 2018b), to multiplets. Accordingly, divide the matrix elements of $D$ into intra- and inter-subspace contributions. For intra-subspace blocks, the functional form is maintained, that is, only intrapair $\alpha\beta$-contributions appear,

[TABLE]

In (19), the Kronecker delta has an obvious meaning, for instance, $\delta_{p\Omega_{g}}=1$ if $p\in\Omega_{g}$ or $\delta_{p\Omega_{g}}=0$ otherwise. Note also that $\Omega{}_{g}\in\Omega_{\mathrm{II}}$. A major change is observed with respect to the $D^{\alpha\beta}$ of a singlet, namely: $D_{pp,pp}^{\alpha\beta}=0$, $\forall p\in\Omega_{\mathrm{I}}$. Actually, there can be no interactions between electrons with opposite spins in a singly occupied orbital, since there is only one electron with $\alpha$ or $\beta$ spin in each $\left|SM\right\rangle$ of the ensemble.

A different ansatz is needed for inter-subspace contributions ($\Omega_{f}\neq\Omega{}_{g}$). The spin-parallel matrix elements remain Hartree-Fock like,

[TABLE]

whereas the spin-antiparallel blocks are

[TABLE]

where $\Phi_{p}=\sqrt{n_{p}(1-n_{p})}$. It should be noted that the second term in Eq. (21) differs from our previous proposition for singlets (Piris, 2017; Mitxelena et al., 2018). As a result, inter-subspace $\alpha\beta$-contributions involving $\Phi$ in subspace $\Omega_{\mathrm{I}}$, lead to exchange integrals instead of exchange-time-inversion integrals.

It is not difficult to verify that the reconstruction (19)-(21) leads to $\mathrm{<}\hat{S}^{2}\mathrm{>}=\frac{\mathrm{N_{I}}}{2}\left(\frac{\mathrm{N_{I}}}{2}+1\right)$. In fact, the expectation value of the operator $\hat{S}^{2}$ reads as (Piris et al., 2009)

[TABLE]

Note that

[TABLE]

According to Eq. (20), the spin-parallel matrix elements of the 2RDM have only inter-subspace contributions ($\Omega_{f}\neq\Omega{}_{g}$). Consequently, the trace of $D^{\sigma\sigma}$ can be cast as

[TABLE]

where we have considered Eq. (15), and

[TABLE]

Taking into account the $\alpha\beta$-blocks of the 2RDM in Eqs. (19) and (21), the sum in the last $\alpha\beta$-term of Eq. (22) can be cast as

[TABLE]

where it has been considered that $\Phi_{p}=\sqrt{n_{p}(1-n_{p})}=n_{p}=1/2$ for $\forall p\in\Omega_{\mathrm{I}}$.

Combining Eq. (22) with Eqs. (23) and (25), one arrives at the ensemble average of the square of the total spin, namely:

[TABLE]

This result agrees with a multiplet of a total spin $S=\mathrm{N_{I}}/2$. The corresponding multiplicity, $2S+1=\mathrm{N_{I}}+1$, is the highest multiplicity for $\mathrm{N_{I}}$ electrons. A mixed state associated to our approximate 2RDM reconstruction will be non-degenerate (Paunez, 2000).

In order to simplify the derivation, let us assume real spatial orbitals, then $\mathcal{L}_{pq}=\mathcal{K}_{pq}$. After a simple algebra, the energy (6) can be written as

[TABLE]

where

[TABLE]

is the energy of a pair of electrons with opposite spins. It reduces to a NOF obtained from ground-state singlet wavefunction, so $E_{g}$ describes accurately two-electron systems (Piris, 2018c). In the last term of Eq. (27), $E_{fg}$ correlates the motion of electrons with parallel and opposite spins belonging to different subspaces ($\Omega_{f}\neq\Omega{}_{g}$):

[TABLE]

The functional (27)-(29) is PNOF7 for multiplets. The solution is established by optimizing the energy (27) with respect to the ONs and to the NOs, separately. The conjugate gradient method is used for performing the optimization of the energy with respect to auxiliary variables that enforce automatically the N-representability bounds of the 1RDM. The self-consistent procedure proposed in (Piris and Ugalde, 2009) yields the NOs by an iterative diagonalization procedure, in which orbitals are not constrained to remain fixed along the orbital optimization process.

III The NOF-MP2 method

The weakness of approach (21) is the absence of the inter-subspace dynamic electron correlation since $\Phi_{p}=\sqrt{n_{p}(1-n_{p})}$ has significant values only when the ONs differ substantially from 1 and 0. Consequently, PNOF7 is able to recover the complete intra-subspace, but only the static inter-subspace correlation. To add the missing dynamic correlation, we can use the NOF-MP2 method, Eq. (1).

The extension of NOF-MP2 to multiplets is straightforward and does not require substantial changes with respect to the method implemented for singlets (Piris, 2018b). PNOF7 provides the reference NOs to form $\tilde{E}_{hf}$. However, one has to take into account the presence of the singly occupied spatial orbitals. In this sense, the zeroth-order Hamiltonian for the modified MP2 is constructed from a closed-shell-like Fock operator that contains a HF density matrix with doubly ($2n_{g}=2$) and singly ($2n_{g}=1$) occupied orbitals. As a result, the reference energy $\tilde{E}_{hf}$ is

[TABLE]

In Eq. (30), the last term eliminates the $\alpha\beta$-contribution of the singly occupied orbitals to the energy because in each pure state $\left|SM\right\rangle$ of the ensemble there is no such interaction.

$E^{sta}$ is the sum of the static intra-space and inter-space correlation energies:

[TABLE]

where $\Lambda_{p}=1-\left|1-2n_{p}\right|$ is the amount of intra-space static correlation in each orbital as a function of its occupancy. Note that $\Lambda_{p}$ goes from zero for empty or fully occupied orbitals to one if the orbital is half occupied.

$E^{dyn}$ is obtained from the second-order correction $E^{\left(2\right)}$ of the MP2 method. The first-order wavefunction is a linear combination of all doubly excited configurations, considering one electron with $\alpha$ or $\beta$ spin in $\Omega_{\mathrm{I}}$. The dynamic energy correction takes the form

[TABLE]

where

[TABLE]

and $N_{B}$ is the number of basis functions. The amplitudes $T_{pq}^{fg}$ are obtained by solving the modified equations for the MP2 residuals (Piris, 2018b).

In order to avoid double counting of the electron correlation, the amount of dynamic correlation in each orbital $p$ is defined by functions $C_{p}$ of its occupancy, namely,

[TABLE]

where $C_{p}$ is divided into intra-space ($C_{p}^{tra}$) and inter-space ($C_{p}^{ter}$) contributions. According to Eq. (34), fully occupied and empty orbitals yield a maximal contribution to dynamic correlation, whereas orbitals with half occupancies contribute nothing. It is worth noting that $C_{p}^{ter}$ is not considered if the orbital is below $\mathrm{N_{\Omega}}$. Using these functions as the case may be (intra-space or inter-space), the modified off-diagonal elements of the Fock matrix ($\tilde{\mathcal{F}}$) are defined as

[TABLE]

as well as modified two-electron integrals:

[TABLE]

where the subspace index $g=1,...,\mathrm{N}_{\Omega}$. This leads to the following linear equation for the modified MP2 residuals:

[TABLE]

[TABLE]

where $i,j,k$ refer to the strong occupied NOs, and $a,b,c$ to weak occupied ones. It should be noted that diagonal elements of the Fock matrix ($\mathcal{F}$) are not modified. By solving this linear system of equations the amplitudes $T_{pq}^{fg}$ are obtained, which are inserted into the Eq. (32) to achieve $E^{dyn}$.

IV Ionization potentials of first-row transition metal atoms

The ability to exist in multiple spin states makes transition-metal (TM) atoms ideal candidates for testing PNOF7 in multiplets. Fig. 2 shows the calculated ionization potentials (IPs) of the TM elements Sc-Zn using the correlation-consistent valence triple-basis set including polarization (cc-pVTZ) (Balabanov and Peterson, 2005). The IPs are calculated by the energy difference between the positive ions ($\mathrm{X^{+}}$) and the neutral atoms (X): $\mathrm{IP}=E(\mathrm{X^{+}})-E(\mathrm{X})$.

For comparison, the IPs calculated using the same basis set at the coupled cluster singles and doubles level of theory with perturbative triples [CCSD(T)/cc-pVTZ], and the experimental data reported in Table IV of Ref. (Balabanov and Peterson, 2005), have been included. The data sets for these graphs can be found in Table 1. It is worth noting that it is not intended to reproduce the experimental data in this work, since it requires large basis sets and the inclusion of relativistic effects.

There are small amounts of spin-contamination in the solution of the CCSD equations (Balabanov and Peterson, 2005). Nevertheless, let us take the CCSD(T) values as benchmark data. The inspection of Fig. 2 reveals that calculated PNOF7 IPs from Sc to Mn are close to the CCSD(T) values, but deviate from those for the late transition metals, reaching deviations of 23 and 19 kcal/mol for Ni and Zn, respectively. Note that the mean absolute error (MAE) with respect to the experiment is quite similar for both methods, 7.9 and 7.6 kcal/mol respectively, although the behavior of PNOF7 is somewhat erratic due to the aforementioned systems.

The results obtained for the IPs of the TM elements Sc-Zn at the PNOF7-MP2 level of theory have also been included in Fig. 2. It is important to note that in the case of the neutral atoms Co, Ni, and Cu, as well as in the Co and Ni ions, there are two solutions. One of these solutions is the most favorable energetically at the PNOF7 level, whereas the other solution corresponds to the lowest in energy for the NOF-MP2 method. In the former, the non-dynamic character of the electron correlation predominates, while the dynamic correlation predominates in the latter. This reflects the complexity of the problem in these systems where the 3d shell is almost filled.

The data reveal a remarkable improvement, especially in the Cr, Ni and Zn atoms. The inclusion of inter-subspace dynamic electron correlation yields a MAE of 6.5 kcal/mol, which is an outstanding result considering the size of the basis sets employed. Large basis sets could improve the MAE to values lower than 1 kcal/mol (Balabanov and Peterson, 2005).

Acknowledgments: Financial support comes from MCIU/AEI/FEDER, UE (PGC2018-097529-B-100). The author thanks for technical and human support provided by IZO-SGI SGIker of UPV/EHU and European funding (ERDF and ESF).