The direct variational optimization of the ground-state two-electron reduced-density matrix (2-RDM) is typically performed under ensemble N-representability conditions. Accordingly, variationally obtained 2-RDMs for degenerate ground states may not represent a pure state. When considering only ground-state energetics, the ensemble nature of the 2-RDM is of little consequence. However, the use of ensemble densities within an extended random phase approximation (ERPA) yields astonishingly poor estimates of excitation energies, even for simple atomic systems [H. van Aggelen et al., Comput. Theor. Chem. 1003, 50–54 (2013)]. Here, we outline an approach for the direct variational optimization of ground-state 2-RDMs that satisfy pure-state N-representability known as generalized Pauli constraints. Within the ERPA, 2-RDMs that satisfy both ensemble conditions and the generalized Pauli constraints yield much more reliable estimates of excitation energies than those that satisfy only ensemble conditions.

The ground-state energy of an N-electron system is a known linear functional of the two-electron reduced-density matrix (2-RDM).1 Because the 2-RDM is a more compact mathematical object than the N-electron wave function, it seems natural to seek a procedure to directly determine the elements of the 2-RDM without the use of the wave function. Unfortunately, such a procedure is complicated by the fact that the 2-RDM must be constrained to ensure that it derives either from a single antisymmetrized N-electron wave function or an ensemble of wave functions. A reduced-density matrix (RDM) that satisfies this property is said to be N-representable.2 For example, an ensembleN-representable 2-RDM (2D) can be obtained from the partial integration of an ensemble N-electron density matrix as

(1)
$2 D ( x 1 , x 2 | x 1 ′ , x 2 ′ ) = N 2 ∫ d x 3 … d x N$
$∑ k w k Ψ k ( x 1 , x 2 , x 3 , … , x N ) Ψ k ∗ ( x 1 ′ , x 2 ′ , x 3 , … , x N ) .$
Here, Ψk represents a normalized and antisymmetrized N-electron wave function, the symbol x1 (x2, etc.) represents the spin and spatial coordinates of electron 1 (electron 2, etc.), and the non-negative weights, wk, sum to one. If only one coefficient in this expansion is non-zero, then the 2-RDM is said to be pure-stateN-representable.

Trivial conditions for the N-representability of a p-body RDM (where p < N) include its Hermiticity and antisymmetry with respect to the exchange of particles or holes. A hierarchy of p-positivity conditions also requires that the various representations of the p-RDM be positive semidefinite and interrelated according to the anticommutation relations of fermionic creation and annihilation operators.3 For example, at the one-particle level, the one-electron RDM (1-RDM, 1D) and the one-hole RDM (1Q) must satisfy

$1 D ⪰ 0 , with 1 D j i = 〈 Ψ | a ˆ i † a ˆ j | Ψ 〉 ,$
(2)
$1 Q ⪰ 0 , with 1 Q i j = 〈 Ψ | a ˆ j a ˆ i † | Ψ 〉 ,$
(3)

and

$1 D j i + 1 Q i j = δ i j .$
(4)

Here, it has become useful to introduce the notation of second quantization; the symbols $a ˆ †$ and $a ˆ$ represent creation and annihilation operators, respectively, and the indices i and j are spin orbital labels. The notation 1D⪰0 indicates that all of the eigenvalues of the 1-RDM are non-negative. Taken together, Eqs. (2)–(4) imply that the eigenvalues of an N-representable 1-RDM must lie between zero and one. Indeed, these bounds, when combined with the trace condition, Tr(1D) = N, constitute Coleman’s necessary and sufficient ensemble N-representability conditions for the 1-RDM.2 Similarly, at the two-particle level, the 2-RDM, the two-hole RDM, and the electron-hole RDM must be positive semidefinite and interrelated. These necessary yet insufficient conditions on the N-representability of the 2-RDM are the “PQG” conditions of Garrod and Percus.4

The variational determination of the ground-state 2-RDM under such ensemble N-representability conditions can routinely be achieved via semidefinite programming (SDP) techniques.3,5–15 On the other hand, the direct determination of excited-state densities and energies is much more challenging. Without the imposition of additional constraints that differentiate the ground- and excited-state 2-RDMs, variational approaches can only yield the 2-RDM for the lowest-energy state of a given spin symmetry. Hence, several methods have been proposed to extract excited-state information from ground-state 2-RDMs, including the Hermitian operator method16–20 and the extended random-phase approximation (ERPA);21,22 both of these approaches seem to work well with variationally determined 2-RDMs for simple closed-shell systems. However, van Aggelen et al.22 demonstrated a striking failure of the ERPA for open-shell second-row atoms with degenerate ground states. In such cases, 2-RDMs optimized under ensemble N-representability conditions do not represent a pure state, and the ensemble nature of the 2-RDMs compromises the quality of the corresponding ERPA spectra. This result is unfortunate, as it suggests that the utility of the ERPA in the context of variational 2-RDM (v2RDM) methods is limited to systems with non-degenerate ground states.

