Natural orbital functional for spin-polarized periodic systems

TL;DR

This paper extends natural orbital functional theory to spin-polarized periodic systems, enabling accurate modeling of strong electron correlation in systems like Hubbard models and hydrogen rings.

Contribution

It introduces an extension of the PNOF5 and PNOF7 models to describe spin-uncompensated systems with explicit handling of high-spin cases.

Findings

PNOF7 accurately describes strong correlation in hydrogen rings.

The reconstructed two-particle density matrix maintains N-representability and spin conservation.

Model results agree well with exact solutions for test systems.

Abstract

Natural orbital functional theory is considered for systems with one or more unpaired electrons. An extension of the Piris natural orbital functional (PNOF) based on electron pairing approach is presented, specifically, we extend the independent pair model, PNOF5, and the interactive pair model PNOF7 to describe spin-uncompensated systems. An explicit form for the two-electron cumulant of high-spin cases is only taken into account, so that singly occupied orbitals with the same spin are solely considered. The rest of electron pairs with opposite spins remain paired. The reconstructed two-particle reduced density matrix fulfills certain N-representability necessary conditions, as well as guarantees the conservation of the total spin. The theory is applied to model systems with strong non-dynamic (static) electron correlation, namely, the one-dimensional Hubbard model with periodic boundary conditions and hydrogen rings. For the latter, PNOF7 compares well with exact diagonalization results so the model presented here is able to provide a correct description of the strong-correlation effects.

Equations

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

$$<{\hat{S}}_z>=$$

$$<{\hat{S}}^2>=$$

$$D_{pq,rt}^{\sigma \sigma ,\sigma \sigma }$$

$$D_{pq,rt}^{\alpha \beta ,\alpha \beta }$$

$$E$$

$$-\frac{1}{2}\sum _{pq}\left(n_q^{\alpha }{n}_p^{\alpha }+n_q^{\beta }{n}_p^{\beta }\right){\mathcal{K}}_{pq}+\sum _{pqrt}{\tilde{\lambda }}_{pq,rt}⟨rt|pq⟩$$

$${\tilde{\lambda }}_{pq,rt}={\lambda }_{pq,rt}^{\alpha \alpha ,\alpha \alpha }+{\lambda }_{pq,rt}^{\alpha \beta ,\alpha \beta }+{\lambda }_{qp,tr}^{\alpha \beta ,\alpha \beta }+{\lambda }_{pq,rt}^{\beta \beta ,\beta \beta }$$

$$<{\hat{S}}^2>=\frac{N^{\alpha }+N^{\beta }}{2}+\frac{{\left(N^{\alpha }-N^{\beta }\right)}^2}{4}-\sum _p{n}_p^{\alpha }{n}_p^{\beta }-2\sum _{pq}{\lambda }_{pq,qp}^{\alpha \beta ,\alpha \beta }$$

$$\sum _p{m}_p=2S,2\sum _p{n}_p=N-2S$$

$$2\sum _q{\lambda }_{pq,qp}^{\alpha \beta ,\alpha \beta }=2\sum _q{\lambda }_{qp,pq}^{\alpha \beta ,\alpha \beta }=n_p-n_p^{\alpha }{n}_p^{\beta }$$

$${\lambda }_{pq,rt}^{\sigma \sigma ,\sigma \sigma }$$

$${\lambda }_{pq,rt}^{\alpha \beta ,\alpha \beta }$$

$$E$$

$$+\frac{1}{2}\sum _{pq}\left[\left(n_q^{\alpha }+n_q^{\beta }\right)\left(n_p^{\alpha }+n_p^{\beta }\right)-{\tilde{\Delta }}_{pq}\right]{\mathcal{J}}_{pq}$$

$$-\frac{1}{2}\sum _{pq}\left[\left(n_q^{\alpha }{n}_p^{\alpha }+n_q^{\beta }{n}_p^{\beta }\right)-\left({\Delta }_{pq}^{\alpha \alpha }+{\Delta }_{pq}^{\beta \beta }\right)\right]{\mathcal{K}}_{pq}$$

$${\tilde{\Delta }}_{pq}={\Delta }_{pq}^{\alpha \alpha }+{\Delta }_{pq}^{\alpha \beta }+{\Delta }_{qp}^{\alpha \beta }+{\Delta }_{pq}^{\beta \beta }$$

$${\Delta }_{pp}^{\alpha \beta }=n_p^{\alpha }{n}_p^{\beta },{\Pi }_{pp}=n_p$$

$$n_p^{\alpha }=n_p^{\beta }=n_p,m_p=0,p\in {\Omega }_P$$

$$\sum _{p\in {\Omega }_g}{n}_p=1,2\sum _{p\in {\Omega }_P}{n}_p=2\sum _{g=1}^{N_P/2}\sum _{p\in {\Omega }_g}{n}_p=N_P$$

