Precise Quantum Chemistry calculations with few Slater Determinants

Domain: physics
Published: 2026-03-16
Editorial outcome: accept

Abstract
Slater determinants have underpinned quantum chemistry for nearly a century, yet their full potential has remained challenging to exploit. In this work, we show that a variational wavefunction composed of a few hundred optimized non-orthogonal determinants can achieve energy accuracies comparable to the state of the art. This is obtained by introducing an optimization method that leverages the quadratic dependence of the variational energy on the orbitals of each determinant, enabling an exact iterative optimization, and uses an efficient tensor-contraction algorithm to evaluate the effective Hamiltonian with a computational cost that scales as the fourth power of the number of basis functions. We benchmark the accuracy of the proposed method with exact full-configuration interaction results where available, and we achieve lower variational energies than coupled cluster (CCSD(T)) for several molecules in the double-zeta basis.

Introduction
Slater Determinants (SDs) have been instrumental in shaping our understanding of quantum chemistry since their introduction nearly a century ago. In 1927, Heitler and London [1] explained the qualitative nature of the covalent bond in H using a wave function composed of a sum of non-orthogonal SDs built from single-atom orbitals. About two decades later, Coulson and Fischer [2] demonstrated that the Heitler-London wave function could be significantly improved by allowing atomic orbitals to hybridize and variationally optimizing them, thereby eliminating the need to introduce ionic structures which appear unphysical for H . Subsequent developments along this line, most notably the introduction of the Generalized Valence Bond Ansatz [3] and the Spin-Coupled wave function [4, 5], marked the foundation of modern Valence Bond theory [6, 7], and played a key role in our current understanding of chemical bonding. In Valence Bond theory, bonds emerge from resonating chemical structures, which are superpositions of basic wavefunctions and are themselves spin-determined linear combinations of non-orthogonal SDs constructed with semi-localized orbitals. The compactness of the wave function is crucial for this interpretation and can only be achieved through simultaneous optimization of both the structures and the orbitals. However, this remains a challenging task compared to molecular orbital based methods [8], which exploit orbital orthogonality to achieve unmatched numerical efficiency, at the expense of interpretability. Achieving chemical accuracy requires too many contributing SDs when built with orthogonal orbitals, while non-orthogonal SDs, potentially more interpretable, pose significant numerical challenges when high accuracy is required.

In view of the above, it is not surprising that many techniques have been proposed to optimize SDs. The Full Configuration Interaction (FCI) expansion, while exact within a chosen finite basis, scales exponentially with basis-set size, rendering it impractical for all but the very smallest systems. Approximate methods can be generally divided into two classes: (i) methods using one or multiple fixed-reference determinants that optimize a linear combination of excitations thereof, and (ii) multi-configuration determinantal methods that optimize the orbitals of the SDs, along with the configuration interaction (CI) coefficients.

A notable example of a method belonging to the first category is truncated CI, where excitations from a reference determinant are limited to a certain excitation level. For instance, only single and double excitations are included in CISD. In such cases, the number of determinants is combinatorial in the excitation level and CI truncation breaks size-consistency. The prohibitive scaling can be mitigated by selected-CI methods, which retain only a subset of relevant excitations. However, identifying the important determinants is generally a nontrivial task. Quadratic Configuration Interaction [9] and Coupled Cluster [10] methods such as CCSD(T) correct the size-consistency by approximately including also higher excitations; however, unlike CI methods, they are non-variational and come at an increased computational cost. The methods discussed so far all use orthogonal Slater determinants. The non-orthogonal configuration interaction technique [11], on the other hand, optimizes a linear combination of given, fixed, non-orthogonal Slater Determinants.

The second class of methods, which is the most relevant in spirit to the topic of this article, consists of a variety of multi-configurational approaches that employ different strategies to optimize both the CI coefficients and the orbitals of the SDs. Thanks to the improved quality of the SDs achieved with the orbital optimization, these CI expansions are typically much more compact than those in the first class, albeit at the cost of increased computational effort. Multi-configurational self-consistent field methods, such as the Complete Active Space Self Consistent Field Method [12], optimize a single orbital rotation applied to all or a subset of the CI determinants, and have recently been generalized to non-orthogonal orbitals [13]. Gradient-based approaches optimizing only the orbitals of a single determinant at a time have been proposed for two determinants [14], and many determinants [15], or Hartree-Fock-Bogoliubov wavefunctions [16, 17]. Resonating Hartree-Fock [18–20] and nonorthogonal multiconfigurational self-consistent field [21] exploit the Thouless parametrization of the SDs [22] to optimize all orbitals simultaneously using gradient-based methods. However, it can be difficult to make the process converge [17, 23, 24]. Valence Bond Self-Consistent Field methods [25, 26] use second-order optimizers to optimize a sum of non-orthogonal determinants. Multi-configuration time-dependent Hartree methods [27], extended to fermions [28] use numerical integration to apply the real or imaginary time evolution to the determinants.

SDs also play a crucial role in many-body methods such as Quantum Monte Carlo. Compact SD expansions serve as efficient trial wavefunctions in Auxiliary-Field Quantum Monte Carlo [29, 30] and Diffusion Monte Carlo [31], and underpin variational ansätze like the Slater-Jastrow form [32] as well as recent variants augmented with neural networks [33–36].

In this work we show that a sum of fewer than a thousand SDs is sufficient to accurately represent the all-electron ground-state wavefunction of small-size molecules, when the SDs are fully optimized without orthogonality constraints. These results are achieved with a computational optimization technique that we introduce in this article, which we call the Exact Iterative Determinant-Orbital Solver (EIDOS). EIDOS iteratively optimizes the variational energy when viewed as a function of a string of orbitals, by selecting one orbital from each determinant and solving a generalized eigenvalue problem at every step. It relies on: (i) the direct parameterization of the molecular orbital coefficients, as opposed to the popular Thouless parameterization of the SDs [22]; (ii) an efficient analytical evaluation of the Hessian of the energy functional, resulting in a favorable O(m4) scaling with the basis set size m (see the Methods Section for details). We assess the accuracy of EIDOS on several molecules in the cc-pVDZ basis, benchmarking with full-CI data where available, and with DMRG and coupled-cluster (CCSD(T)) otherwise. We find that our variational energies are consistently below coupled-cluster, and within chemical accuracy of full-CI and DMRG, demonstrating that EIDOS is competitive with state-of-the-art methods, but with significantly more favorable scaling. We additionally discuss how the number of determinants scales with bond order, and we demonstrate convergence toward the complete-basis-set limit for LiH. To test the method’s robustness, we analyze the strongly-correlated regime along the N dissociation curve, and we verify that it correctly reproduces the triplet ground state of O .

Results
Non-orthogonal Slater-Determinant Ansatz optimized by the EIDOS algorithm. The wave function we consider consists of a sum of unrestricted unnormalized non-orthogonal SDs

| Φ¯ ⟩ := (cid:88) ND I=1 Φˆ(I↑)Φˆ(I↓) | 0 ⟩ , (1)

where N is the number of determinants, | 0 ⟩ is the vacuum state, and

SΦˆ(I,σ) := (cid:89) nσ i=1 (cid:88) m µ=1 ϕ(I,σ) i,µ cˆ† µ,σ , (2)

where ϕ(I,σ) are independent orbitals, n is the number of electrons of spin σ, m is the number of basis orbitals, and cˆ† is the creation operator of an electron in basis orbital µ with spin σ.

We optimize the average energy of the non-orthogonal Slater determinant wave-function ansatz of Eq. (1) using our novel EIDOS algorithm. Below, we briefly outline the main ideas underlying EIDOS, and refer the reader to the Methods section for the details. The efficiency of EIDOS relies on two key aspects: (i) the average energy of the wavefunction of Eq. (1) is the ratio of two quadratic functions of ϕ(1,σ),...,ϕ(ND,σ) at fixed electron label i = 1, which implies that these parameters can be optimized by solving a series of generalized eigenvalue problems, where the procedure is repeated after a random quantum-state-preserving electron-label transformation in each determinant; (ii) the needed effective quadratic forms can be computed with a total computational cost of O(m4) with an efficient tensor-contraction strategy.

Benchmark with CCSD(T), DMRG and FCI for equilibrium geometries in the cc-pVDZ basis. To validate the accuracy of the proposed EIDOS method, we compare it against CCSD(T), and, where available, DMRG and FCI, for the all-electron ground-state energy of several molecules at equilibrium in the cc-pVDZ basis set. In our calculations, the electron number varies from 4 (LiH) to 20 (for LiCl), except for the Krypton atom with 36 electrons in 27 orbitals, where we perform a particle-hole transformation to reduce the effective number of interacting particles to 18. As reported in Fig. 1, the EIDOS optimized energies from up to N = 768 Slater determinants (N = 2048 is obtained with the extension described in Supplementary Note 3) are well within chemical accuracy (1 kcal/mol) with respect to FCI, including those at the limit of what is still tractable by current distributed FCI codes [37]. Furthermore, our variational energies are consistently lower than those obtained using the coupled cluster method at the CCSD(T) level. The latter is considered the gold standard for molecules near equilibrium with minimal multireference character, where dynamic electron correlation dominates, which is the case for a large part of the molecules we consider here. These results are significant because CCSD(T) has an asymptotic computational scaling of O(m7), while our method scales like O(m4).

For all molecules we find the expected spin symmetry for the ground state, with a maximal spin contamination at 768 determinants of ⟨ Sˆ2 ⟩ S(S+1) 5 10−3 ( 9 10−4 with 2048 determinants) for BN, and less for the other molecules. We also note that 768 determinants consistently improve the results of 512 determinants, and that the optimization steps of the EIDOS algorithm to achieve convergence are roughly independent of the number of determinants (see Supplementary Note 1).

Scaling of the number of SDs with bond order. In general, we expect the number of SDs needed to precisely model the ground state to increase with correlations and, in particular, we expect a dependence on the number of valence electrons. To verify this picture in a controlled setting, we consider the three molecules LiF, BeO and N , which all have a similar number of electrons, n 12,14 , the same number of basis orbitals in the cc-pVDZ basis set, m = 28, but different bond order. In Fig. 2 we report the EIDOS energy error w.r.t. FCI as a function of the number of determinants, ranging from N = 1 up to N = 768, observing a smooth and consistent scaling. Moreover, the log-log scale we use in the plot graphically shows that the decay of the energy error is certainly not slower than a power-law of the inverse number of SDs. The slope of the decay clearly decreases with the bond order/number of valence electrons, meaning that for N we need more than an order of magnitude more determinants to attain chemical accuracy than what is needed for LiF, with BeO lying in between. The single-bonded BH shows a similar scaling as LiF, while it is interesting to remark that the curves for C and BN, which are considered to be doubly-bonded molecules, nearly collapse on the curve of the triple-bonded N .