Necessary and sufficient pure-state N-representability conditions for the 1-RDM, also called generalized Pauli constraints (GPC), are much more complicated than the ensemble conditions, and, until a short while ago, such conditions were only known for special cases.23 Recently, Altunbulak and Klyachko devised an algorithm to systematically identify all GPC for systems composed of an arbitrary numbers of electrons distributed among a finite number of natural spin orbitals.24 Now that expressions for the GPC are known, several recent studies have explored where pure N-representable solutions to the many-electron problem lie relative to the boundary of the polytope formed by the constraints.25–28 There is some debate as to whether ground states are exactly pinned or only quasi-pinned to the boundary of the polytope, yet excited states can clearly be buried deep within it.27 Beyond these geometric studies, Theophilou et al.29 have explored the role of the GPC within reduced-density-matrix functional theory (RDMFT).30 For the ground states of several three electron systems, multiple GPC were consistently violated, regardless of the choice of RDMFT functional. Further, the quality the RDMFT energy improved when the computations enforced the GPC. We do not necessarily expect the imposition of GPC to similarly impact the quality of the ground-state energy within a v2RDM optimization, but the failures of the ERPA demonstrated in Ref. 22 suggest that the constraints could play an important role in that context.

In this work, we develop a procedure to impose the generalized Pauli constraints within a ground-state v2RDM optimization. For many-electron systems with degenerate ground states, we demonstrate that the ERPA excitation energies from variationally obtained 2-RDMs that satisfy these pure-state N-representability conditions are in good agreement with excitation energies obtained from complete active space self-consistent field (CASSCF) computations. On the other hand, the ERPA spectra for RDMs that satisfy only ensemble N-representability conditions are unreliable. In this case, agreement between the ERPA and CASSCF is sometimes poor, and the ERPA spectra vary wildly, depending on the initial v2RDM optimization conditions.

The ground-state electronic energy is an exact functional of the one- and two-electron reduced-density matrices,

$E = 1 2 ∑ i j k l 2 D k l i j ( i k | j l ) + ∑ i j 1 D j i h i j .$
(5)

Here, the symbol (ik|jl) represents a two-electron repulsion integral in Mulliken notation, hij represents the sum of kinetic and electron-nucleus potential energy integrals, and the elements of the 2-RDM are defined in second quantized notation as

$2 D k l i j = 〈 Ψ | a ˆ i † a ˆ j † a ˆ l a ˆ k | Ψ 〉 .$
(6)

The indices i, j, k, and l run over all spin orbitals.

The minimization of the electronic energy given by Eq. (5) subject to known N-representability conditions constitutes a semidefinite optimization problem. Specifically, it is the primal formulation of the problem, expressed in the language of SDP as

$minimize E primal = c T ⋅ x such that Ax = b and M ( x ) ⪰ 0 .$
(7)

Here, x represents the primal solution vector, Eprimal represents the energy of the primal solution, and the matrix A and vector b offer a compact representation of the N-representability constraints (see below). When enforcing the two-particle (PQG) N-representability conditions, the mapping M(x) maps the primal solution vector onto the RDMs as

$M ( x ) = 1 D 0 0 0 0 0 1 Q 0 0 0 0 0 2 D 0 0 0 0 0 2 Q 0 0 0 0 0 2 G .$
(8)

Equation (8) introduces two additional 2-body RDMs; the symbols 2Q and 2G represent the two-hole RDM and the electron-hole RDM, respectively, the elements of which are

$2 Q k l i j = 〈 Ψ | a ˆ i a ˆ j a ˆ l † a ˆ k † | Ψ 〉$
(9)

and

$2 G k l i j = 〈 Ψ | a ˆ i † a ˆ j a ˆ l † a ˆ k | Ψ 〉 .$
(10)

The vector c similarly maps onto the one- and two-electron integrals such that Eprimal and the energy given by Eq. (5) are equivalent.