$$\begin{array}[]{c}\Delta_{qp}=n_{p}^{2}\delta_{qp}+n_{q}n_{p}\left(1-\delta_{qp}\right)\delta_{q\Omega_{g}}\delta_{p\Omega_{g}}\\ \\ \Pi_{qp}=n_{p}\delta_{qp}+\Pi_{qp}^{g}\left(1-\delta_{qp}\right)\delta_{q\Omega_{g}}\delta_{p\Omega_{g}}\\ \\ \Pi_{qp}^{g}=\left\{\begin{array}[]{cc}-\sqrt{n_{q}n_{p}}\,,&p=g\textrm{ or }q=g\\ +\sqrt{n_{q}n_{p}}\,,&otherwise\end{array}\right.\\ \\ \delta_{q\Omega_{g}}=\begin{cases}1,&q\in\Omega_{g}\\ 0,&q\notin\Omega_{g}\end{cases};\quad g=1,2,\ldots,N_{P}/2\\ \\ \end{array}$$

$$\begin{array}[]{c}E_{P}^{pnof5}=\sum\limits_{g=1}^{N_{P}/2}E_{g}+\sum\limits_{f\neq g}^{N_{P}/2}E_{fg}\\ \\ E_{g}=\sum\limits_{p\in\Omega_{g}}n_{p}\left(2\mathcal{H}_{pp}+\mathcal{J}_{pp}\right)+\sum\limits_{q,p\in\Omega_{g},q\neq p}\Pi_{qp}^{g}\mathcal{L}_{pq}\\ \\ E_{fg}=\sum\limits_{p\in\Omega_{f}}\sum\limits_{q\in\Omega_{g}}\left[n_{q}n_{p}\left(2\mathcal{J}_{pq}-\mathcal{K}_{pq}\right)\right]\end{array}$$

$$E_P^{pnof7}=E_P^{pnof5}+\sum _{f\ne g}^{N_P/2}\sum _{p\in {\Omega }_f}\sum _{q\in {\Omega }_g}{\Pi }_{qp}^{\Phi }{\mathcal{L}}_{pq}$$

$$\sum _{p\in {\Omega }_U}{m}_p=2S=N_U$$

$$\begin{array}[]{c}E^{pnofi}=E_{P}^{pnofi}+E_{PU}+E_{U}\quad(i=5,7)\\ \\ E_{PU}=\sum\limits_{p\in\Omega_{P}}\sum\limits_{q\in\Omega_{U}}n_{p}\left(2\mathcal{J}_{pq}-\mathcal{K}_{pq}\right)\\ \\ E_{U}=\sum\limits_{p\in\Omega_{U}}\mathcal{H}_{pp}+\frac{1}{2}\sum\limits_{p,q\in\Omega_{U}}\left(\mathcal{J}_{pq}-\mathcal{K}_{pq}\right)\end{array}$$

$$\begin{array}[]{c}{\displaystyle\hat{H}=-t\sum_{\left\langle\mu,\upsilon\right\rangle,\sigma}\left(\hat{a}_{\mu,\sigma}^{\dagger}\hat{a}_{\upsilon,\sigma}+\hat{a}_{\mu,\sigma}\hat{a}_{\upsilon,\sigma}^{\dagger}\right)+U\sum_{\mu}\hat{n}_{\mu,\alpha}\hat{n}_{\mu,\beta}}\end{array}$$

Full text

The ground-state energy of any electronic system could be computed using the first-order reduced density matrix (1-RDM) according to Gibert’s theorem (Gilbert, 1975) along with the works of Donnelly and Parr (Donnelly and Parr, 1978), Levy (Levy, 1979) and Valone (Valone, 1980). In most calculations, the spectral decomposition of the 1-RDM is used to approximate the energy in terms of the naturals orbitals $\left\{\phi_{i}\right\}$ and their occupation numbers $\left\{n_{i}\right\}$. This formulation is called natural orbital functional (NOF) theory (NOFT). It is important to note that attaining the exact energy in density functional theory requires a greater effort than in NOFT, since the non-interacting part of the electronic Hamiltonian is actually a one-particle operator, so it has an explicit dependence on the 1-RDM. A detailed account of the state of the art of the NOFT can be found elsewhere (Piris, 2007; Piris and Ugalde, 2014; Pernal and Giesbertz, 2016).

So far, the exact reconstruction of the functional has not been achieved, therefore, we have to settle for making approximations. The most practical approach is to reconstruct solely the two-particle RDM (2-RDM) and employ the well-known energy expression in terms of it. Approximating the energy in that way implies that the functional still depends on the 2-RDM (Donnelly, 1979), hence the functional $N$-representability problem arises (Ludeña et al., 2013; Piris, 2018a). We have to impose the $N$-representability conditions (Mazziotti, 2012) to the 2-RDM reconstructed in terms of 1-RDM to guarantee that our approximate ground-state energy has physical sense. As far as we know, only the functionals developed by Piris and collaborators (PNOFi, i=1-7) (Piris, 2013a, 2014, 2017) rely on the reconstruction of the 2-RDM subject to necessary $N$-representability conditions.

Within NOFT, few attempts have been made to describe systems with spin polarization. The first calculations in this type of systems were probably reported in (Goedecker and Umrigar, 1998) by Goedecker and Umrigar, who used a NOF that can be traced back to Müller (M) (Muller, 1984) and Buijse and Baerends (BB) (Buijse, 1991) proposals. A formulation of this functional considering spin-independent natural orbitals (NOs) but spin-dependent occupation numbers (ONs) was applied by Lathiotakis et al. (Lathiotakis et al., 2005) to the first-row atoms. The major shortcoming of this approach is that it does not conserve, in general, the total spin. Rohr and Penal (Rohr and Pernal, 2011) showed that MBB functional has a non-positive fractional spin error, with the exception of one-electron systems. In fact, extensive violations of $N$-representability conditions for this functional, and others related to it, have been reported (Rodríguez-Mayorga et al., 2017).

The discovery (Klyachko, 2006; Altunbulak and Klyachko, 2008) of a systematic way to derive pure-state $N$-representability conditions for the 1-RDM, also known as generalized Pauli constraints (GPCs), has lately allowed to gain new insights on spin-uncompensated systems (Theophilou et al., 2015), as well as to open a new way to develop NOFs (Benavides-Riveros and Marques, 2018). Indeed, necessary pure $N$-representability conditions on the 2-RDM can be derived by applying the GPCs to a family of effective 1-RDMs (Mazziotti, 2016). It should be noted that the application of GPCs restricts the 1-RDM variational space leading to improvements in energy, but it does not improve the reconstruction of the approximate functional per se. A 1-RDM that represents a pure state does not guarantee that the reconstructed electron-electron potential energy will be pure-state $N$-representable, except if the reconstruction is the exact one. The functional $N$-representability problem continues to exist for pure 1-RDMs when we make approximations for the functional, and it is still related to the $N$-representability of the 2-RDM.

The first extension of the Piris reconstruction functional to high-spin multiplet states was reported in the past for PNOF1 (Leiva and Piris, 2007). Later, a necessary condition to ensure the conservation of the total spin was obtained for the two-particle cumulant matrix (Piris et al., 2009). PNOF3 showed exceptional performance for atoms and molecules without spin compensation (Piris et al., 2010a). Unfortunately, closer analysis of the dissociation curves for various diatomics revealed that PNOF3 overestimates the amount of electron correlation, when the non-dynamic electron correlation becomes important. It was demonstrated that this ill behavior was related to the violation of the $N$-representability conditions, in particular, to the violation of the G positivity conditions for the 2-RDM (Piris et al., 2010b).

In this work, we focus on two recent approximations that are based on pairs of electrons with opposite spins (Piris, 2018b). The first is the independent pair model PNOF5 (Piris et al., 2011, 2013), which remarkably turned out to be strictly pure $N$-representable (Pernal, 2013; Piris, 2013b). PNOF5 takes into account the important part of the dynamical electron correlation corresponding to the intra-pair interactions, and most of the non-dynamical effects. However, no inter pair electron correlation is accounted. To include the missing interactions, PNOF6 (Piris, 2014) and PNOF7 (Piris, 2017) were proposed. The latter was recently (Mitxelena et al., 2018a) improved by an adequate choice of sign factors for the inter-pair interactions. Since PNOF7 provides a robust description of non-dynamic correlation effects, we will limit ourselves to this interacting pair model.

Our goal is to extend PNOF5 and PNOF7 to describe systems with spin polarization. We follow the spin-restricted formulation proposed in (Leiva and Piris, 2007) which has the virtue of avoiding spin contamination effects. Our hypothesis is to take the occupation numbers of singly occupied orbitals that contribute to the spin of the system equal to one, whereas the rest of electron pairs with opposite spins are distributed in obital subspaces, with one pair per subspace. In the latter, the orbital occupancies are fractional, which account for the electron correlation. The resulting functionals are tested in energy calculations for model systems with strong non-dynamic electron correlation, namely, the one-dimensional Hubbard model with periodic boundary conditions and hydrogen rings up to 16 atoms, in order to compare with exact diagonalization values.

I Theory

In NOFT, the electronic energy is given in terms of the NOs and their ONs, namely,

[TABLE]

where $\mathcal{H}_{ii}$ denotes the diagonal elements of the one-particle part of the Hamiltonian involving the kinetic energy and the external potential operators, $<kl|ij>$ are the matrix elements of the two-particle interaction, and $D[n_{i},n_{j},n_{k},n_{l}]$ represents the reconstructed 2-RDM from the ONs. Restriction of the ONs to the range $0\leq n_{i}\leq 1$ represents a necessary and sufficient condition for ensemble $N$-representability of the 1-RDM (Coleman, 1963) under the normalization condition $\sum_{i}n_{i}=N$.

The $N$-electron Hamiltonian commonly used in electronic calculations does not contain any spin coordinates. Consequently, the eigenfunctions of the Hamiltonian are also eigenfunctions of both spin operators $\hat{S}_{z}$ and $\hat{S}^{2}$. For $\hat{S}_{z}$ eigenstates, only density matrix blocks that conserve the number of each spin type are non-vanishing. Specifically, the 1-RDM has two nonzero blocks $\alpha$ and $\beta$, whereas the 2-RDM has three independent nonzero blocks: $D^{\alpha\alpha}$, $D^{\alpha\beta}$, and $D^{\beta\beta}$.

Let us divide the spin-orbital set $\left\{\phi_{i}\left(\mathbf{x}\right)\right\}$ into two subsets: $\left\{\varphi_{p}^{\alpha}\left(\mathbf{r}\right)\alpha\left(\mathbf{s}\right)\right\}$ and $\left\{\varphi_{p}^{\beta}\left(\mathbf{r}\right)\beta\left(\mathbf{s}\right)\right\}$. In order to avoid spin contamination effects, the spin-restricted theory is employed, in which a single set of orbitals is used for $\alpha$ and $\beta$ spins: $\varphi_{p}^{\alpha}\left(\mathbf{r}\right)=\varphi_{p}^{\beta}\left(\mathbf{r}\right)=\varphi_{p}\left(\mathbf{r}\right)$. Accordingly, the expectation values of the spin operators read as (Piris, 2007)

[TABLE]

The matrix elements of the 2-RDM can conveniently be expressed in terms of the cumulant expansion:

[TABLE]

where $\mathbf{\lambda}$ is the cumulant matrix (Mazziotti, 1998). Using this expansion, the energy (1) reads

[TABLE]

where $\mathcal{J}_{pq}=\left\langle pq|pq\right\rangle$ and $\mathcal{K}_{pq}=\left\langle pq|qp\right\rangle$ are the usual Coulomb and exchange integrals, respectively. $\widetilde{\lambda}_{pq,rt}$ denotes the spinless cumulant matrix,

[TABLE]

From the sum rules of the 2-RDM blocks, it can be easily shown (Piris, 2007) that the expectation value of $\widehat{S}^{2}$ is

[TABLE]

For a given value $S$, there are $\left(2S+1\right)$ energy degenerate spin-multiplet states: $|SM_{S}>$, $M_{S}=-S,...,S$. Let us focus on the high-spin multiplet state ($M_{S}=S$) and assume that $n_{p}^{\alpha}=n_{p}+m_{p\,},\;n_{p}^{\beta}=n_{p}$ so that $n_{p}$ and $m_{p}$ fulfill the following constrains:

[TABLE]

Then, $<\widehat{S}^{2}>=S(S+1)$ implies (Piris et al., 2009) the following sum rule for the cumulant $\alpha\beta$-block:

[TABLE]

In general, the cumulant $\lambda$ has a dependence of four indexes, so it is expensive from the computational point of view to use such quantities. We shall use the reconstruction functional proposed in Ref. (Piris, 2006), in which the two-particle cumulant is explicitly constructed in terms of two-index matrices, $\Delta$ and $\Pi$. The latter are selected to satisfy necessary $N$-representability conditions and sum rules by the 2-RDM. A systematic application of $N$-representability conditions in the reconstruction of $\lambda$ has led to the PNOF series (Piris, 2013a, 2014, 2017). The latter has the following spin structure:

[TABLE]

which leads to the energy functional

[TABLE]

Here, $\Delta^{\sigma\sigma^{\prime}}$ are real symmetric matrices, whereas $\Pi$ is a spin-independent Hermitian matrix. These matrices satisfy several constraints imposed by the sum rules and $N$-representability conditions of the two-particle cumulant. $\widetilde{\Delta}_{pq}$ denotes the spinless $\Delta$ matrix,

[TABLE]

whereas the new integral $\mathcal{L}_{rp}=\left\langle rr|pp\right\rangle$ arises from the correlation between particles with opposite spins and is called the exchange and time-inversion integral (Piris, 1999).

By combining the sum rule (9) with the ansatz (11), one arrives at the following diagonal elements (Piris et al., 2009)

[TABLE]

that guarantee the conservation of the total spin.

In this work, we restrict ourselves to pairing-based approximations PNOF5 and PNOF7, that have been successfully implemented for singlets. Here, we extend both NOFs to situations in which a system could have one or more unpaired electrons, so we can have spin-polarized systems.

Consider $N_{U}$ unpaired electrons, and $N_{P}$ paired electrons, so that $N_{U}+N_{P}=N$. Similarly, divide the orbital space $\Omega$ into two subspaces: $\Omega=\Omega_{U}\oplus\Omega_{P}$. In $\Omega_{P}$, the spatial orbitals are double occupied ($n_{p}\neq 0$), whereas singly occupied orbitals ($n_{p}=0$) can only be found in $\Omega_{U}$. We assume further that all spins corresponding to the $N_{P}$ electrons are coupled as a singlet, so the occupancies for particles with $\alpha$ and $\beta$ spin are equal:

[TABLE]

According to the electron-pairing approach in NOFT (Piris, 2018b), the orbital space $\Omega_{P}$ is in turn divided into $N_{P}/2$ mutually disjoint subspaces $\Omega_{P}=\Omega_{1}\oplus\Omega_{2}\oplus...\oplus\Omega_{N_{P}/2}$. Each subspace $\Omega{}_{g}$ contains one orbital $g$ below the level $N_{P}/2$, and $N_{g}$ orbitals above it. In what follows, we consider $N_{g}$ equal to a fixed number that corresponds to the maximum value allowed by the basis set used. Taking into account the spin, the total occupancy for a given subspace $\Omega{}_{g}$ is 2, which is reflected in additional sum rules for the ONs, namely,

[TABLE]

The simplest way to satisfy the constraints imposed on the two-particle cumulant leads to PNOF5 (Piris et al., 2011, 2013):

[TABLE]

It is worth noting that $\Delta$ is a spin-independent matrix ($\Delta_{qp}=\Delta_{qp}^{\alpha\alpha}=\Delta_{qp}^{\alpha\beta}=\Delta_{qp}^{\beta\beta}$). In addition, $\mathrm{\Delta_{\mathit{qp}}}$ and $\mathrm{\Pi_{\mathit{qp}}}$ are zero between orbitals belonging to different subspaces, so the 2-RDM reconstruction of PNOF5 corresponds to an independent pair model. The resulting energy for the $N_{P}$ electrons is

[TABLE]

where $\mathcal{J}_{pp}=<pp|pp>$ is the Coulomb interaction between two electrons with opposite spins at the spatial orbital $p$.

To go beyond the independent-pair approximation, $\mathrm{\Delta_{\mathit{qp}}=0}$ was kept, while nonzero $\Pi$-elements were considered between orbitals belonging to different subspaces. Accordingly, new elements $\Pi_{qp}^{\Phi}=-\Phi_{q}\Phi_{p}$ with $\Phi_{q}=\sqrt{n_{q}(1-n_{q})}$ were included for orbitals $q$ and $p$ belonging to different pairs (subspaces) in PNOF7 (Piris, 2017; Mitxelena et al., 2018a). Thus, the energy of the $N_{P}$ electrons becomes

[TABLE]

On the other hand, taking into account Eq. (16), it follows from (8) that

[TABLE]

To attain a high-spin multiplet state, let us take $m_{p}=1$ and dim$\left\{\Omega_{U}\right\}=N_{U}$. This assumption is trivial for a doublet, but it is more restrictive for higher multiplets because an underestimation of the energy can occur. Indeed, the correlation $\Delta$ and $\Pi$ matrices, which determine the PNOF two-electron cumulant, are null at the boundary values of the ONs (Piris, 2013a). Therefore, the single-occupied orbitals ($n_{p}^{\alpha}=1,\;n_{p}^{\beta}=0,\;p\in\Omega_{U}$) are not allowed, by construction, to contribute to the correlation energy. This assumption is the most critical for spin-polarized systems.

After some algebra, the energy (12) for a spin-uncompensated system can be written as

[TABLE]

The solution is established by optimizing the energy functional (21) 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 1-RDM. In general, orbitals below $N_{P}/2$ are almost fully occupied, and those that are above the single-occupied orbitals are almost empty, however, in a system with a strong non-dynamic electron correlation, several NOs can have ONs equal or close to 1/2.

The self-consistent procedure proposed in (Piris and Ugalde, 2009) yields the NOs by the iterative diagonalization of a Hermitian matrix F. The off-diagonal elements of F are determined explicitly by the hermiticity of the Lagrange multipliers. The first-order perturbation theory applying to each cycle of the diagonalization process provides an aufbau principle for determining the diagonal elements of F.

In the following sections, we analyze the 1D Hubbard model with different number of sites, and hydrogen chains with different ring sizes, for different spin values, in order to test the NOFs given by Eq. (21) in strong non-dynamic correlation regimes.

II Hubbard model

The Hubbard model is an ideal candidate for the study of electron correlation in solid state physics and quantum chemistry. In second quantization notation, the one-dimensional (1D) Hubbard Hamiltonian takes the form (Baeriswyl et al., 1995)

[TABLE]

where Greek indices $\mu$ and $\upsilon$ denote sites of the model, $\left\langle\mu,\upsilon\right\rangle$ indicates only near-neighbors interactions, $t$ is the hopping parameter, $U$ is the on-site inter-electron repulsion parameter, and $\hat{n}_{\mu,\sigma}=\hat{a}_{\mu,\sigma}^{\dagger}\hat{a}_{\mu,\sigma}$ where $\hat{a}_{\mu,\sigma}^{\dagger}\left(\hat{a}_{\mu,\sigma}\right)$ corresponds to the fermionic creation(annihilation) operator. The first term models the kinetic energy of the electrons hopping between atoms, whereas the second term is the potential energy that emerges from the on-site interaction $U$.

It is well-known that the Hartree-Fock (HF) approximation retrieves the exact solution for the 1D Hubbard model at half-filling if $U=0$, which in this limit corresponds to the tigh-binding model (Slater and Koster, 1954). Conversely, in the $U/t\rightarrow\infty$ limit the model becomes equivalent to the spin-1/2 Heisenberg model (Timm, 2015). Recently, the performance of commonly used NOFs for singlet states in the 1D Hubbard model has been assessed (Mitxelena et al., 2017; Mitxelena et al., 2018b, a), showing that PNOF7 is in good agreement with exact diagonalization (ED) results for the Hubbard model at half-filling on even number of sites. Since the purpose of this work is to introduce an approach for studying spin-uncompensated systems with the NOFs developed in our group, in this work we analyze the performance of spin-uncompensated PNOFs by using the Hubbard model. Actually, our main goal is to study many multiplicities in order to prove the advantages and shortcomings of our new approximation.

Many spin multiplicities are of interest in electronic structure theory. Thus, a complete test should include multiplicities ranging from singlets to octets, the latter observed, for example, in Gadolinium. In the following, we study the performance of spin-uncompensated PNOF5 and PNOF7 by using the Hubbard model, and employing the $U/t$ ratio to cover all correlation regimes. The maximum number of sites considered is 14 in order to compare with ED calculations. PNOF calculations have been carried out using DoNOF code developed by M. Piris and coworkers based on the iterative diagonalization method (Piris and Ugalde, 2009), whereas ED results have been computed using a modified version of the code developed by Knowles and Handy (Knowles and Handy, 1984, 1989).

In Figs. 1-3, we show the energy differences obtained with PNOF5 and PNOF7 with respect to the ED values varying $U/t$. In general, we can observe that the error increases to a maximum located at $U/t\thickapprox 5$, and then seems to decrease to a non-zero minimum at the $U/t\rightarrow\infty$ limit. In comparison with the performance obtained for spin-compensated systems (Mitxelena et al., 2017; Mitxelena et al., 2018a, b), the errors shown for spin-polarized systems are slightly larger for all correlation regimes. Recall that PNOF5 takes into account the whole intra-pair electron correlation, while PNOF7 also includes non-dynamic correlation between electron pairs, so it is expected that errors with respect to ED are greater for PNOF5 than for PNOF7.

At the $U/t\rightarrow\infty$ limit, also known as the strong correlation limit, electrons try to keep away one from each other by half-filling the sites, which corresponds to the Mott-Hubbard regime. In this case, there is no kinetic energy, all ONs are exactly equal to $1/2$, which corresponds to a pure non-dynamic correlation regime. Therefore, the fact that the energy difference in the $U/t\gg 1$ region is different from zero is related to the lack of non-dynamic electron correlation. Even using PNOF7, our model only takes into account part of this correlation between electrons with opposite spins in the orbital subspace $\Omega_{P}$, that is, the non-dynamic correlation of the $N_{P}$ paired electrons. As can be seen from Eq. (21), the $E_{PU}$ term is of HF type interactions hence our approximate NOFs do not include correlation between electrons with opposite spins involving the subspace $\Omega_{U}$. Singly occupied orbitals do not contribute to the electron correlation (vide supra).

On the other hand, both functionals PNOF5 and PNOF7 underestimate dynamic inter-pair correlation effects, which turns out to be important at low and intermediate $U/t$ values, namely $U/t\approx 5$. Recently (Piris, 2017, 2018c), one of us (M.P.) has proposed a new method to recover the missing dynamic correlation in spin-compensated systems by means of perturbation theory. Corrections of this type are not within the scope of this work and will be incorporated for spin-polarized systems in the future. In what follows, we focus mainly on the region where the effects of non-dynamic correlation predominate, that is, for large $U/t$.

In Fig. 1, we plot the results for systems with odd number of sites (3, 5, 11 and 13) and spin $S=1/2$ ($N_{U}=1$). Note that there are no differences between PNOF5 and PNOF7 in the 3 sites system with $S=1/2$ because there is only one electron pair with opposite spins considered ($N_{P}=2$). We also note that errors with respect to ED results increase by rising sites from 5 ($N_{P}=4$) to 13 ($N_{P}=12$). Indeed, the electron correlation is increasingly ignored by accounting for more $N_{P}$ pairs for one unpaired electron. The introduction in Eq. (19) of the term to account for non-dynamic inter-pair correlation makes PNOF7 more robust with the increase in the number of electron pairs, which is consistent with the results obtained for singlet states (Mitxelena et al., 2018a).

In Fig. 2, the triplet states ($N_{U}=2$) are considered for systems with an even number of sites. It is worth mentioning that for 4 sites, the PNOF5 and PNOF7 values coincide due to having only one pair with opposite spins ($N_{P}=2$). We observe that the trend is similar to that obtained for $S=1/2$, however, the errors are higher in the $S=1$ case (notice that the scale of the $\Delta$E axis is different between figures 1 and 2). There are two unpaired electrons, whose correlation effects are neglected, so that errors increase for all systems with $S=1$ with respect to the corresponding systems with $S=1/2$.

Figures 1 and 2 show how our approximation behaves for a given $N_{U}$ by varying the value of $N_{P}$. In both figures, we observe that the curves approximate each other by increasing the number of sites for a given value of $U/t$, especially in the region of high non-dynamic correlation. Consequently, the error with respect to the ED values tends to stabilize.

We are also interested in the performance of our approach for a fixed $N_{P}$. Consequently, we decided to analyze which effects cause the increase of the system size and its spin multiplicity at the same time. In Fig. 3, we start with the 11 sites Hubbard model ($N_{U}=7,\,S=7/2$) with only two electron pairs ($N_{P}=4$). Then, we keep fixed $N_{P}=4$ and add an additional site and an extra unpaired electron consecutively ($N_{U}=8,9,10,\,S=4,9/2,5$). Neither PNOF5 nor PNOF7 yield significant differences for $U/t<5$, however, for larger $U/t$ values, where non-dynamic correlation regimes prevail, PNOF7 performs better than PNOF5 as expected. This result indicates that the error of our functionals with respect to the ED results seems to stabilize for the variation of the spin multiplicity together with the size for a given number $N_{P}$ of electron pairs as well.

In Fig. 4, we have collected the results obtained with PNOF7 for different multiplicities and a given number of electrons (14 sites Hubbard model). The best results are obtained for minimum and maximum multiplicities. For $S=0$, a great amount of electron correlation is recovered by PNOF7 for any correlation regime. For large spin values, the HF energy components dominate up to the full polarization (S = 7) in which there is no correlation energy. For the latter, the exact result is reduced to the HF energy, which is also recovered with our approach. For the intermediate spin values, the interactions related to the $E_{PU}$ energy component in Eq. (21) should be improved to reduce the relatively large errors in Fig. 4.

III Hydrogen rings

The lack of long-range inter-electronic interactions may be the most important limitation of the Hubbard model. In this section, we focus on rings composed of hydrogen atoms, which model the non-dynamic electronic correlation in the presence of long-range interaction effects, to examine whether the results obtained in the previous section are still valid. Here, $<kl|ij>$ in Eq. (1) are the matrix elements of the bare Coulomb interaction. These systems have previously been used to perform benchmarking of many-body approximations (Boguslawski et al., 2014; Mitxelena et al., 2017). Theoretical investigations on thermodynamic stabilities (Yamaguchi et al., 1983; Wright and DiLabio, 1992) and aromaticity (Jiao et al., 1996) of hydrogen rings have also been carried out.

In Fig. 5, we plot total energy differences with respect to the ED results obtained by using PNOF5, PNOF7 and unrestricted HF (UHF) for $S=1/2(a),1(b),3/2(c)$ and $2(d)$. The maximum number of hydrogens considered was 16 in order to compare with ED calculations. Total energies corresponding to ED and UHF have been carried out by means of the Psi4 suite of programs (Parrish et al., 2017). All calculations were done with the minimal STO-3G basis set (Hehre et al., 1969). The inter-nuclear distance between hydrogen atoms was fixed to $R_{H-H}=2.0\>\textrm{\r{A}}$ in order to have mostly non-dynamic electron correlation (Sinitskiy et al., 2010). Interestingly, PNOF5 and UHF behave similarly for any spin multiplicity, however, UHF retrieves electron correlation at the expense of breaking the total spin symmetry, while PNOF5 affords the correct value of $<\hat{S}^{2}>$.

As already seen for the Hubbard model, the error obtained when using PNOF7 is significantly reduced with respect to PNOF5. Taking into account the curves shown in Fig. 5, we can observe that the difference $E^{pnof7}-E^{ED}$ seems to tend earlier to a constant value than the energies of PNOF5 and UHF, which produces larger errors as the number of hydrogens increases, reaching the asymptote presumably for much larger rings. This represents an important step forward in the description of spin-polarized systems with strong correlation effects, since PNOF7 shows a significant improvement over the widely used mean-field UHF method.

Finally, in Fig. 6 we report the total energies ($E(H_{n})$) divided by the number ($n$) of hydrogens as a function of the ring size using PNOF5, PNOF7, UHF and ED. Comparing $E(H_{n})/n$ obtained with the different approximations with respect to the ED values, we observe that PNOF7 is in good agreement with the latter, although the error increases slightly as S increases.

IV Closing Remarks

The extension of PNOF5 and PNOF7 has been achieved for systems with spin polarization. The obtained natural orbital functionals provide an adequate value of $<\hat{S}^{2}>$ and $<\hat{S}_{z}>$, so that the spin contamination effects are not present. As in the case of singlet states, N-representability conditions on the reconstructed two-particle reduced density matrix are satisfied for high-spin multiplet state.

In the approach proposed here, we have considered a subspace of fully-singly occupied orbitals with the same spin, while the rest of the electron pairs with opposite spins are distributed in the remaining orbital subspace. In this way, the electron pairs form a singlet, and the unpaired electrons are responsible for the spin of the system. This simplification works well for a doublet, but it is more restrictive for higher multiplets because the single-occupied orbitals are not allowed to contribute to the correlation energy.

Despite its simplicity, the model presented here is able to provide a qualitatively correct description of the strong-correlation effects that appear in the one-dimensional Hubbard model and the hydrogen rings. A behavior similar to that obtained for singlets, especially for doublets and triplets, was achieved for these systems, so it has been shown that the performance of PNOF5 and PNOF7 does not deteriorate when they are extended to spin-polarized systems.

PNOF5 fully takes into account the intra-pair electron correlation, while PNOF7 also includes non-dynamic correlation between electron pairs, so the errors with respect to exact diagonalization are greater for PNOF5 than for PNOF7. It was shown that errors seem to stabilize as the system size increases regardless of the value of the spin. In addition, PNOF7 showed a significant improvement over the widely used unrestricted Hartree-Fock method for spin-polarized systems. Therefore, PNOF7 is an ideal candidate to be employed in the study of spin-uncompensated systems composed of a large number of atoms.

The best results were obtained for minimum and maximum multiplicities. For the intermediate spin values, the interactions related to the paired-unpaired energy component should be improved to reduce the relatively errors. In this vein, an explicit form of the two-electron cumulant of the reconstructed two-particle reduced matrix capable to describe fractional occupancies of the unpaired electrons is still missing. A work in this direction is underway.

V Acknowledgments

Financial support comes from Ministerio de Economía y Competitividad (Ref. CTQ2015-67608-P). The authors thank for technical and human support provided by IZO-SGI SGIker of UPV/EHU and European funding (ERDF and ESF). R. Q-M. is grateful to Conacyt for the grant “Fronteras” (Project 867) and the PhD grant 587950. I. M. is grateful to Vice-Rectory for research of the UPV/EHU for the PhD. grant (PIF//15/043).

?refname?

• Gilbert (1975)

T. L. Gilbert, Phys. Rev. B 12, 2111 (1975).

• Donnelly and Parr (1978)

R. A. Donnelly and R. G. Parr, J. Chem. Phys. 69, 4431 (1978).

• Levy (1979)

M. Levy, Proc. Natl. Acad. Sci. USA 76, 6062 (1979).

• Valone (1980)

S. M. Valone, J. Chem. Phys. 73, 1344 (1980).

• Piris (2007)

M. Piris, in Reduced-Density-Matrix Mechanics: with applications to many-electron atoms and molecules, edited by D. A. Mazziotti (John Wiley and Sons, Hoboken, New Jersey, USA, 2007), chap. 14, pp. 387–427.

• Piris and Ugalde (2014)

M. Piris and J. M. Ugalde, Int. J. Quantum Chem. 114, 1169 (2014).

• Pernal and Giesbertz (2016)

K. Pernal and K. J. H. Giesbertz, Top Curr Chem 368, 125 (2016).

• Donnelly (1979)

R. A. Donnelly, J. Chem. Phys. 71, 2874 (1979).

• Ludeña et al. (2013)

E. V. Ludeña, F. J. Torres, and C. Costa, J. Mod. Phys. 04, 391 (2013).

• Piris (2018a)

M. Piris, in Many-body approaches at different scales: a tribute to N. H. March on the occasion of his 90th birthday, edited by G. G. N. Angilella and C. Amovilli (Springer, New York, 2018a), chap. 22, pp. 283–300.

• Mazziotti (2012)

D. A. Mazziotti, Phys. Rev. Lett. 108, 263002 (2012).

• Piris (2013a)

M. Piris, Int. J. Quantum Chem. 113, 620 (2013a).

• Piris (2014)

M. Piris, J. Chem. Phys. 141, 044107 (2014).

• Piris (2017)

M. Piris, Phys. Rev. Lett. 119, 063002 (2017).

• Goedecker and Umrigar (1998)

S. Goedecker and C. Umrigar, Phys. Rev. Lett. 81, 866 (1998).

• Muller (1984)

A. M. K. Muller, Phys. Lett. 105A, 446 (1984).

• Buijse (1991)

M. A. Buijse, Thesis (Vrije Universiteit, Amsterdam, The Netherlands, 1991).

• Lathiotakis et al. (2005)

N. Lathiotakis, N. Helbig, and E. Gross, Phys. Rev. A 72, 030501(R) (2005).

• Rohr and Pernal (2011)

D. R. Rohr and K. Pernal, J. Chem. Phys. 135, 074104 (2011).

• Rodríguez-Mayorga et al. (2017)

M. Rodríguez-Mayorga, E. Ramos-Cordoba, M. Via-Nadal, M. Piris, and E. Matito, Phys. Chem. Chem. Phys 19, 24029 (2017).

• Klyachko (2006)

A. a. Klyachko, J. Phys. Conf. Ser. 36, 72 (2006).

• Altunbulak and Klyachko (2008)

M. Altunbulak and A. Klyachko, Commun. Math. Phys. 282, 287 (2008).

• Theophilou et al. (2015)

I. Theophilou, N. N. Lathiotakis, M. a. L. Marques, and N. Helbig, J. Chem. Phys. 142, 154108 (2015).

• Benavides-Riveros and Marques (2018)

C. L. Benavides-Riveros and M. A. L. Marques, Eur. Phys. J. B 91, 133 (2018).

• Mazziotti (2016)

D. A. Mazziotti, Phys. Rev. A 94, 032516 (2016).

• Leiva and Piris (2007)

P. Leiva and M. Piris, Int. J. Quantum Chem. 107, 1 (2007).

• Piris et al. (2009)

M. Piris, J. M. Matxain, X. Lopez, and J. M. Ugalde, J. Chem. Phys. 131, 021102 (2009).

• Piris et al. (2010a)

M. Piris, J. M. Matxain, X. Lopez, and J. M. Ugalde, J. Chem. Phys. 132, 031103 (2010a).

• Piris et al. (2010b)

M. Piris, J. M. Matxain, X. Lopez, and J. M. Ugalde, J. Chem. Phys. 133, 111101 (2010b).

• Piris (2018b)

M. Piris, in Theoretical and Quantum Chemistry at the Dawn of the 21st Century, edited by T. Chakraborty and R. Carbó-Dorca (Apple Academic Press, New Jersey, 2018b), chap. 22, pp. 593–620.

• Piris et al. (2011)

M. Piris, X. Lopez, F. Ruipérez, J. M. Matxain, and J. M. Ugalde, J. Chem. Phys. 134, 164102 (2011).

• Piris et al. (2013)

M. Piris, J. M. Matxain, and X. Lopez, J. Chem. Phys. 139, 234109 (2013).

• Pernal (2013)

K. Pernal, Comp. Theor. Chem. 1003, 127 (2013).

• Piris (2013b)

M. Piris, J. Chem. Phys. 139, 064111 (2013b).

• Mitxelena et al. (2018a)

I. Mitxelena, M. Rodríguez-Mayorga, and M. Piris, Eur. Phys. J. B 91, 109 (2018a).

• Coleman (1963)

A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963).

• Mazziotti (1998)

D. A. Mazziotti, Chem. Phys. Lett. 289, 419 (1998).

• Piris (2006)

M. Piris, Int. J. Quantum Chem. 106, 1093 (2006).

• Piris (1999)

M. Piris, J. Math. Chem. 25, 47 (1999).

• Piris and Ugalde (2009)

M. Piris and J. M. Ugalde, J. Comput. Chem. 30, 2078 (2009).

• Baeriswyl et al. (1995)

D. Baeriswyl, D. K. Campbell, J. M. P. Carmelo, F. Guinea, and E. Louis, eds., The Hubbard Model: Its Physics and Mathematical Physics (Nato Science Series B), vol. 343 (Springer US, 1995).

• Slater and Koster (1954)

J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).