Convergence to the large-basis-set limit. In Fig. 3 we consider the ground-state energy of the LiH molecule at equilibrium distance obtained with the EIDOS technique for various basis sets of sizes ranging from m = 6 to m=146 orbitals, the largest for which FCI energies are still available [37]. The approximate curve collapse for the energy error in a given basis set as a function of the ratio of the basis set size m and the number of SDs (shown in the inset of Fig. 3) suggests that the number of determinants needed to get a fixed error in a given basis scales linearly with the basis set size. We conjecture that this might be due to the short-distance cusp conditions of the continuum-space wave function, which require having an increasing number of determinants in the approach to the continuum limit. Nevertheless, as we show in the main plot of Fig. 3, the absolute energy decreases with increasing basis set size m for a fixed number of SDs N , as does the error with respect to the true ground-state energy in the infinite-basis-set limit.

Capturing static correlations. As a prototypical example of bond breaking, we consider the Nitrogen molecule. As its bond is stretched, the wavefunction acquires an inherent multireference character. Single-reference methods such as CISD and CCSD(T) have serious difficulties there, converging more slowly or not converging at all [38]. In Fig. 4 we show the EIDOS dissociation curve of the N molecule in comparison to UCISD (which uses 30724 orthogonal determinants), UCCSD(T) and DMRG. We outperform UCISD and UCCSD(T) already with a few hundred independent Slater determinants. As can be seen in the inset, close to equilibrium, we obtain energies within chemical accuracy from DMRG, while deviations become larger at extended bond lengths, where static-correlation effects emerge. These effects are, however, less severe than in CCSD, which inherits them from the UHF reference, and we do not rule out that they could be further alleviated with an appropriate choice of the initial guess.

Triplet and singlet dioxygen states. As an additional application, we discuss the low-energy states of molecular oxygen for different sectors of the M and S quantum numbers. While ⟨ Sˆ2 ⟩ is directly tuned to the value of S in this calculation, we can optionally add a λSˆ2 penalty term to the Hamiltonian at no extra cost during the optimization. In Fig. 5 we study the convergence of the EIDOS method to the O ground state, both in the M =0 and M =1 spin sectors, which is found to be a triplet in both sectors, as expected. For the M =1 sector we need fewer determinants to accurately represent the wavefunction due to the smaller size of the subspace. By adding a λSˆ2 penalty term, with λ = 0.1, we can converge to a singlet excited state. The inset of Fig. 5 shows the absence of spin contamination, apart from the lowest number of determinants in the M =0 sector.

Methods
Notation. We use capital letters for many-body quantities: Φ¯ , Φ(I) denote n-electron wavefunctions, N is the number of determinants, and I,J,K, 1,...,N are determinant labels. For single-electron quantities we use lowercase letters: n is the number of electrons, i,j, 1,...,n are electron labels, m is the number of basis functions, and µ,ν, 1,...,m are bare basis-orbital labels. cˆ† is the creation operator of an electron in the µ-th basis state. For simplicity, we assume the basis to be orthonormal; this implies, in particular, that the usual anticommutation relations { cˆ µ, cˆ ν } = { cˆ† µ, cˆ† ν } = 0, { cˆ µ, cˆ† ν } = δ µ,ν hold. We consider the following wave-function ansatz

Φ¯ = (cid:88) ND I=1 Φ(I) , (3)

where { | Φ(1) ⟩ , | Φ(2) ⟩ ,..., | Φ(ND) ⟩ } is a set of non-orthogonal SDs, i.e. ⟨ Φ(I) | Φ(J) ⟩ for I = J is not necessarily zero. Each Slater determinant | Φ(I) ⟩ is fully specified by its own sequence of orbitals ( ϕ(I) 1 , ϕ(I) 2 ,..., ϕ(I) n ) , given as a column vector of coefficients, ϕ(I) i ∈ Cm. We can then build each SD | Φ(I) ⟩ from the orbital sequence as

| Φ(I) ⟩ = (cid:32) (cid:88) m µ1=1 ϕ(I) 1,µ1 cˆ† µ1 (cid:33) ··· (cid:32) (cid:88) m µn=1 ϕ(I) n,µn cˆ† µn (cid:33) | 0 ⟩ . (4)

We stress that the orbitals ϕ(I) i ∈ Cm are completely arbitrary: we do not impose collinearity or orthogonality within a single determinant or with respect to the orbitals of different determinants. In this parametrization, in contrast to the Thouless parametrization, the determinant is multilinear in the orbital coefficient vectors ϕ(I), and therefore linear when viewed as a function of one of them. Singling out an orbital, which we choose to be the first one for notational simplicity and without loss of generality, we can write

| Φ(I) ⟩ = (cid:88) m µ1=1 ϕ(I) 1,µ1 cˆ† µ1 | Φ(I) 1+ ⟩ , (5)

where “hole-SD” | Φ(I) 1+ ⟩ := (cid:80) µ2...µn ϕ(I) 2,µ2 cˆ† µ2 ··· ϕ(I) n,µn cˆ† µn | 0 ⟩ . Similarly, ⟨ Φ(I) | Φ(J) ⟩ = 0, i.e. the hole-SDs are defined with the first electron removed. We can write the full wavefunction in terms of hole-SDs as

Φ¯ = (cid:88) ND I=1 (cid:88) m µ1=1 ϕ(I) 1,µ1 cˆ† µ1 | Φ(I) 1+ ⟩ . (6)

The Exact Iterative Determinant-Orbital Solver (EIDOS). In this section, we present our method to variationally optimize the average energy of an ansatz consisting of a sum of SDs by solving at each step a linear algebra problem. We first show why it is enough to solve a generalized eigenvalue problem, and afterwards how to efficiently compute the terms appearing in the linear algebra problem.

Our first finding is that the energy expectation value with the wave-function ansatz of Eq. (3), when regarded as a function of a single orbital of each SD, is simply a ratio of two quadratic forms

E = ⟨ Φ¯ | Hˆ | Φ¯ ⟩ / ⟨ Φ¯ | Φ¯ ⟩ = v† H v / v† S v , (7)

where v ∈ C(ND×m), defined by

v Iµ := ϕ(I) 1,µ , (8)

is the “vector” with flattened index Iµ, for I ∈ { 1,...,N } and µ ∈ { 1,...,m } , containing the orbitals that were factored out from each determinant, and H , S ∈ C(ND×m)×(ND×m) are the “effective matrices” for respectively the Hamiltonian and the norm of the wave-function, and they are given by

H Iµ,Jν = ⟨ Φ(I) 1+ | cˆ µ Hˆcˆ† ν | Φ(J) 1+ ⟩ , (9)

S Iµ,Jν = ⟨ Φ(I) 1+ | cˆ µ cˆ† ν | Φ(J) 1+ ⟩ , (10)

where I,J ∈ { 1,...,N } and µ,ν ∈ { 1,...,m } . The hole-wavefunction matrix elements ⟨ Φ(I) 1+ | cˆ µ Hˆcˆ† ν | Φ(J) 1+ ⟩ needed for H can be calculated in a numerically-exact way at an asymptotic cost of O(N2 m 4 ), including the singular case where the hole-SDs are linearly dependent. Similarly, ⟨ Φ(I) 1+ | cˆ µ cˆ† ν | Φ(J) 1+ ⟩ can be calculated at a cost of O(N2 m 2 ). Therefore, asymptotically, calculating H is not more expensive than the energy expectation value.

The energy Eq. (7) takes the form of a generalized Rayleigh quotient, and can be minimized exactly w.r.t. v by solving the generalized eigenvalue problem (GEV)

H v = ϵ S v , (11)

with a computational cost of O(N D 3m3) operations. We remark that S is singular, with null space spanned by the orbitals of the hole-SDs | Φ(I) 1+ ⟩ . By projecting out the kernel of S , the problem can be reduced from N D m to a smaller eigenvalue problem of size N D (m n+1).

The procedure described so far minimizes the energy with respect to a single orbital of each SD. The simplest algorithm that one can devise is an iterative one, where one selects a different orbital at every step. After optimizing an orbital, the wavefunction reaches a local minimum with respect to the coefficients of this orbital. However, as soon as another orbital is optimized, this is in general no longer true for the previously optimized one, requiring re-optimization of all orbitals until a stationary point is reached. The selection of orbitals for this can be done deterministically through sweeps, or stochastically by independently selecting the orbitals of each determinant in e.g. a uniform way. Here we consider a more general stochastic approach where a random linear combination of the orbitals of the previous step is optimized. The SD representation ( ϕ(I) 1 , ϕ(I) 2 ,..., ϕ(I) n ) in Eq. (4) is not unique, and the different redundant representations are related through the special linear group SL(n,C), i.e. the group of square matrices with determinant 1. We can make use of this by stochastically sampling V(I) ∈ SL(n,C) independently for each determinant I and transforming the orbitals as

ϕ(I) i (cid:55)→ (cid:88) n j=1 V(I) ij ϕ(I) j , (12)

which leaves the quantum state unchanged. Simple sweeps and uniform selection are included as special cases. The complete iterative algorithm is sketched in Box 1.

We discuss an efficient generalization for the case when the wavefunction is an eigenstate of Sˆ z (Eq. (1) in the main text) in Supplementary Note 2. For a single determinant this is equivalent to an unrestricted Hartree-Fock calculation. We remark that the method cannot be used to optimize a restricted Ansatz (Φˆ(I↑) ≡ Φˆ(I↓)), since the resulting wave function is not multilinear. In Supplementary Note 3 we discuss an extension of the method to reduce the size of the eigenvalue problem based on the selective optimization of orbitals of only a subset of all determinants at every step, which is used to obtain the 2048-determinant results of Fig. 1.

Efficient calculation of the effective matrices. We briefly discuss how to compute the matrices H and S appearing in the EIDOS loop in an efficient and numerically stable way. For notational simplicity, in this section, we drop the determinant index I, and we denote by | Φ ⟩ := | Φ(I) ⟩ , | Ψ ⟩ := | Φ(J) ⟩ two generic Slater determinants, and we accordingly define ϕ := ϕ(I) j and ψ := ϕ(J) j . As it is well known, we can write a generalized Wick’s theorem for the overlap of two SDs [44]

⟨ Φ | Ψ ⟩ = det[A] , (13)