The notation M(x)⪰0 indicates that the eigenvalues of each RDM are non-negative. In addition to the positivity of the RDMs, necessary N-representability conditions require that the blocks of M(x) correctly map to one another according to the anticommutation relations of fermionic creation and annihilation operators. These linear equality constraints are represented within the SDP problem by the action of the constraint matrix, A, on the primal solution vector: Ax = b. For example, if the dimension of the one-electron basis is k, then the relationship between 1D and 1Q given in Eq. (4) can be mapped onto the action of k2 rows of the matrix A on x; the right-hand side of Eq. (4) can similarly be mapped onto k2 elements of the constraint vector b. Additional equality constraints, including spin, trace, and contraction conditions,31,32 can also be represented in this way.

Alternatively, the semidefinite optimization problem can be cast in terms of its dual formulation,

$maximize E dual = b T ⋅ y such that z = c − A T y and M ( z ) ⪰ 0 ,$
(11)

where y and z are the dual solution vectors, Edual represents the energy of the dual solution, and M(z) must also be non-negative. The symbols bT and AT represent the transpose of the constraint vector and matrix described above.

To solve the SDP problem, we use a boundary-point optimization, similar to those outlined in Refs. 15 and 33–35, that consists of the following iterative two-step procedure:

1. Solve AATy = A(cz) + μ(bAx) for y.

2. Update M(x) = μ−1Ω+ and M(z) = − Ω, where Ω = M(μx + ATyc).

Here, μ is a parameter that drives the optimization toward either the primal or dual solution; μ can be updated in the course of the optimization according to the prescription in Ref. 35. Step 1 is solved by conjugate gradient methods, and the positive/negative components of Ω are separated by diagonalization. Steps 1 and 2 are repeated until the primal error ‖Axb‖, the dual error ‖ATyc + z‖, and the difference between the primal and dual energies are sufficiently small.

The ensemble N-representability conditions require that the natural orbital occupation numbers (the eigenvalues of the 1-RDM, λi) lie between zero and one,

$0 ≤ λ i ≤ 1 .$
(12)

These bounds are enforced by the procedure outlined above. The pure-state N-representability conditions place additional equality/inequality constraints on the natural orbital occupation numbers. For example, consider a system of three-electrons distributed among six spin orbitals. The Hilbert space in which the electrons reside is represented by $H 6$, and the family of antisymmetrized three-electron wave functions in this Hilbert space is denoted by $∧ 3 H 6$. Table I provides the GPC satisfied by a pure-state N-representable 1-RDM for the family $∧ 3 H 6$. Here, the natural orbitals are ordered such that λiλi+1. Note that constraint 1 in Table I is usually expressed as an inequality constraint (λ4λ5λ6 ≤ 0), but we have recast it as an equality constraint via the introduction of an additional variable quantity, g1, which must be non-negative. For a system composed of a larger numbers of electrons or described by a Hilbert space of a larger rank, the GPC include many more inequality constraints. Hence, the general case requires that we introduce many gp, where p indicates a particular constraint. The constraints can be written such that each gp must be non-negative, and, in the course of the optimization, the positivity of gp can be enforced along with that of the RDMs. That is, we can define the modified mapping

$M ( x ) = 1 D 0 0 0 0 0 0 1 Q 0 0 0 0 0 0 2 D 0 0 0 0 0 0 2 Q 0 0 0 0 0 0 2 G 0 0 0 0 0 0 g ,$
(13)

where g represents a diagonal matrix with gp on the diagonal.

TABLE I.

Generalized Pauli constraints for $∧ 3 H 6$.

No. Constraint
g1 + λ4λ5λ6 = 0
λ1 + λ6 = 1
λ2 + λ5 = 1
λ3 + λ4 = 1
No. Constraint
g1 + λ4λ5λ6 = 0
λ1 + λ6 = 1
λ2 + λ5 = 1
λ3 + λ4 = 1

The GPC in Table I are expressed in terms of the natural orbital occupation numbers. However, v2RDM optimizations in general are not performed in the natural orbital basis. For example, the orbitals in a v2RDM-driven CASSCF computation are chosen to minimize the electronic energy, rather than to diagonalize the 1-RDM. We instead express the GPC in terms of the nondiagonal 1-RDM and its eigenvectors. Then, for example, constraint 1 in Table I becomes

$g 1 + ∑ i j [ U i 4 U j 4 − U i 5 U j 5 − U i 6 U j 6 ] 1 D j i = 0 ,$
(14)