• Timm (2015)

C. Timm, Theory of Magnetism: International Max Planck Research School for Dynamical Processes in Atoms, Molecules and Solids (Technische Univ. Dresden, Institute for Theoretical Physics, Dresden, Germany, 2015).

• Mitxelena et al. (2017)

I. Mitxelena, M. Piris, and M. Rodríguez-Mayorga, J. Phys. Condens. Matter 29, 425602 (2017).

• Mitxelena et al. (2018b)

I. Mitxelena, M. Piris, and M. Rodríguez-Mayorga, J. Phys. Condens. Matter 30, 089501 (2018b).

• Knowles and Handy (1984)

P. Knowles and N. Handy, Chem. Phys. Lett. 111, 315 (1984), ISSN 0009-2614.

• Knowles and Handy (1989)

P. J. Knowles and N. C. Handy, J. Chem. Phys. 91, 2396 (1989).

• Piris (2018c)

M. Piris, Phys. Rev. A 98, 022504 (2018c).

• Boguslawski et al. (2014)

K. Boguslawski, P. Tecmer, P. W. Ayers, P. Bultinck, S. De Baerdemacker, and D. Van Neck, Phys. Rev. B 89, 201106 (2014).

• Yamaguchi et al. (1983)

Y. Yamaguchi, J. F. Gaw, and H. F. Schaefer, The Journal of Chemical Physics 78, 4074 (1983).

• Wright and DiLabio (1992)

J. S. Wright and G. A. DiLabio, The Journal of Physical Chemistry 96, 10793 (1992).

• Jiao et al. (1996)

H. Jiao, P. v. R. Schleyer, and M. N. Glukhovtsev, The Journal of Physical Chemistry 100, 12299 (1996).

• Parrish et al. (2017)

R. M. Parrish, L. A. Burns, D. G. Smith, A. C. Simmonett, A. E. DePrince, E. G. Hohenstein, U. Bozkaya, A. Y. Sokolov, R. Di Remigio, R. M. Richard, et al., J. Chem. Theory Comput. 13, 3185 (2017).

• Hehre et al. (1969)

W. J. Hehre, R. F. Stewart, and J. A. Pople, J. Chem. Phys. 51, 2657 (1969).

• Sinitskiy et al. (2010)

A. V. Sinitskiy, L. Greenman, and D. A. Mazziotti, J. Chem. Phys. 133, 014104 (2010).