where A = ϕ†ψ = (cid:80)m µ=1 ϕ∗ i,µ ψ j,µ is the matrix of pairwise overlaps of the orbitals. We consider a generic number-conserving normal-ordered K-body operator

Wˆ := (cid:88) ξ1...ξK ζ1...ζK wζ1...ζK cˆ† ζK ...cˆ† ζ1 cˆ ξ1 ...cˆ ξK , (14)

and compute the matrix element

⟨ Φ | Wˆ | Ψ ⟩ = (cid:88) ξ1...ξK ζ1...ζK wζ1...ζK ρζ ξ 1 1 ... ζ ξ K K , (15)

where

ρζ ξ 1 1 ... ζ ξ K K := ⟨ Φ | cˆ† ζK ...cˆ† ζ1 cˆ ξ1 ...cˆ ξK | Ψ ⟩ (16)

is the K-body reduced transition density matrix. Analytical expressions for the one, two and K-body reduced (transition) density matrices and matrix elements of non-orthogonal Slater determinants are well-documented in the literature [11, 15, 18, 24, 44–55]. Here, we derive the versions needed for our purpose.

We start with the simplifying assumption that the overlap ⟨ Φ | Ψ ⟩ is nonzero (the matrix A is non-singular). In this case, we can write down a generalized Wick’s theorem for the K-body reduced transition density matrix,

ρζ ξ 1 1 ... ζ ξ K K = det[A] det [ Bξ1 ζ1 ... Bξ1 ζK ... BξK ζ1 ... BξK ζK ] , (17)

where Bξ ζ = (cid:80) i,j ψ ∗ i (A−1) ij ϕ j,ζ . We remark that a more general expression for Hartree-Fock-Bogoliubov wavefunctions has been derived [56, 57], which reduces to Eq. (17) for the special case of SD. For completeness we include a derivation of Eq. (17) in Supplementary Note 4. Eq. (17) has an overall computational scaling of O(K3m2K) when each determinant of the K × K matrices containing the B’s is computed explicitly (see e.g. Ref. [53, 55]).

Efficiently-contractable formula. In this work we only need to consider K ∈ { 0,1,2,3 } , for which it is computationally feasible to expand the determinants of Eq. (17) into their K! terms, and directly contract them with w,

⟨ Φ | Wˆ | Ψ ⟩ = (cid:88) ξ1...ξK ζ1...ζK wζ1...ζK (cid:88) σ∈SK sgn(σ) (cid:89) K t=1 Bξt ζσ[t] , (18)

where S K is the symmetric group and σ[t] is the t-th element of the permutation σ. Whenever wζ1...ζK ξ1...ξK admits a compact tensor decomposition, we can expect the application of Eq. (18) to be advantageous. To give an example of this fact, let us consider a factorizable 3-body operator given by wζ1ζ2ζ3 ξ1ξ2ξ3 = wζ1ζ2 ξ1ξ2 wζ3 ξ3 , which is exactly the form we encounter in Eq. (20). In this case, an application of the efficiently-contractable formula, Eq. (18), requires only O(m4) operations, whereas using Eqs. (15) and (16) needs O(m6) operations.

For the zero-overlap case ⟨ Φ | Ψ ⟩ = det[A] = 0, the standard approach taken is to do a singular value decomposition (SVD) of A [11, 15, 21, 24, 53–55, 57, 58], either generically or using optimized equations for each value of the rank of A. Here we follow the former approach. Taking A = UΛVH where Λ is diagonal and U,V are unitary, we write down the equation of the matrix element of a generic normal-ordered K-body operator, given by

⟨ Φ | Wˆ | Ψ ⟩ = det[U]det [ VH ] (cid:88) ξ1...ξK ζ1...ζK wζ1...ζK (cid:88) σ∈SK sgn(σ) (cid:89) K t=1 X ξt,kt Y kt,ζσ[t] , (19)

where X ξ,k = (cid:80) i ψ i ∗ ,ξ U i,k and Y k,ζ = (cid:80) j V j ∗ ,k ϕ j,ζ . The tensor T k1...kK is the product of all singular values which are not indexed, defined as T k1...kK = (cid:81) i∈{1...n} \ {k1...kK} Λ i . We remark that this equation is correct regardless of the singular values. The contraction of Eq. (19) has a computational cost larger by a factor of O(n) compared to Eq. (18).

For the matrices H and S , defined in Eq. (9) and Eq. (10), we need to calculate matrix elements of the form ⟨ Φ | cˆ µ Wˆcˆ† ν | Ψ ⟩ for all µ,ν ∈ { 1,...,m } where Wˆ is a K-body operator (see Eq. (14)) and thus cˆ µ Wˆcˆ† ν is a K+1-body operator. This explains why for a two-body Hamiltonian we only need to consider K ∈ { 0,1,2,3 } . In the following, we show that the calculation of ⟨ Φ | cˆ µ Wˆcˆ† ν | Ψ ⟩ for ΦWˆ ∈ Ψ { . } µ σ ν µ σ ν σ ξ ′ ζ ⟨ | | ⟩ can be brought into normal order, resulting in a sum of K+1, K and K 1-body operators. For K ≥ 1, in the nonzero-overlap case ⟨ Φ | Ψ ⟩ ̸ = 0 for which we can use Eq. (18), we have that

⟨ Φ | cˆ µ Wˆcˆ† ν | Ψ ⟩ = ⟨ Φ | Ψ ⟩ (cid:88) σ∈SK+1 sgn(σ) (cid:89) K+1 t=1 Bξt ζσ[t] + δ µ,ν (cid:88) σ∈SK sgn(σ) (cid:89) K t=1 Bξt ζσ[t] + (cid:88) u∈{1,...,K} (cid:88) σ∈SK sgn(σ) (cid:89) K t=1 Bξt ζσ[t] + (cid:88) v∈{1,...,K} (cid:88) σ∈SK−1 sgn(σ) (cid:89) K−1 t=1 Bξt ζσ[t] + (cid:88) u,v∈{1,...,K} (cid:88) σ∈SK−1 sgn(σ)( 1)u+v Bξt ζσ[t] , (20)

where we take ξ K+1 ≡ µ and ζ K+1 ≡ ν for the first term, ζ u ≡ ν for the third and ξ v ≡ µ in the fourth, and use Einstein summation for the remaining ξ and ζ indices. The tensors in the equation can be contracted in such a way that the computational cost is O((K+1)!m2K). We point out that for the special case K = 2 this scaling has been shown using low-rank perturbations by Ref. [19]. It can be shown that also in the zero-overlap case, using Eq. (19), the tensors can be contracted in O((K+1)!m2Kn). We provide the full equation in Supplementary Note 6.

For the special case K = 0, relevant for the calculation of the S matrix, we get

⟨ Φ | cˆ µ cˆ† ν | Ψ ⟩ = ⟨ Φ | Ψ ⟩ (Bµ ν + δ µ,ν) , (21)

which can be computed in O(m2).

We use Eq. (20) and Eq. (21) to compute the matrices Eq. (9) and Eq. (10). Given a generic two-body Hamiltonian Hˆ, this results in a computational cost of O(N2m4), whenever the overlap is nonzero, and O(N2m4n) otherwise.

Implementation details We work with the second quantized formulation of the electronic molecular Hamiltonian in a finite basis set of size m, given by

Hˆ = (cid:88) µν h µν cˆ† µσ cˆ νσ + 1 2 (cid:88) µνξζ h µνξζ cˆ† µσ cˆ† ξσ′ cˆ ζσ′cˆ νσ +E nuc , (22)

where µ,ν,ξ,ζ ∈ { 1,...,m } and σ,σ′ ∈ { ↑,↓ } . We get the values for the one- and two-electron integrals, h µν , h µνξζ , and E nuc , from PySCF [59]. We choose the molecular orbital basis obtained from a Hartree-Fock calculation as an orthonormal basis, and use complex-valued Slater determinants in double precision unless otherwise specified. We remark that the O(m5) scaling of transforming the two-electron integral explicitly into an orthogonal basis could be straightforwardly avoided by including it in Eq. (18) and Eq. (19) and by contracting the R tensors in O(m4).