where U is the unitary matrix that diagonalizes 1D. As linear equality constraints, the GPC can now be incorporated into a modified boundary-point semidefinite optimization in which the constraint matrix, A, changes as the natural orbitals evolve. For a given 1-RDM, we determine the natural orbitals and solve for the dual solution according to the expression in Step 1 above. The primal solution is updated in Step 2, new natural orbitals are determined, and the optimization continues. We have incorporated the GPC into a development version of our SDP solver for the ground-state v2RDM problem.15 The solver is implemented as a plugin to the Psi4 electronic structure package.36

The semidefinite optimization procedure outlined above yields 1- and 2-RDMs for the ground state of a given spin symmetry. Excited-state information can then be determined using an extended random phase approximation. Within the ERPA, the nth excited electronic state, |Ψn〉, is defined by single excitations out of the ground state:

$| Ψ n 〉 = ∑ i j c i j n a ˆ j † a ˆ i | Ψ 〉 .$
(15)

The expansion coefficients, $c i j n$, and the excitation energies, ωn, are determined by solving the generalized eigenvalue problem

$∑ i j c i j n 〈 Ψ | [ a ˆ k † a ˆ l , [ H ˆ , a ˆ j † a ˆ i ] ] | Ψ 〉 = ω n ∑ i j c i j n 〈 Ψ | [ a ˆ k † a ˆ l , a ˆ j † a ˆ i ] | Ψ 〉 .$
(16)

The left- and right-hand sides of Eq. (16) can expressed in terms of the ground-state 1- and 2-RDM, which can be obtained from a v2RDM computation as described above. In all computations presented below, the summation over spin-orbitals i and j is restricted to include only those orbitals that are active in the underlying v2RDM optimization. We have implemented the ERPA as a plugin to Psi4.

Consider the ground state of a boron atom (2P). This state is triply degenerate, so a v2RDM optimization under ensemble N-representability conditions is unlikely to yield pure-state RDMs. As discussed in Ref. 22, the ensemble nature of the RDMs has a significant impact on the quality of excitation energies extracted using the ERPA; the error in the 4P←2P excitation energy, relative to that from full configuration interaction, was reported in that work to be 2.674 eV. Here, we explore the role of the GPC in the ability of the ERPA to describe this transition in a series of five-electron atoms and ions. We compare the 4P←2P excitation energies computed by the v2RDM/ERPA approach to those obtained at the state-averaged (SA) CASSCF level of theory. Both the 2P and 4P states are triply degenerate, so the SA-CASSCF procedure involved all six states. In the CASSCF and v2RDM optimizations, the 1ag, 2ag, 1b1u, 1b2u, and 1b3u orbitals were active; the boron isoelectronic series thus belongs to the $∧ 5 H 10$ family. All computations were performed within the cc-pVQZ basis set, and the v2RDM procedure used SA-CASSCF-optimized orbitals. The v2RDM optimizations enforce the ensemble two-particle (PQG) N-representability conditions,4 as well as constraints on the expectation value of $S ˆ 2$.31,32

Figure 1 illustrates violations in the 160 GPC for $∧ 5 H 10$ systems in units of electrons; expressions for these constraints can be found in the supplementary material of Ref. 24. The data were generated using two different initial choices for the 1- and 2-RDM. First, the v2RDM optimizations were seeded with Hartree-Fock RDMs [Fig. 1(a)]. For all species, many GPC are violated, and the greatest violations range from 0.3 electrons for Be to 5.9 electrons for Ne5+. Despite these large violations, in all cases, the v2RDM ground-state energies agree with those from SA-CASSCF to within 0.6 mEh. The violations in the GPC are even more apparent when seeding the v2RDM optimization with random positive semidefinite RDMs [Fig. 1(b)]. For all systems, some GPC are violated by at least 3.0 electrons, and the greatest error for Ne5+ is almost 7.0 electrons. It would be interesting to consider what physical meaning can be ascribed to GPC that are consistently violated here. However, the form of the GPC is so complicated that it is difficult to attach any physical meaning to them whatsoever. For example, the constraint that is violated to the greatest degree in Fig. 1(b) takes the form

$− 7 λ 1 + 3 λ 2 + 13 λ 3 + 13 λ 4 + 23 λ 5 + 33 λ 6 − 7 λ 7 − 27 λ 8 − 27 λ 9 − 17 λ 10 ≤ 45 .$
(17)

It is not immediately clear what useful physical meaning, if any, such a complex expression contains.

FIG. 1.

Violations (in units of electrons) in the generalized Pauli constraints for $∧ 5 H 10$ when seeding an ensemble v2RDM optimization with (a) Hartree-Fock and (b) random positive semidefinite RDMs. Constraints that are satisfied appear as black.

FIG. 1.

Violations (in units of electrons) in the generalized Pauli constraints for $∧ 5 H 10$ when seeding an ensemble v2RDM optimization with (a) Hartree-Fock and (b) random positive semidefinite RDMs. Constraints that are satisfied appear as black.

Close modal

Despite the sizeable differences in the degree to which the GPC are violated when using different initial conditions, ground-state energies obtained from v2RDM optimizations seeded with Hartree-Fock or random RDMs agree to within 0.01 mEh, which is consistent with the energy convergence criterion used throughout this work (|EprimalEdual| < 0.1 mEh). At first pass, it may come as a surprise that such strikingly different 1-RDMs can be associated with essentially identical energies, and it is tempting to attribute the differences to an incomplete or unsuccessful optimization. However, the differences simply reflect the ensemble nature of the optimized 1-RDMs. Recall that an ensemble N-representable 1-RDM can be obtained from the partial integration of an ensemble N-electron density matrix, similar to the integration given in Eq. (1). In the case of the triply degenerate 2P ground state for the boron series, an ensemble N-electron density matrix and the corresponding 1-RDM contain contributions from all three states. Changes to the initial conditions in the v2RDM optimization alter the relative weights of the states that contribute to the optimized ensemble 1-RDM.

Table II illustrates the 4P←2P excitation energies for the boron series computed at the SA-CASSCF and v2RDM/ERPA levels of theory, using the same set of active orbitals utilized in the ground-state computations. In describing the 4P states, the SA-CASSCF procedure considers the low-spin ($M S = 1 2$) 4P states, while the ERPA considers the high-spin ($M S = 3 2$) 4P states. That is, the ERPA expansion of the wave function in Eq. (15) was restricted to include only spin-forbidden excitations that increase the value of MS for the wave function. The rationale behind this discrepancy is practical. The SA-CASSCF procedure in Psi4 requires that each state has the same MS value. However, the low-spin 4P states will not be correctly described by the ERPA when considering only single excitations out of the ground state; the low-spin 4P states have significant double-excitation character. Similarly, even for the high-spin 4P states, we can only expect the ERPA to properly describe two of the three degenerate states. Consider the Hartree-Fock reference configuration for the b1u-symmetry 2P state of boron [(1ag)2(2ag)2(1b1u)1]. The b2g- and b3g-symmetry 4P states can be generated by single excitations out of this reference, but the b1g-symmetry 4P state differs from the reference by two electrons. Hence, Table II includes only two excitation energies from v2RDM/ERPA.

TABLE II.

4P←2P excitation energies (eV) for the boron isoelectronic series. Ensemble v2RDM optimizations were seeded with either (1) Hartree-Fock or (2) random positive semidefinite RDMs. Pure-state v2RDM optimizations were seeded with ensemble RDMs initially optimized from (3) Hartree-Fock or (4) random positive semidefinite RDMs.

ERPA
Ensemble Pure
1 2 3 4 SA-CASSCF
1.409  1.552  1.384  1.384
Be  1.409  1.622  1.384  1.384  1.352
3.107  3.312  3.039  3.030
3.107  3.413  3.039  3.030  2.991
4.889  5.126  4.735  4.747
C+  4.889  5.228  4.735  4.747  4.694
6.704  6.950  6.453  6.466
N2+  6.704  7.047  6.453  6.466  6.410
8.555  8.791  8.178  8.191
O3+  8.555  8.870  8.179  8.191  8.131
10.445  10.641  9.901  9.912
F4+  10.445  10.705  9.902  9.912  9.849
12.386  12.522  11.621  11.620
Ne5+  12.386  12.544  11.621  11.620  11.564
Maximum error  0.822  0.980  0.057  0.063
ERPA
Ensemble Pure
1 2 3 4 SA-CASSCF
1.409  1.552  1.384  1.384
Be  1.409  1.622  1.384  1.384  1.352
3.107  3.312  3.039  3.030
3.107  3.413  3.039  3.030  2.991
4.889  5.126  4.735  4.747
C+  4.889  5.228  4.735  4.747  4.694
6.704  6.950  6.453  6.466
N2+  6.704  7.047  6.453  6.466  6.410
8.555  8.791  8.178  8.191
O3+  8.555  8.870  8.179  8.191  8.131
10.445  10.641  9.901  9.912
F4+  10.445  10.705  9.902  9.912  9.849
12.386  12.522  11.621  11.620
Ne5+  12.386  12.544  11.621  11.620  11.564
Maximum error  0.822  0.980  0.057  0.063