To improve the overall numerical stability of EIDOS, we normalize the determinants | Ψ(I) ⟩ for computing N and H S . We remark that this can always be done without loss of generality, as the “CI-like” linear variational problem of optimizing the coefficients d I of | Ψ¯ ⟩ = (cid:80) d I | Ψ(I) ⟩ is a linear variational problem, and rescaling the determinants has no effect on its solution. Furthermore, for additional numerical stability, also without loss of generality, we orthogonalize the orbitals of every determinant | Ψ(I) ⟩ separately at every step, using a QR decomposition of the non-square matrices ϕ(I) 1 ϕ(I) 2 ... ϕ(I) n ∈ Cm×n, leaving the state invariant. This is also a transformation that does not affect the solution of the generalized eigenvalue problem at each EIDOS step.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript «Precise Quantum Chemistry calculations with few Slater Determinants» presents a variational wave function composed of nonorthogonal determinants which can compactly represent electron correlation effects. In general, I think the use of nonorthogonal Slater determinants is interesting, given that several previous works also indicates that nonorthogonal determinants can provide compact representations of the wave function. The code accompanying the paper was nicely set up and easy to use. However, I have some general and specific remarks regarding the paper.\nGeneral remarks:\n1) The idea of using nonorthogonal determinants is not new, as also indicated by the authors since they cite some of the literature. Except that the authors seem to have an efficient implementation, I can’t see that this approach offers something that the other methods do not offer. If so, I think this needs to be adressed more clearly in the paper.\n2) Although the boundary between dynamic and static correlation is a little bit blurry (especially along a dissociation curve) it is difficult to assess to which extent it captures static and/or dynamic correlation. In particular, in Fig 1 everything is compared to CCSD(T), indicating that it is dynamical correlation which is the target, since CCSD(T) is not suitable for static correlation. For some of the molecules in Fig 1, the results are closer to CCSD(T) than DMRG and FCI. However, the paper only talks about static correlation.\n3) With respect to static correlation, I don’t understand why CCSD(T) (as highlighted in the abstract and discussion) is used as a reference for comparison, It is well known that low order truncations of CC are only suitable for capturing dynamical correlation and is not suitable for bond breaking (static correlation). Hence, since «Capturing static correlations» is mentioned explicitly, using CCSD(T) as a reference is the wrong choice. CCSD(T) contain dynamical correlation, which is not suitable for the stretching of the bond whereas DMRG only contain static correlation. Hence, for Fig 4, when the error wrt DMRG is already at milliHartree level (Nsd=768) that is just the error in the static correlation, and it does not include the dynamical correlation error (which I suppose can be large for stretching such an electron rich bond).\nSpecific remarks:\n1) The random rotation of the orbitals: Is it correct to assume that there is no initial guess of the determinants, since the orbitals are randomly rotated? And since it is called a rotation, I suppose the transformation matrix in eq 12 is unitary. If so, this should be specified clearly. Also, I cannot find it clearly explained why one does the random rotation of the orbitals, and the magnitude of such a rotation (given it is unitary, there is some in the rotations). Further, given that nonorthogonal orbitals are used, the metric of the start guess orbitals is relevant. If the transformation matrix is unitary, the metric is will be conserved.\n2) The molecular orbitals: I am not sure why the authors use quotation marks when writing molecular orbitals. As far as I can see, they are just normal nonorthogonal molecular orbitals (unless I am mistaken?), and there is no need for the use of quotation marks. Also, this is not important at all: but MO stands for «molecular orbitals», so there is no need to talk about «MO orbitals», one can merely use the term «MOs».\n3) On page 2 it is written «EIDOS exactly optimizes the variational energy..». I don’t understand what is meant by exactly. Is it that the algorithm will only converge to minima (as opposed to saddle points)? And if so, how is this ensured?\n4) In the next sentence on page 2 it says «EIDOS relies on the direct parametrization of the molecular orbital coefficients- as opposed to the popular Thouless parametrization of the SDs». I think this needs some clarfication. Acting with a Thouless transformation (e.g. exponentially parametrized unitary orbital rotation operator) on a Slater determinant yields a parametrization of the orbitals coefficients directly (equal to that of eq 12 in manscript).\n5) Regarding Thouless parametrization of the states, there is also a statement on page 2 regarding using gradient-based methods and difficulties to converge. I think it the second order approaches of wave function optimization (also using Thouless parametrizations) could be mentioned.\n6) Under «Methods» the authors write «(assumed here, for simplicity, to be orthonormal for the usual anticommutation relations to hold)» It is unclear why this assumption is needed. If it is for the purpose of rearranging creation operators entering the SD definitions, this assumption is not needed. Creation. operators always anticommute among themselves.", "reviewer_2": "The authors introduce their EIDOS method, a (as far as I know) novel algorithm for variational optimization of non-orthogonal Slater determinant expansions, variationally optimizing both the orbitals and CI amplitudes simultaneously in a greedy quasi-randomized approach. They demonstrate that using a few hundred such fully optimized determinants yields ground-state energies for small molecules with accuracy rivalling or surpassing standard methods (e.g., CCSD(T)), at a polynomial computational cost (for a fixed determinant number). The history of non-orthogonal CI expansions is long, and there have been many approaches suggested previously. The high accuracy of these expansions for small molecules is well known, although they are largely niche approaches due to the lack of size extensivity and (exponential?) scaling with number of determinants.\nThe main success of the work to be celebrated is the development of the algorithm for a full quasi-randomized variational optimization of the orbitals of the determinants. It is a neat approach which linearizes the determinants based on random rotations of the variations of a single degree of freedom. This is a significant achievement, and the ideas of this optimization approach could find broader impact.\nHowever, the results are only presented up to a maximum of 20 electrons, and while accurate, does not currently enable calculations which cannot be achieved with alternate state of the art methods (e.g. DMRG or selected-CI). This is a high bar for the initial introduction of a new algorithm, and would still be happy to see it published as a novel algorithm, but with some more results detailing the scope of systems which the method could treat, and a wider discussion about the scaling of the algorithm with system size. In particular, while the lack of size consistency is noted, this major limitation (and the related lack of size extensivity) is largely glossed over. I would suggest this needs further scrutiny for a more honest appraisal of the potential of the method and its prospects. While the SI has a plot of size consistency for dissociated N2, more instructive would be the scaling of the error with respect to system size (i.e. characterizing the size extensivity error) – do we expect to require an exponential number of determinants for fixed relative error? A formal exponential number of determinants is not a disqualifying characteristic either (see e.g. DMRG and sCI) – if this has a small prefactor/exponent, then this can still be a useful method, and at least clarify the domain of applicability of this tool.\nI would suggest e.g. a series of increasing length chain of Hydrogen atoms at a fixed interatomic spacing, and to observe the error compared to DMRG for fixed determinant number as additional plots in the main text (and at changing bond lengths). Presumably the correlation energy will converge to zero, but at what rate? How does this compare to CISD or CCSD (the latter being size extensive)? How does the number of determinants need to scale to counteract this? This would be more instructive than the bond order or basis set plots, the latter of which in particular is not expected to scale exponentially. This plot may also clarify the potential of the approach to extend to results beyond what is possible with DMRG or other systematically improvable methods, even if not possible with the current implementation.\nOther recommended changes:\n- Fig. 1 should compare results to DMRG rather than CCSD(T). All systems presented should be able to be solved essentially exactly with DMRG, and so would make a far better reference.\n- Only a couple of sentences are given in the main text (not in the methods section) to describe the algorithm, and this is rather opaque. The methods section is good for describing the technical details, but given that the algorithm itself is the main novelty of the approach and major scientific advance, I would hope for a little more in the main text to describe it qualitatively.\n- Furthermore, it would be good to have more insights into the convergence of the algorithm, especially in tricky cases. For instance, this is a quasi-randomized greedy optimization for a non-convex optimization problem (characterized as ‘partial-but-exact’ in the main text, which I found particularly uninformative!), so presumably there are cases where it finds local minima, exhibits a slow-down in convergence, or is found particularly sensitive to the initial conditions (perhaps e.g. where there are low-lying excited states?). It would be good to include some further insights into the numerical performance of the algorithm beyond Fig. B1, also commenting on the rate of convergence for different cases and starting point (e.g. keep the same random rotations, but reorder the initial orbitals?). How does this also vary with system size? Where is this algorithm expected to struggle to optimize to the true solutions compared to alternate iterative methods – are there any insights here?\n- It would be good to compare the approach to modern approaches to (generalized) non-orthogonal CI, or fragment-localized optimization algorithms.\n- There are a few minor typographical errors and inconsistent notation (e.g., missing or mismatched parentheses and indices) remain in the text—careful proofreading is advised.\nOverall, the manuscript introduces an interesting and innovative algorithm for optimizing NOCI expansions. However, it does not clearly demonstrate, either numerically or a theoretical path, towards how the approach can be extended to larger systems or achieve higher accuracy (either in absolute terms or as a function of wall time) compared to existing alternatives. I would be pleased to see this work published in Nature Communications if these points can be addressed. Specifically, a more thorough and candid assessment of the algorithm’s scaling with system size is needed, both in terms of computational cost and accuracy with respect to number of determinants. This will help clarify the potential niche and practical applicability of this method, and justify its relevance for publication in a broad-interest scientific journal.", "reviewer_3": "The concept of releasing the orthogonality restrictions to obtain more compact wave functions is a very attractive one and has received quite some attention over the last decade. In the here presented manuscript, a method is described that efficiently optimizes a wave function expressed as a linear combination of non-orthogonal Slater determinants. The method has the potential to be ground-breaking and the article can become a landmark paper in the field. There are a few things that need to be looked at to put the paper in the context of previous work (point 1); to discuss the stability of the algorithm when larger, more realistic systems are targeted (point 2); to get a feeling of how important non-orthogonality is (point 3); to address the importance of full orbital relaxation (point 4) and to clarify the way in which the orbitals are picked for optimization (point 5).\n1. A comparison with the work of Olsen is in place (JCP 143, 114102). In this study (and some follow-up papers) a procedure is described in which the wave function is optimized with respect to both the CI and orbital coefficients in which the configuration state functions are expressed in non-orthogonal orbitals.\n2. The use of iterative solvers seems to be an unavoidable step to surpass the memory bottleneck. Although not further worked out in the manuscript, this can be a critical step in the further development of the method and make it capable to use for other than the model molecules discussed here. With the increased number of degrees of freedoms due to the loss of orthogonality and the possible linear dependencies, the use of iterative solvers may degradate the stability of the optimization process. More information on the prospects of iterative solvers could be very helpful for other researchers in the field to judge the suitability of further developing this promising method.\n3. What is the degree of non-orthogonality among the different determinants? Maximum overlap, percentage of (nearly) orthogonal determinants, to name a few indicators. An overview of some numerics on this point is more than welcome.\n4. The method minimizes the average energy of the wave function expressed as a linear combination of non-orthogonal Slater Determinants. This is different than the approaches that adopt a \"different orbitals for different electronic configurations\" strategy. The latter includes full orbital relaxation for each electronic configuration (configuration state function, Slater determinant) in the wave function, while the here presented EIDOS variant only partially recovers orbital relaxation. This is especially relevant for core-hole states or for electronic states where charge transfer excitations contribute considerably to the wave function. In such cases a much larger number of SD will be needed to obtain a balanced description of the different configurations. Is is possible to design a strategy to perform independent optimizations of the most important electronic configurations and subsequently account for the interaction between these fully relaxed (non-orthogonal) states?\n5. What is the reason to randomly select the orbital that is singled out for optimization? Does this imply special advantages compared to a simple linear selection procedure? How can we be sure that all orbitals in the SD are at least optimized once? Isn't this procedure prone to get stuck in a local minimum? Furthermore, it is not at all clear how steep will be the increase in the number of iterations in this loop when bigger systems (with much more occupied orbitals) are studied."}, "author_response": "NCOMMS-25-34613-T 1 of 18\nResponse to the Reviewer 1\nComment 1a — Themanuscript“PreciseQuantumChemistrycalculationswithfewSlater\nDeterminants” presents a variational wave function composed of nonorthogonal determinants\nwhich can compactly represent electron correlation e ects. In general, I think the use of\nnonorthogonal Slater determinants is interesting, given that several previous works also in-\ndicates that nonorthogonal determinants can provide compact representations of the wave\nfunction. The code accompanying the paper was nicely set up and easy to use.\nReply: We thank the Reviewer for understanding the importance of our work.\nComment 1b — However, I have some general and specific remarks regarding the paper.\nGeneral remarks:\n1) The idea of using nonorthogonal determinants is not new, as also indicated by the authors\nsince they cite some of the literature. Except that the authors seem to have an e cient\nimplementation,Ican’tseethatthisapproacho erssomethingthattheothermethodsdonot\no er. If so, I think this needs to be adressed more clearly in the paper.\nReply: We agree with the Reviewer that Slater determinants have a fundamental role in\nQuantum Chemistry and, accordingly, many works have been dedicated to their study. Ob-\ntaining compact representations of the wave function with few Slater determinants is a long\nstanding problem. In particular, optimizing sums of non-orthogonal SDs has a long history\ndatingbackalmostacentury. Ournovelalgorithm,EIDOS,allowstoachievethislong-sought\ngoal in a systematic way.\nWe changed the introduction of the article to emphasize more that the novelty is not in the\nvariationalwavefunctionansatz,butinthemethodtooptimizeit(seeItem1oftheSummary\nof Changes below).\nAs the Reviewer correctly points out, given the centrality of the problem and its long his-\ntory, we discuss in the article classes of methods that have been proposed to optimize non-\northogonal SDs, and we point out their limitations.\nThe method that comes closest in goals to EIDOS is ResHF. However, we are not aware of\nanyworkthatmanagestofullyoptimizemorethanahandfulofnon-orthogonaldeterminants\natthesametime. Ontheotherhand,thereistheFEDapproach,optimizingonedeterminant\nat a time, while keeping the rest fixed. In Fig. 1 of this Reply we present our EIDOS results\nwith the alternative approaches of J. Chem. Phys. 139, 204102 (2013) (Fig. 1, cited as Ref.\n17inthecurrentversionofthemanuscript),whichisthestate-of-the-artforResHFandFED\nimplementations. Comparingthecurvesof“SD(UHF)cart”inFig.1ofthisReplyto“UHF”\n(black triangles) in Fig. 1 of J. Chem. Phys. 139, 204102 (2013), it can be seen that the\nwavefunction optimized with EIDOS has a significantly lower energy at the same number of\ndeterminants compared to the FED approach, or, in other words, that it can obtain a more\ncompact representation using fewer determinants. EIDOS can also systematically converge\nfor a much large number of determinants.\nNCOMMS-25-34613-T 2 of 18\nWe revised the introduction of the article to emphasize more that the key novelty lies in the\nway we optimize the determinants, departing from existing approaches (see Item 1).\nFigure1: Ground-stateenergyoftheN moleculeat1and1.5timestheequilibriumdistance\n2\nin the cc-pVDZ basis as a function of the number of non-orthogonal determinants, for both\ncartesian (cart) and spherical (sph) representation of the gaussian basis set, without any\nsymmetry projections. We compare to the spherical-basis FED results of Fig. 1 of J. Chem.\nPhys. 139, 204102 (2013), from where we have taken the geometries and coupled-cluster\nenergies, and from which we have extracted the black datapoints.\nComment 1c — 2) Although the boundary between dynamic and static correlation is a\nlittle bit blurry (especially along a dissociation curve) it is di cult to assess to which extent\nit captures static and/or dynamic correlation. In particular, in Fig 1 everything is compared\nto CCSD(T), indicating that it is dynamical correlation which is the target, since CCSD(T)\nis not suitable for static correlation. For some of the molecules in Fig 1, the results are closer\nto CCSD(T) than DMRG and FCI. However, the paper only talks about static correlation.\nReply: WeagreewiththeReviewerthatthereisnoprecisemathematicaldefinitionofstatic\nand dynamic correlation, and therefore distinguishing between them quantitatively is highly\nproblematic.\nNevertheless, at a qualitative level, using the Reviewer’s arguments on dynamic correlation,\nFig. 1 of the manuscript actually shows that the method can capture dynamic correlation.\nIndeeed, for molecules at equilibrium where single-reference methods such as CCSD(T) are\naccurate(whichcorrespondstoallmoleculesinFig. 1, withtheexceptionofperhapsC , N ,\n2 2\nO and BeO), we get excellent results.\nreference configurations, neither does the EIDOS optimization algorithm we use to optimize\nit. It is fully unbiased and captures as much correlation as possible, subject only to the\nconstraint on the number of determinants, regardless of whether the correlation is static or\nNCOMMS-25-34613-T 3 of 18\ndynamic.\nInordertospecificallytesttheabilitytooftheAnsatztocapturestaticcorrelation, itisbest\nto consider out of equilibrium geometries, as we do in Fig. 4 and in the section “Capturing\nstatic correlations.”, which is discussed in detail in the next point raised by the Reviewer.\nWe thank the Reviewer for stressing this reasonable interpretation of dynamic correlation,\nand we have added a discussion about it in the article (see Item 12). We have also decided\nto use as reference FCI and DMRG in Fig. 1 to avoid possible confusion about the ability of\nEIDOS to capture dynamic and static correlation.\nComment 1d — 3) With respect to static correlation, I don’t understand why CCSD(T)\n(as highlighted in the abstract and discussion) is used as a reference for comparison, It is\nwell known that low order truncations of CC are only suitable for capturing dynamical cor-\nrelation and is not suitable for bond breaking (static correlation). Hence, since “Capturing\nstatic correlations” is mentioned explicitly, using CCSD(T) as a reference is the wrong choice.\nCCSD(T) contain dynamical correlation, which is not suitable for the stretching of the bond\nwhereas DMRG only contain static correlation. Hence, for Fig 4, when the error wrt DMRG\nisalreadyatmilliHartreelevel(Nsd=768)thatisjusttheerrorinthestaticcorrelation,andit\ndoes not include the dynamical correlation error (which I suppose can be large for stretching\nsuch an electron rich bond).\nReply: We agree with the Reviewer that the bond breaking of the N molecule is relevant\n2\nto understand the role of static correlation. This is why in Fig. 4 we not only compare with\nCCSD(T) (more precisely, UCCSD(T)) as stated by the Reviewer. On the contrary, in the\ninsetofFig. 4,whenweanalyzethebond-breakingregime,whichistheonemostrelevantfor\nstaticcorrelation,weuseDMRGasreference,amethodthattheReviewercorrectlyidentifies\nas able to capture static correlation.\nWe do not fully agree with the Reviewer’s statement that “DMRG only contains static cor-\nrelation”. While we agree that DMRG might not be the most e cient or optimal method\nfor dynamic correlation, for large-enough bond dimension it must capture all correlation en-\nergy, including dynamic correlation, as it is a numerically-exact technique, just like EIDOS.\nBy pushing a DMRG calculation it is possible to get arbitrarily close to the true answer in\nprinciple, or at least verify that there is a problem in the convergence of the algorithm.\nIn the inset in Fig. 4 of the manuscript we compare with the DMRG results with a bond\ndimension of 4000 from Table III. of J. Chem. Phys. 121, 6110–6116 (2004) (Ref. 38 in the\ncurrent version of manuscript). These variational-energy results are consistently lower than\nMulti-reference CI and coupled-cluster calculations, which should be able to capture both\nstatic and dynamic correlation (see Fig. 4 of the same reference), and compatible with the\nhighest-order coupled cluster calculations.\nFollowing the Reviewer’s requests for clarifications, we updated the caption of Fig. 4 of\nthe article to emphasize the crucial role of the bond-breaking regime shown in the inset in\ndiscussing static correlation, and the comparison with DMRG (see Item 10).\nNCOMMS-25-34613-T 4 of 18\nComment 1e — 1) The random rotation of the orbitals: Is it correct to assume that there\nis no initial guess of the determinants, since the orbitals are randomly rotated? And since it\nis called a rotation, I suppose the transformation matrix in eq 12 is unitary. If so, this should\nbe specified clearly. Also, I cannot find it clearly explained why one does the random rotation\nof the orbitals, and the magnitude of such a rotation (given it is unitary, there is some in the\nrotations). Further, given that nonorthogonal orbitals are used, the metric of the start guess\norbitals is relevant. If the transformation matrix is unitary, the metric is will be conserved.\nReply: We thank the Reviewer for raising this question. The Reviewer correctly points out\nthat there is no initial guess of the determinants, and they are randomly initialized with a\ngaussian distribution. The use of the word “rotation” in Algorithm 1 was not optimal since,\nas we write in the main text, one can use any determinant-one matrix, and we thank the\nReviewerforspottingthistypo. WealsochangedtheletterofthematricesinEq. 12toavoid\nconfusion, rewrote the entire paragraph, and updated Algorithm I to make this step more\nclear (see Item 2).\nComment 1f — 2) The molecular orbitals: I am not sure why the authors use quotation\nmarkswhenwritingmolecularorbitals. AsfarasIcansee,theyarejustnormalnonorthogonal\nmolecularorbitals(unlessIammistaken?),andthereisnoneedfortheuseofquotationmarks.\nAlso, this is not important at all: but MO stands for “molecular orbitals”, so there is no need\nto talk about “MO orbitals”, one can merely use the term “MOs”.\nReply: We thank the Reviewer for raising this point. We use the term molecular orbitals\nin quotation marks to refer to the vectors\n„(I)\nwhich, as noted by the Reviewer, represent\ni\nlinear combinations of the atomic orbitals, mutually non-orthogonal both within the same\ndeterminant and for di erent determinants. We changed the wording throughout the main\ntext to refer to them simply as orbitals (see Item 11).\nComment 1g — 3) On page 2 it is written “EIDOS exactly optimizes the variational\nenergy..”. I don’t understand what is meant by exactly. Is it that the algorithm will only\nconverge to minima (as opposed to saddle points)? And if so, how is this ensured?\nReply: We report here, for completeness, the full sentence cited by the Reviewer: “EIDOS\nexactly optimizes the variational energy when viewed as a function of a string of orbitals,\nobtained by selecting one orbital from each determinant.” What we meant here is that, if\nwe consider all but one orbital per determinant to be fixed, then we can find the optimal\nparameters for the remaining orbital of each determinant exactly by solving a generalized\neigenvalue problem.\nConcerning global convergence: while we do not have any of theoretical guarantee apart\nfrom the fact that the variational energy does not increase after an optimization step, we\nempirically found that, by iteratively optimizing di erent orbitals (and freezing all the rest)\nwe typically converge to states which are very close to the ground state. We expect the\nrandom SL transformation of the orbitals which we introduced to the optimization process\nNCOMMS-25-34613-T 5 of 18\n(which includes random orbital rotations) to help move the system away from any saddle\npoints and eventually lead it towards a minimum. However, the presence of local minima in\naddition to the global minimum remains a possibility. With regards to this, we can add that\nwhilestudyingtheN moleculeoutofequilibriumwherewewouldobservespincontamination\n2\nfor small number of determinants. Please refer to our answer to Comment 2e of Reviewer 2\nfor more details.\nWe rewrote the aforementioned article sentence to make it more clear (see Item 6).\nComment 1h — 4) In the next sentence on page 2 it says “EIDOS relies on the direct\nparametrization of the molecular orbital coe cients- as opposed to the popular Thouless\nparametrization of the SDs”. I think this needs some clarfication. Acting with a Thouless\ntransformation(e.g. exponentiallyparametrizedunitaryorbitalrotationoperator)onaSlater\ndeterminant yields a parametrization of the orbitals coe cients directly (equal to that of eq\n12 in manscript).\nReply: In terms of variational expressivity, we agree that the Thouless representation is\nequivalenttothedirectorbitalparametrizationusedinthearticle(oneisthematrixlogarith-\nm/exponentialoftheother). However,thecrucialpointforwhatconcernstheoptimizationis\nthat the Thouless parametrization does not lead to a multilinear wave function in its param-\neters. This fundamental di erence, while not formally a ecting the variational expressivity,\nrenderstheexactsingle-orbitaloptimizationtechniqueweproposeinthearticle,EIDOS,sim-\nply impossible within the Thouless parametrization. As the Thouless parametrization has a\npredominant role in the literature, this explains why, despite a long history, only inexact first\norder and second order algorithmic optimization strategies were considered before our work.\nThis also gives some more historical and technical context for the question of Comment 1i.\nWe have added a comment highlighting the di erence w.r.t. the Thouless representation in\nthe main text of the article (see Item 7).\nComment 1i — 5) Regarding Thouless parametrization of the states, there is also a state-\nmentonpage2regardingusinggradient-basedmethodsanddi cultiestoconverge. Ithinkit\nthe second order approaches of wave function optimization (also using Thouless parametriza-\ntions) could be mentioned.\nReply: We already cite several works employing second-order and quasi-second-order op-\ntimization methods in the Introduction. We would be happy to consider additional, more\nspecific suggestions, as the relevant literature is quite extensive. See also Comment 1h for\nfurther context on the Thouless parametrization.\nComment 1j — 6) Under “Methods” the authors write “(assumed here, for simplicity,\nto be orthonormal for the usual anticommutation relations to hold)” It is unclear why this\nassumptionisneeded. IfitisforthepurposeofrearrangingcreationoperatorsenteringtheSD\ndefinitions, this assumption is not needed. Creation. operators always anticommute among\nNCOMMS-25-34613-T 6 of 18\nthemselves.\nReply: The Reviewer rightly points out that the orthonormality assumption is not needed\nin general. The anticommutator of general non-orthogonal creation operators is given by\n{cˆi,cˆ†j} = Sij\nwhere Sij is the overlap matrix of the orbitals, instead of ”ij for orthogonal\noneswhichwouldcomplicatetheequations. Weremarkthatthemoregeneralnon-orthogonal\nequations can be recovered by mapping to an auxiliary ortogonal basis, inserting S ≠ 1 2 into\ntheequationsthroughout. Howeverwedecidednottowriteitdowninthismoregeneralform\nforclarityofexposition. Alternatively, andequivalently, onecanalsodefineabi-orthonormal\n(dual) basis associated with the original non-orthogonal orbitals, together with the corre-\nsponding set of dual creation and annihilation operators.\nResponse to Reviewer 2\nComment 2a — The authors introduce their EIDOS method, a (as far as I know) novel\nalgorithm for variational optimization of non-orthogonal Slater determinant expansions, vari-\nationally optimizing both the orbitals and CI amplitudes simultaneously in a greedy quasi-\nrandomizedapproach. Theydemonstratethatusingafewhundredsuchfullyoptimizeddeter-\nminants yields ground-state energies for small molecules with accuracy rivalling or surpassing\nstandard methods (e.g., CCSD(T)), at a polynomial computational cost (for a fixed determi-\nnantnumber). Thehistoryofnon-orthogonalCIexpansionsislong,andtherehavebeenmany\napproaches suggested previously. The high accuracy of these expansions for small molecules\nis well known, although they are largely niche approaches due to the lack of size extensivity\nand (exponential?) scaling with number of determinants.\nThe main success of the work to be celebrated is the development of the algorithm for a full\nquasi-randomized variational optimization of the orbitals of the determinants. It is a neat\napproach which linearizes the determinants based on random rotations of the variations of a\nsingle degree of freedom. This is a significant achievement, and the ideas of this optimization\napproach could find broader impact.\nReply: We thank the Reviewer for recognizing the achievements of our article for this long-\nstanding problem.\nComment 2b — However, the results are only presented up to a maximum of 20 electrons,\nand while accurate, does not currently enable calculations which cannot be achieved with\nalternate state of the art methods (e.g. DMRG or selected-CI). This is a high bar for the\ninitial introduction of a new algorithm, and would still be happy to see it published as a\nnovel algorithm, but with some more results detailing the scope of systems which the method\ncould treat, and a wider discussion about the scaling of the algorithm with system size. In\nparticular, while the lack of size consistency is noted, this major limitation (and the related\nlack of size extensivity) is largely glossed over. I would suggest this needs further scrutiny for\nan more honest appraisal of the potential of the method and its prospects. While the SI has a\nNCOMMS-25-34613-T 7 of 18\nplot of size consistency for dissociated N2, more instructive would be the scaling of the error\nwith respect to system size (i.e. characterizing the size extensivity error) – do we expect to\nrequire an exponential number of determinants for fixed relative error? A formal exponential\nnumber of determinants is not a disqualifying characteristic either (see e.g. DMRG and sCI)\n– if this has a small prefactor/exponent, then this can still be a useful method, and at least\nclarify the domain of applicability of this tool.\nI would suggest e.g. a series of increasing length chain of Hydrogen atoms at a fixed inter-\natomicspacing,andtoobservetheerrorcomparedtoDMRGforfixeddeterminantnumberas\nadditional plots in the main text (and at changing bond lengths). Presumably the correlation\nenergywillconvergetozero,butatwhatrate? HowdoesthiscomparetoCISDorCCSD(the\nlatterbeingsizeextensive)? Howdoesthenumberofdeterminantsneedtoscaletocounteract\nthis? Thiswouldbemoreinstructivethanthebondorderorbasissetplots,thelatterofwhich\nin particular is not expected to scale exponentially. This plot may also clarify the potential of\ntheapproachtoextendtoresultsbeyondwhatispossiblewithDMRGorothersystematically\nimprovable methods, even if not possible with the current implementation.\nReply: We thank the Reviewer for suggesting this comparison, which we have performed\nin detail. In Fig. 2 of this Reply and Fig. 3 of this Reply we study Hydrogen chains in the\nSTO-6G basis, with up to 102 electrons.\nNCOMMS-25-34613-T 8 of 18\nFigure 2: Energy per atom for uniform hydrogen chains H in the STO-6G basis set for\nN\ndi erent values of the atomic spacing . Shown is the energy of di erent numbers of non-\nR\northogonal determinants n compared to DMRG.\nD\nNCOMMS-25-34613-T 9 of 18\nFigure3: Numberofdeterminantsnecessarytocaptureacertainpercentageofthecorrelation\nenergy as a function of system size. The system studied are uniform hydrogen chains in the\nSTO-6Gbasisforfourdi erentvaluesofinter-atomicspacing . pointsontheplotgenerated\nR\nthorugh fitting and interpolating the number of detrminants needed from the data provided\nin Fig. 2 of this Reply.\nThis detailed analysis shows that the Ansatz is not size extensive, as expected. While it is\nnotpossibletodeterminetheexactasymptoticnumberofdeterminantsforafixedcorrelation\nenergy relative error, the data is not inconsistent with an exponential increase, even if there\nappears e ects of saturation at large system sizes. This is certainly true at large separation\nbetween atoms, where the wave function is a product state of single-atom wave function,\nwhichcreatesanexponentialnumberofdeterminantsinordertomaintainafixedcorrelation\nenergy. Therefore, we strongly suspect that the lack of size extensivity has the same origin\nas the lack of size consistency, as expected. This was the motivation to only discuss size\nconsistency in the previous version of the manuscript.\nWehaveaddedthismaterialtothenewSectionHoftheSupplementaryInformation(seeItem3\nof the Summary of Changes below).\nNCOMMS-25-34613-T 10 of 18\nComment 2c — Fig. 1shouldcompareresultstoDMRGratherthanCCSD(T).Allsystems\npresented should be able to be solved essentially exactly with DMRG, and so would make a\nfar better reference.\nReply: We thank the Reviewer for this suggestion, which we have implemented. We ran\nadditionalDMRGcalculationsfortheremainingmoleculesandupdatedthefigureaccordingly\n(see Item 4).\nComment 2d — Only a couple of sentences are given in the main text (not in the methods\nsection) to describe the algorithm, and this is rather opaque. The methods section is good\nfor describing the technical details, but given that the algorithm itself is the main novelty of\nthe approach and major scientific advance, I would hope for a little more in the main text to\ndescribe it qualitatively.\nReply: We thank the Reviewer for the feedback. We have expanded the explanation of the\noptimization in the Results section of the Article (see Item 8).\nComment 2e — Furthermore, it would be good to have more insights into the convergence\nof the algorithm, especially in tricky cases. For instance, this is a quasi-randomized greedy\noptimizationforanon-convexoptimizationproblem(characterizedas‘partial-but-exact’inthe\nmain text, which I found particularly uninformative!), so presumably there are cases where\nit finds local minima, exhibits a slow-down in convergence, or is found particularly sensitive\nto the initial conditions (perhaps e.g. where there are low-lying excited states?). It would be\ngoodtoincludesomefurtherinsightsintothenumericalperformanceofthealgorithmbeyond\nFig. B1, also commenting on the rate of convergence for di erent cases and starting point\n(e.g. keep the same random rotations, but reorder the initial orbitals?). How does this also\nvary with system size? Where is this algorithm expected to struggle to optimize to the true\nsolutionscomparedtoalternateiterativemethods–arethereanyinsightshere? -Itwouldbe\ngood to compare the approach to modern approaches to (generalized) non-orthogonal CI, or\nfragment-localized optimization algorithms.\nReply: In general, we have found the algorithm to converge quite consistently, assuming\nall the necessary steps have been taken to ensure numerical stability (see Section A in the\nSupplementary Information). One case where we did observe local minima was for the N\n2\nmolecule distances larger than the equilibrium when the bond was breaking, but only when\nusing a relatively small number of determinants. Using a larger number of determinants\nusually eliminates this problem.\nMore specifically, in Fig. 4 of this Reply we include a plot of the optimization at a distance\nof 4.2 a.u. with 64 determinants, 92 simulations from scratch, with di erent seeds, and the\nsame using 5 complete optimizations with 512 determinants. As the number of determinants\nis increased they all seem to converge to similar states.\nNCOMMS-25-34613-T 11 of 18\nFigure4: OptimizationofN atbonddistanceof4.2a.u. inthecc-pVDZbasis. Completeop\n2\ntimization from scratch starting from randomly initialized determinants with di erent seeds.\nLeft column: absolute energy error w.r.t. DMRG as a function of the number of iterative\nsteps. Right column: Energy error and total spin of the final wavefunction of each optimiza\ntion. Top row: 92 simulations with 64 determinants. Bottom row: 5 simulations with 512\ndeterminants.\nCompared to alternative iterative optimization approaches for the same Ansatz, we are not\naware of any particular disadvantages of the approach presented here. In fact, it optimizes\nbetter than comparable approaches such as FED for this exact system, see our reply to\nComment 1b of Reviewer #1.\nRegarding the fragment-localized approach, a procedure where orbitals are localized on frag\nments of molecules, it seems to us that this approach is targeting very large large molecules,\nwhich are not the scope of the method introduced in our article.\nComment2f— -Thereareafewminortypographicalerrorsandinconsistentnotation(e.g.,\nmissing or mismatched parentheses and indices) remain in the text—careful proofreading is\nNCOMMS-25-34613-T 12 of 18\nadvised.\nReply: Wehavefixedseveraltyposthroughoutthetext(seeItem9),andwewouldbehappy\nto know if we missed some.\nComment 2g — Overall,themanuscriptintroducesaninterestingandinnovativealgorithm\nforoptimizingNOCIexpansions. However,itdoesnotclearlydemonstrate,eithernumerically\nor a theoretical path, towards how the approach can be extended to larger systems or achieve\nhigher accuracy (either in absolute terms or as a function of wall time) compared to existing\nalternatives. IwouldbepleasedtoseethisworkpublishedinNatureCommunicationsifthese\npoints can be addressed. Specifically, a more thorough and candid assessment of the algo\nrithm’s scaling with system size is needed, both in terms of computational cost and accuracy\nwith respect to number of determinants. This will help clarify the potential niche and prac\ntical applicability of this method, and justify its relevance for publication in a broad-interest\nscientific journal.\nReply: We thank the Reviewer for the suggestions and we are confident that the additional\nextensive results prepared for the Reply clarify the points raised.\nResponse to Reviewer 3\nComment 3a — The concept of releasing the orthogonality restrictions to obtain more\ncompactwavefunctionsisaveryattractiveoneandhasreceivedquitesomeattentionoverthe\nlastdecade. Intheherepresentedmanuscript,amethodisdescribedthate cientlyoptimizes\nawavefunctionexpressedasalinearcombinationofnon-orthogonalSlaterdeterminants. The\nmethodhasthepotentialtobeground-breakingandthearticlecanbecomealandmarkpaper\nin the field.\nReply: We thank the Reviewer for recognizing the significance of our work.\nComment 3b — There are a few things that need to be looked at to put the paper in the\ncontext of previous work (point 1); to discuss the stability of the algorithm when larger, more\nrealistic systems are targeted (point 2); to get a feeling of how important non-orthogonality is\n(point 3); to address the importance of full orbital relaxation (point 4) and to clarify the way\nin which the orbitals are picked for optimization (point 5).\n1. A comparison with the work of Olsen is in place (JCP 143, 114102). In this study (and\nsome follow-up papers) a procedure is described in which the wave function is optimized with\nrespect to both the CI and orbital coe cients in which the configuration state functions are\nexpressed in non-orthogonal orbitals.\nNCOMMS-25-34613-T 13 of 18\nReply: The cited article by J. Olsen describes a CASSCF-like MCSCF procedure where the\norbitals in the active space are allowed to be non-orthogonal, and provides the equations for\ngradients and Hessian in the Thouless representation, as well as an e cient implementation\nof this based on an iterative solver and direct hessian-matrix-vector products. What we\npropose here is also based on non-orthogonal orbitals; however, it is more flexible, in the\nsense that every determinant has its own independent set of independent orbitals, which are\nall optimized, i.e. no orbitals are shared between determinants, and we optimize all orbitals,\nnot only those of an active space.\nRegarding direct hessian-matrix-vector products, this is what is required for the use of it\nerative solvers as we proposed in the discussion section, and certainly represents a future\ndirection worthwhile of further explorations.\nWeaddedacitationofthearticleofJ.Olsentotheintroduction(seeItem13oftheSummary\nof Changes below).\nComment 3c — 2. The use of iterative solvers seems to be an unavoidable step to surpass\nthe memory bottleneck. Although not further worked out in the manuscript, this can be a\ncritical step in the further development of the method and make it capable to use for other\nthan the model molecules discussed here. With the increased number of degrees of freedoms\nduetothelossoforthogonalityandthepossiblelineardependencies,theuseofiterativesolvers\nmaydegradatethestabilityoftheoptimizationprocess. Moreinformationontheprospectsof\niterative solvers could be very helpful for other researchers in the field to judge the suitability\nof further developing this promising method.\nReply: In our brief initial testing we found it di cult to converge with iterative solvers.\nThis is an important technical point for further developments that will necessitate additional\ninvestigation in future work.\nAn alternative straightforward approach for scaling the method to larger number of deter\nminants, as we have mentioned in the discussion section, is to optimize an orbital of a only\na subset of the determinants at every step. This functionality is already present in the code\ncapsule submitted with the initial version of the article.\nWe have updated Fig. 1 of the Article to include results obtained with the partial- optimiza\ntionapproach, andwehaveaddedsometechnicaldetailsonthisextensionintheSupplemen\ntary Information (see Item 5). This has allowed to converge with up to 2048 determinants\nand significantly improve the variational energy of Fig. 1.\nComment 3d — 3. What is the degree of non-orthogonality among the di erent deter\nminants? Maximum overlap, percentage of (nearly) orthogonal determinants, to name a few\nindicators. An overview of some numerics on this point is more than welcome.\nReply: To answer this question, in the following we present the simulation results of the C\n2\nmolecule in the cc-pVDZ basis with 768 non-orthogonal determinants. For completeness, we\ninclude in Fig. 5 of this Reply a plot of the norms of the determinants, sorted in decreasing\nNCOMMS-25-34613-T 14 of 18\norder, showing that all determinants contribute to the wavefunction.\nFigure 5: Norm   =     of the non-orthogonal determinants in the optimized\ni i i\nÎ Î È | Í\nwavefunction of the C molecule in the cc-pvdz basis for 768 determinants. Ordered w.r.t.\n2 \ndecreasing norm.\nIn Fig. 6 of this Reply we show the histogram of the distribution of the maximum fidelity of\nevery determinant with any other one.\nFigure6: Histogramofthemaximumfidelityofeachofthe768determinantswithanyother\none. Values range from 0.018 to 0.987.\nNCOMMS-25-34613-T 15 of 18\nAs can be seen from Fig. 6, for every determinant there are other ones which are non\northogonal to it. The maximal fidelity of a determinant with all the others is always larger\nthan 0.018, and not smaller than 0.1 for all but 2% of the determinants.\nFinally, in Fig. 7 of this Reply we plot of the fraction of pairs of determinants with a fidelity\nlower than a certain cuto , showing how many are close to orthogonal to each other.\nFigure 7: Fraction of entries of the pairwise fidelity matrix F ij = È  i |È |     i i Í |   È   j j Í| | 2  jÍ below cuto .\n  and   goes over all elements of the set of 768 optimized non-orthogonal determinants for\ni j\nthe C2 molecule in the cc-pVDZ basis.\nThis shows, for example, that around half of all pairs of determinants have a fidelity of less\nthan 0.01, and 10% of less than 10 4.\n≠\nComment3e— 4. Themethodminimizestheaverageenergyofthewavefunctionexpressed\nas a linear combination of non-orthogonal Slater Determinants. This is di erent than the\napproaches that adopt a \"di erent orbitals for di erent electronic configurations\" strategy.\nThelatterincludesfullorbitalrelaxationforeachelectronicconfiguration(configurationstate\nfunction, Slater determinant) in the wave function, while the here presented EIDOS variant\nonly partially recovers orbital relaxation. This is especially relevant for core-hole states or for\nelectronicstateswherechargetransferexcitationscontributeconsiderablytothewavefunction.\nIn such cases a much larger number of SD will be needed to obtain a balanced description\nof the di erent configurations. Is is possible to design a strategy to perform independent\noptimizations of the most important electronic configurations and subsequently account for\nthe interaction between these fully relaxed (non-orthogonal) states?\nNCOMMS-25-34613-T 16 of 18\nReply:\nAs long as ground-state problems are concerned (the scope of the present manuscript), our\nalgorithmensures,inastochasticsense,thateveryorbitalisfullyoptimizedintheend. Thus,\nthemethodcorrespondsto‘di erentorbitalsforeachSlaterdeterminant’,withnoconstraints\nin the optimization problem; every orbital is fully relaxed. However, the point raised here\nby the Reviewer on the possibility of using the method to tackle excited states appearing,\ne.g., in reactions and spectroscopical experiments, is a valid future direction. In general, we\nwouldliketoremarkthattheonlypropertywerequireisthatthequantityminimizedcanbe\nexpressed as a ratio of two quadratic forms, which is the case for the energy of a single state,\nfidelity. For example, a state-averaging approach would not be possible due to the di erent\ndenominators in the sum of the average energies. Regarding the final question, it should\nalways be possible to optimize two di erent states (e.g. two di erent geometries) separately,\nand then do a NOCI calculation of the union of all individual determinants.\nComment 3f — 5. What is the reason to randomly select the orbital that is singled out\nfor optimization? Does this imply special advantages compared to a simple linear selection\nprocedure? How can we be sure that all orbitals in the SD are at least optimized once? Isn’t\nthis procedure prone to get stuck in a local minimum? Furthermore, it is not at all clear how\nsteep will be the increase in the number of iterations in this loop when bigger systems (with\nmuch more occupied orbitals) are studied.\nReply:\nRegarding the choice of optimization strategy, after doing additional testing on a simple LiH\nmolecule we could not find major di erences in convergence between simple sweeps, uniform\nslelection and e.g. a more advanced strategy based on a unitary rotation matrix, at least for\nthis simple case. However, we remark that, in some cases, it was necessary to orthogonalize\nthe orbitals to ensure numerical stability, which itself is already a SL transformation, and\nin this case simple sweeps are not possible to realize. Given enough iterations are done, a\nrandom selection will eventually optimize all of them at least once, this is an instance of the\ncouponcollectorsproblem,andsoO(nlogn)stepsshouldsu cefornorbitals. However,using\nrandom unitary or SL transformations we expect this to be achieved faster, since at every\nstep a linear combination of all the orbitals of the previous step are optimized. Furthermore,\nwhen using UHF-type determinants it is necessary to choose between optimizing the spin-up\nand spin-down parts of each determinant. At most, we would expect this could potentially\nget stuck if we do not mix up and down optimization between di erent determinants.\nIn general, we have found the algorithm to be quite robust. The initialization of random\ndeterminantsalsodoesnotposeproblems,astheyaregenerallyoptimizedalreadyafterafew\nsteps, as we show below.\nForthenitrogenmoleculeinthebond-breakingregimewedidobservethealgorithmconverg\ning to di erent local minima when using small numbers of determinants, please refer to our\nreply to comment 2e of Reviewer 2 in which we explain this in more detail.\nTo give more insight about the convergence for larger systems we show the optimization of\nNCOMMS-25-34613-T 17 of 18\nthe additional results for the hydrogen chain discussed in the answer to comment 2b, where\nwedoublethenumberofdeterminantsevery10000optimizationsteps. Plottedistherelative\nerror in the correlation energy compared to DMRG.\nFigure 8: Convergence of the relative correlation energy error w.r.t. DMRG as a function\nof the number of steps. for di erent system sizes and increasing number of determinants\n(initialized from the previous smaller simulation, adding randomly initalized determinants).\nAs can be seen from the plots, here we do not observe a significant di erence in the amount\nof iterations needed to converge for a fixed number of determinants, upon increasing the size\nof the system.\nWe rewrote the relevant section in the article for the optimization to make it more clear (see\nItem 2).\nNCOMMS-25-34613-T 18 of 18\nSummary of changes\n1. We have emphasized more in the introduction that the novelty is the optimization strategy,\nnot the Ansatz which has been widely used before.\n2. We have rewritten the section describing the random transformation of the orbitals in the\noptimization.\n3. We have added results for the 1D hydrogen chain to the Supplementary Information.\n4. We have changed Fig. 1, where with additional DMRG calculations we can now compare to\nFCI and DMRG instead of coupled cluster on the y axis.\n5. We have added to Fig. 1 results obtained with a partial optimization of the determinants,\nwhich is needed to obtain results with 2048 determinants.\n6. We have adapted the wording “exact” in the explanation of the method to render it more\nprecise.\n7. We have added a comment clarifying that in the Thouless representation the ansatz is not\nmultilinear.\n8. We have expanded explanation of the method when it is introduced in the Results section.\n9. We have fixed a few typos throughout the article, as can be seen in more detail from the\nincluded latexdi  file.\n10. We have adapted the caption of Fig. 4 to emphasize more the importance of the inset.\n11. We have replaced occurrences of MO orbitals with orbitals.\n12. We have added a sentence discussing dynamic correlation in the discussion of Fig. 1.\n13. We have added a citation for JCP 143, 114102 in the introduction.\nPlease see the included latexdi  file for the exact changes.\nNCOMMS-25-34613-T 1 of 4\nResponse to the Reviewer 1\nComment 1a — The authors have addressed most comments provided in a satisfactory\nmanner, except for two points. Which should be addressed.\nReply: We thank the reviewer for his considerations. We address the remaining two points\nbelow.\nComment 1b — 1)Theautorsstate(Item1ofSummaryofchanges)theyhaveemphasized\nmore in the introduction that the novelty lies in the optimization strategy rather than the\nansatz. However, I cannot find the changes which does this. I only find minor updates of\nsentence structures, moving paragraphs etc.\nReply: We acknowledge that the updates that we have performed in the introduction were\nminor as we believe that the last paragraph is clear about the fact that the novelty is in the\noptimization strategy rather than the ansatz. In particular, we are of the opinion that the\nsecond phrase of the paragraph is explicit about this fact: “These results are achieved with\na computational optimization technique that we introduce in this article, which we call the\nExact Iterative Determinant-Orbital Solver (EIDOS).”\nComment 1c — 2) In comment 1j, I ask about if the assumption of antcommutativity is\nneeded, since I understood it to be related to the ordering of creation operators in eq 3 and 4.\nCreationoperatorsalwaysanticommuteamongthemselves,whetherthebasisisorthogonalor\nnot.\nHowever, the author reply discusses an anticommutator of an creation and an annihilation\noperator, which of course relates to the orbital overlap.\nThis raises new questions for me. Since they are discussing anticommutation of creation and\nannihilationoperators,Itakeitthismeansthattheyapproximatetheanticommutativitywhen\nexpressing the expectation values of the Hamiltonian (eq 7)? If this is the case, then I am\nwonderinghowupperboundisensured, ifthematrixelementsarenotcorrect(i.e., ifderiving\nthe matrix elements uses an orthogonality which is not there).\nHowever, if you are indeed talking about the ordering of the creation operators, this is not an\nassumption but a fact, both for orthogonal and nonorthogonal bases.\nReply: The last paragraph, which we highlight, clarifies the misunderstanding. Indeed, we\nhave not assumed that the canonical anticommutation relations (see Eq. (1) of this Reply)\nhold; instead, we have assumed that the basis set is orthonormal, which implies that the\ncanonical anticommutation relations (see Eq. (1) of this Reply) hold. For completeness, we\nshall reproduce the full phrase in the article: “For simplicity, we assume the basis to be\northonormal; this implies, in particular, that the usual anticommutation relations hold”.\nWe have updated this sentence in the manuscript to include the explicit anticommutation\nNCOMMS-25-34613-T 2 of 4\nrelations equations from Eq. (1) of this Reply (see below), to avoid any misunderstandings.\nIn an e!ort to be as complete as possible in our response, as an optional reading, we pro\nvide below an answer to the first part of the comment 1c, even if, as discussed above, it\nis not necessary as the emphasized paragraph of comment 1c clarifies the triviality of the\nmisunderstanding.\nComment 1c - first paragraph — 2) In comment 1j, I ask about if the assumption of\nantcommutativity is needed, since I understood it to be related to the ordering of creation\noperators in eq 3 and 4. Creation operators always anticommute among themselves, whether\nthe basis is orthogonal or not.\nReply: We want toclarifythat inthis work we use standardsecondquantization formalism.\nInparticular,tobefullyprecise,thefermionicoperatorsdonotanticommute,buttheysatisfy\nthe canonical anticommutation relations given by\n{cˆµ,cˆω} = {cˆ†µ ,cˆ†ω} =0,\n{cˆµ,cˆ†ω} =ωµ,ω. (1)\nComment 1c - second paragraph — However, the author reply discusses an anticom\nmutator of an creation and an annihilation operator, which of course relates to the orbital\noverlap.\nReply: In our article we use canonical anticommutation relations, see Eq. (1) of this Reply.\nInthepreviousreplytotheReferee,wewerediscussinganextensionoftheformalism,which,\nfor instance, allows to directly deal with a non-orthonormal basis set of orbitals. This would\nhave a computational advantage for a very large basis, as this would avoid to transform the\nHamiltonian into the orthogonal basis. This extension of the formalism is not needed to\nunderstand the present article.\nComment 1c - third paragraph — This raises new questions for me. Since they are\ndiscussing anticommutation of creation and annihilation operators, I take it this means that\ntheyapproximatetheanticommutativitywhenexpressingtheexpectationvaluesoftheHamil\ntonian (eq 7)? If this is the case, then I am wondering how upper bound is ensured, if the\nmatrix elements are not correct (i.e., if deriving the matrix elements uses an orthogonality\nwhich is not there).\nReply: No approximations are made in this article. Every expression and datapoint is\nnumerically exact up to the finite numerical precision permitted by the hardware (we use\ndouble precision for all our calculations). In particular, we use Eq. (1) of this reply, without\nany approximations. The energies that we provide are strictly variational.\nNCOMMS-25-34613-T 3 of 4\nResponse to the Reviewer 2\nComment 2a — Iamhappythatthequestionshavebeenlargelyaddressedintheresponse\nand manuscript modifications, and that I can recommend publication.\nReply: We thank the reviewer for evaluating our response and for recommending the article\nfor publication.\nComment 2b — The only further request (that was requested in the original response),\nis that CCSD (size extensive) and CISD (non size extensive) results are added to the new\nplots in H2. This will also show the value added from the non-orthogonality in compactifying\nthe determinant expansion by comparison to CISD (as well whether the non-extensive energy\ncomponent scales the same way). Furthermore, the comparison to CCSD will quantify the\nmagnitude of the size extensivity error, and when this dominates over the e!ective state.\nReply:\nWe thank the reviewer for this suggestion, which we discuss in the following.\nInFig.1weshowtherelevantplotofthecurrentversionofthearticle,whereweaddedCISD\nand CCSD.\nTheadditionallinesinFig.1ofthisReplyarevisuallyveryhardtodistinguish. Accordingly,\nwe have decided not to update the old Figure. Instead, in the supplementary material of\nthe article we include a zoomed version of the plot for a particular value of , see Fig. 2 of\nR\nthis Reply. We also scale the number of determinants linearly with the system size N, as\nwe believe provides a more appropriate basis for comparison with CISD (whose number of\northogonal determinants scales like N4), and makes the figure less crowded.\nAs we write in the updated version of the article: “We note that EIDOS provides the best\nvariationalenergies,withthesoleexceptionofDMRG,whichisknowntobeoptimalforthese\nlinearchains. ItoutperformsCISDwhileusingordersofmagnitudefewerdeterminants,which\nhighlights the compactness of the representation made possible by the use of non-orthogonal\ndeterminants.”\nResponse to the Reviewer 3\nComment 3 — The authors did a fine job revising the manuscript, provided a lot of extra\ninformation and clarified most of the points raised by the reviewers in a convincing way. This\narticle is a nice starting point for further elaboration in the near future. I will certainly be\nfollowingthedevelopmentofthemethodandverycurioustoseehowitevolves(excitedstates,\ngradients, properties, etc.)\nReply: We would like to thank the reviewer for appreciating the article and the potential of\nthe method.\nNCOMMS-25-34613-T 4 of 4\nFigure 1: Energy per atom for uniform hydrogen chains H in the STO-6G basis set for di!erent\nN\nvalues of the atomic spacing . Shown is the energy of di!erent numbers of non-orthogonal deter\nR\nminants n compared to DMRG, UCISD and UCCSD.\nD\nFigure 2: Energy per atom for uniform hydrogen chains H in the STO-6G basis set with a fixed\nN\nthe atomic spacing of = 1.8. We scale the number of non-orthogonal determinants N linearly\nD\nR\nwith system size N and compare the resulting energy against DMRG, UCISD, and UCCSD. For\nUCISDwealsoreportthe(rounded)numberoforthogonaldeterminantsintheexpansionasalabel\nnext to the datapoints."}]

Decision letter