The ERPA excitation energies were generated from four different solutions to the ground-state v2RDM problem. First, the 1- and 2-RDM were optimized under the usual PQG N-representability conditions with an initial guess for the RDMs of either the Hartree-Fock solution or random positive semidefinite matrices (labeled “ensemble”). The RDMs were also optimized subject to the PQG conditions and the 160 GPC for $∧ 5 H 10$ starting from either of the ensemble solutions (labeled “pure”). Ensemble RDMs optimized from a Hartree-Fock guess yield excitation energies that differ substantially (by up to 0.814 eV) from those from SA-CASSCF, but the ERPA at least correctly predicts that the 4P states are degenerate. These large errors are troubling in and of themselves, but what is more problematic is the degree to which the excitation energies change when seeding the v2RDM optimization with a different guess. In the case of a random seed, the errors are as large as 0.98 eV, and the ERPA no longer predicts that the 4P states are degenerate. This result is not too surprising, considering the degree to which violations of the GPC increase with the random initial seed [see Fig. 1(b)]. With the application of pure-state N-representability conditions, the errors in the ERPA excitation energies, relative to those from SA-CASSCF, fall to only 0.007-0.057 eV, and the 4P states appear as degenerate states in the spectra (to within 1 meV, in two cases). The ERPA excitation energies from pure-state v2RDM optimizations seeded with different initial RDMs differ by at most 13 meV (0.5 mEh); this difference can be traced to the convergence thresholds in the ground-state v2RDM optimization. Throughout this work, v2RDM optimizations were considered converged when ‖Axb‖ < 10−4, ‖ATyc + z‖ < 10−4, and |EprimalEdual| < 10−4Eh. With these thresholds, the 4P←2P excitation energies for boron using the ERPA with pure-state RDMs differ by 8 meV, depending on the initial optimization conditions. By decreasing the primal and dual convergence thresholds to 10−5 and the primal/dual energy gap to 10−5Eh, this difference reduces to only 1 meV.

We have also explored the role of the GPC in describing the 5S←3P excitation in the carbon (six-electron) isoelectronic series. Table III provides excitation energies for the series computed at the SA-CASSCF and v2RDM/ERPA levels of theory, using the cc-pVQZ basis set. The 3P state is triply degenerate, the 5S state is nondegenerate, and the SA-CASSCF procedure involved all four states. Again, the 1ag, 2ag, 1b1u, 1b2u, and 1b3u orbitals were active, and the v2RDM procedure used SA-CASSCF-optimized orbitals. As with the boron series, the ERPA considers the high-spin (MS = 2) 5S state, while the SA-CASSCF procedure considers the MS = 1 state. The carbon series belongs to the $∧ 6 H 10$ family. By particle-hole equivalence the constraints for $∧ 4 H 10$ were used, except that, in this case, the GPC apply to the eigenvalues of the one-hole RDM, 1Q, rather than those of the 1-RDM. Expressions for the 124 GPC for $∧ 4 H 10$ were taken from the supplementary material of Ref. 24. As in the boron series, ensemble RDMs display large errors, relative to SA-CASSCF (as large as 1.560 eV), and the quality of the excitation energies varies dramatically depending on the initial optimization conditions. This failure is particularly striking, considering the simplicity of the structure of the high-spin 5S state; with only 10 active spin orbitals, the wave function is a single determinant. When the RDMs are optimized under pure-state N-representability conditions, the errors in the ERPA excitation energies fall to 0.001 eV or less.

TABLE III.

5S←3P excitation energies (eV) for the carbon isoelectronic series. Ensemble v2RDM optimizations were seeded with either (1) Hartree-Fock or (2) random positive semidefinite RDMs. Pure-state v2RDM optimizations were seeded with ensemble RDMs initially optimized from (3) Hartree-Fock or (4) random positive semidefinite RDMs.

ERPA
Ensemble Pure
1 2 3 4 SA-CASSCF
Be2−  0.530  0.809  0.485  0.485  0.485
B  1.731  2.217  1.628  1.628  1.628
3.149  3.795  2.956  2.956  2.956
N+  4.772  5.486  4.465  4.464  4.463
O2+  6.510  7.244  6.050  6.050  6.049
F3+  8.331  9.043  7.670  7.669  7.669
Ne4+  10.238  10.864  9.304  9.304  9.303
Maximum error  0.934  1.560  0.001  0.001
ERPA
Ensemble Pure
1 2 3 4 SA-CASSCF
Be2−  0.530  0.809  0.485  0.485  0.485
B  1.731  2.217  1.628  1.628  1.628
3.149  3.795  2.956  2.956  2.956
N+  4.772  5.486  4.465  4.464  4.463
O2+  6.510  7.244  6.050  6.050  6.049
F3+  8.331  9.043  7.670  7.669  7.669
Ne4+  10.238  10.864  9.304  9.304  9.303
Maximum error  0.934  1.560  0.001  0.001

The groundbreaking work of Altunbulak and Klyachko24 defined a systematic procedure to derive pure-state N-representability conditions for arbitrary many-electron systems. This breakthrough has far reaching implications for quantum chemistry, facilitating fundamental investigations into the nature of electron correlation and providing, for the first time, some guidance as to how one could directly determine 1- and 2-RDMs for degenerate states that satisfy pure-state N-representability conditions. We have presented a procedure that incorporates these generalized Pauli constraints into the variational optimization of the ground-state 2-RDM, and the utility of pure RDMs over ensemble ones is clearly visible in the context of the ERPA. For many-electron systems with degenerate ground states, excitation energies obtained from ensemble RDMs are quite poor. However, when enforcing the generalized Pauli constraints, the resulting pure-state N-representable RDMs can safely be used within the ERPA to obtain reliable excited-state information.

The present work constitutes numerical proof that it is indeed possible to enforce the GPC within a v2RDM optimization. However, two extensions are necessary before the GPC can be applied to general systems. The GPC have only been tabulated for systems with 10 or fewer active orbitals;24 future work in this area must include the generation and tabulation of GPC for other systems with larger numbers of active electrons and orbitals. The number of GPC will increase sharply for larger systems, but this number will be significantly less than the number of constraints associated with other necessary N-representability conditions, such as the two-particle PQG conditions. Second, the application of the GPC results in a substantial increase in the complexity and time-to-solution of the present SDP algorithm. On average, the number of macroiterations (cycles of Steps 1 and 2 outlined in Sec. II) increases by a factor of 24 for the boron series ($∧ 5 H 10$) and 8 for the carbon series ($∧ 6 H 10$). This increase in the time-to-solution is not consistent across either series. For example, the total number of macroiterations for the Be atom increased by a factor of ≈100, whereas the number of iterations remained essentially constant for Ne5+. Clearly, the cost of application of pure-state conditions could be prohibitive for more complex systems. Future work must consider strategies to enhance the convergence of the v2RDM optimization procedure. One possible strategy would be to decouple the enforcement of the ensemble PQG and pure-state conditions.

This material is based upon work supported by the National Science Foundation under Grant No. CHE-1554354.

1.
P.-O.
Löwdin
,
Phys. Rev.
97
,
1474
(
1955
).
2.
A. J.
Coleman
,
Rev. Mod. Phys.
35
,
668
(
1963
).
3.
D. A.
Mazziotti
and
R. M.
Erdahl
,
Phys. Rev. A
63
,
042113
(
2001
).
4.
C.
Garrod
and
J. K.
Percus
,
J. Math. Phys.
5
,
1756
(
1964
).
5.
R.
Erdahl
,
Rep. Math. Phys.
15
,
147
(
1979
).
6.
M.
Nakata
,
H.
Nakatsuji
,
M.
Ehara
,
M.
Fukuda
,
K.
Nakata
, and
K.
Fujisawa
,
J. Chem. Phys.
114
,
8282
(
2001
).
7.
D. A.
Mazziotti
,
Phys. Rev. A
65
,
062511
(
2002
).
8.
D. A.
Mazziotti
,
Phys. Rev. A
74
,
032501
(
2006
).
9.
Z.
Zhao
,
B. J.
Braams
,
M.
Fukuda
,
M. L.
Overton
, and
J. K.
Percus
,
J. Chem. Phys.
120
,
2095
(
2004
).
10.
M.
Fukuda
,
B. J.
Braams
,
M.
Nakata
,
M. L.
Overton
,
J. K.
Percus
,
M.
Yamashita
, and
Z.
Zhao
,
Math. Program.
109
,
553
(
2007
).
11.
E.
Cancès
,
G.
Stoltz
, and
M.
Lewin
,
J. Chem. Phys.
125
,
064101
(
2006
).
12.
B.
Verstichel
,
H.
van Aggelen
,
D.
Van Neck
,
P. W.
Ayers
, and
P.
Bultinck
,
Phys. Rev. A
80
,
032508
(
2009
).
13.
B.
Verstichel
,
H.
van Aggelen
,
D. V.
Neck
,
P.
Bultinck
, and
S. D.
Baerdemacker
,
Comput. Phys. Commun.
182
,
1235
(
2011
).
14.
J.
Fosso-Tande
,
D. R.
Nascimento
, and
A. E.
DePrince
III
,
Mol. Phys.
114
,
423
(
2016
).
15.
J.
Fosso-Tande
,
T.-S.
Nguyen
,
G.
Gidofalvi
, and
A. E.
DePrince
III
,
J. Chem. Theory Comput.
12
,
2260
(
2016
).
16.
M.
Bouten
,
P.
van Leuven
,
M.
Mihailovi
, and
M.
Rosina
,
Nucl. Phys. A
221
,
173
(
1974
).
17.
M.
Rosina
,
Int. J. Quantum Chem.
13
,
737
(
1978
).
18.
C.
Valdemoro
,
D. R.
Alcoba
,
O. B.
Oña
,
L. M.
Tel
, and
E.
Pérez-Romero
,
J. Math. Chem.
50
,
492
(
2012
).
19.
L.
Greenman
and
D. A.
Mazziotti
,
J. Chem. Phys.
128
,
114109
(
2008
).
20.
D. A.
Mazziotti
,
Phys. Rev. A
68
,
052501
(
2003
).
21.
K.
Chatterjee
and
K.
Pernal
,
J. Chem. Phys.
137
,
204109
(
2012
).
22.
H.
van Aggelen
,
B.
Verstichel
,
G.
Acke
,
M.
Degroote
,
P.
Bultinck
,
P. W.
Ayers
, and
D. V.
Neck
,
Comput. Theor. Chem.
1003
,
50
54
(
2013
).
23.
R. E.
Borland
and
K.
Dennis
,
J. Phys. B: At. Mol. Phys.
5
,
7
(
1972
).
24.
M.
Altunbulak
and
A.
Klyachko
,
Commun. Math. Phys.
282
,
287
(
2008
).
25.
C.
Schilling
,
D.
Gross
, and
M.
Christandl
,
Phys. Rev. Lett.
110
,
040404
(
2013
).
26.
C. L.
Benavides-Riveros
,
J. M.
Gracia-Bondía
, and
M.
Springborg
,
Phys. Rev. A
88
,
022508
(
2013
).
27.
R.
Chakraborty
and
D. A.
Mazziotti
,
Phys. Rev. A
89
,
042505
(
2014
).
28.
C.
Schilling
,
Phys. Rev. A
91
,
022105
(
2015
).
29.
I.
Theophilou
,
N. N.
Lathiotakis
,
M. A. L.
Marques
, and
N.
Helbig
,
J. Chem. Phys.
142
,
154108
(
2015
).
30.
T. L.
Gilbert
,
Phys. Rev. B
12
,
2111
(
1975
).
31.
E.
Pérez-Romero
,
L. M.
Tel
, and
C.
Valdemoro
,
Int. J. Quantum Chem.
61
,
55
(
1997
).
32.
G.
Gidofalvi
and
D. A.
Mazziotti
,
Phys. Rev. A
72
,
052505
(
2005
).
33.
J.
Povh
,
F.
Rendl
, and
A.
Wiegele
,
Computing
78
,
277
(
2006
).
34.
J.
Malick
,
J.
Povh
,
F.
Rendl
, and
A.
Wiegele
,
SIAM J. Optim.
20
,
336
(
2009
).
35.
D. A.
Mazziotti
,
Phys. Rev. Lett.
106
,
083001
(
2011
).
36.
J. M.
Turney
,
A. C.
Simmonett
,
R. M.
Parrish
,
E. G.
Hohenstein
,
F. A.
Evangelista
,
J. T.
Fermann
,
B. J.
Mintz
,
L. A.
Burns
,
J. J.
Wilke
,
M. L.
Abrams
et al,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
556
(
2012
).