A Tight Lower Bound to the Ground-State Energy | Journal of

Jun 3, 2019 - ... ACS Catalysis, ACS Central Science, ACS Chemical Biology, ACS ..... This lower bound expression is not yet practical since one needs...
0 downloads 0 Views 778KB Size
Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. 2019, 15, 4079−4087

A Tight Lower Bound to the Ground-State Energy Eli Pollak*

Downloaded via NOTTINGHAM TRENT UNIV on August 14, 2019 at 01:44:03 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Chemical and Biological Physics Department, Weizmann Institute of Science, 76100, Rehovot, Israel ABSTRACT: Ninety years ago Temple (Proc. R. Soc. (London) 1928, A119, 276) derived a lower bound for the ground-state energy. The bound was tested and invariably found to be poor as compared to the upper bound obtained through the Rayleigh Ritz procedure due to the fact that it is based also on the second moment of the Hamiltonian. In this paper we (a) improve upon Temple’s lower bound estimate for the overlap squared of the true ground-state wave function with the approximate one and (b) describe in detail and generalize our recent improvement on the Temple lower bound based on utilization of higher-order basis functions derived by the Arnoldi algorithm. Both improvements combined lead to a lower bound on the ground-state energy whose accuracy is better than that of the Temple lower bound. This is exemplified by considering the ground-state energy of a quartic potential where one finds that the improvements lead to a lower bound whose quality is comparable to that of the upper bound. The applicability of the method to atoms and molecules is discussed.

I. INTRODUCTION Obtaining upper bounds to the ground-state energy is a welldefined and straightforward task. The Rayleigh Ritz variational principle assures us that if we use any finite basis set to represent the Hamiltonian operator Ĥ (or any other hermitian operator for that matter), then diagonalization of the matrix will provide an upper bound to the ground-state energy (or equivalently the lowest eigenvalue). As one increases the dimensionality of the basis set, so will the resulting upper bound come ever closer to the true ground-state energy. This procedure has been used to obtain the ground-state energy of the helium atom with an accuracy of 41 significant figures.1 When applying the Temple lower bound2 to the same computation, the best accuracy was for only 32 significant figures. This huge difference between lower and upper bounds has troubled the community of lower bound practitioners,3−6 yet no meaningful method was developed thus far which could compete with the quality of the upper bound. Why is not the Temple class of lower bounds tight? It and all related formulas such as the Weinstein lower bound,7,8 rely among others on the variance of the energy 2 (σΦ = ⟨Φ|Ĥ |Φ⟩ − ⟨Φ|Ĥ |Φ⟩2 ) obtained from the trial ground-state wave function |Φ⟩. For the exact ground-state this variance vanishes. For a trial wave function it is finite and becomes smaller as the quality of the trial wave function improves. Yet, as noted by Nakashima and Nakatsuji,1 when evaluating the ground-state energy as ⟨Φ|Ĥ |Φ⟩, the integrand Φ(x) [HΦ̂(x)] is not necessarily positive definite so that positive and negative inaccuracies may cancel each other out. When considering the variance one is integrating over the positive definite function [HΦ̂(x)]2 and all the inaccuracies in the integrand add up. It is well understood that minimizing ⟨Φ| Ĥ 2|Φ⟩ is not an efficient way of computing energy eigenvalues. The need for good lower bounds is almost obvious. If one could obtain both upper and lower bounds of the same quality © 2019 American Chemical Society

one would have a reliable way of estimating the accuracy of a given computation, without having to converge it. This for example would be useful for renormalization group computations9−11 where one of the central questions always is whether the renormalization procedure did not lose anything of vital importance along the way. But even for more ”standard” problems, knowledge of a tight lower bound would not only provide an estimate of the accuracy but could also be used to improve upon previous estimates by simple averaging of the upper and lower bounds. This has motivated many different efforts,4,5,12−18 but a real breakthrough did not come about. In a recent letter, I showed how one can in principle systematically improve upon the Temple formula, using information obtained from estimating higher moments of the energy.19 In this present paper, further improvements are presented and shown at least for the quartic oscillator analyzed in detail to lead to tight lower bounds whose accuracy is comparable to the accuracy of the upper bound. The resulting method is based on a number of ingredients. One of them is Temple’s lower bound, which enables us to estimate a lower bound for the overlap squared ⟨φ0|Φ⟩2 of the true ground-state wave function |φ0⟩ with the trial function. The second ingredient is Weinstein’s lower bound, which enables one to obtain a lower bound to the lowest lying excited state energy needed for implementation of Temple’s formula. The third ingredient is a systematic application of the Arnoldi basis set20 which utilizes ever increasing moments of the Hamiltonian operator. The derivation of Temple’s lower bound and an improvement of its estimate to the overlap squared ⟨φ0|Φ⟩2 is presented in section II. The improvement of Temple’s method described in ref 19 is generalized to N dimensional Arnoldi Received: April 8, 2019 Published: June 3, 2019 4079

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087

Article

Journal of Chemical Theory and Computation

Rayleigh−Ritz theorem E0 ≥ ε0, and rearranging give the lower bound estimate for the ground-state

basis sets in section III. To demonstrate the efficacy of the suggested algorithm, the quartic oscillator Hamiltonian is used as a numerical testing ground in section IV. The methodology as presented in sections II and III makes extensive use of the computation of higher moments of the Hamiltonian. It is well-known that the Coulomb potential causes the third and higher moments of the Hamiltonian for atoms and molecules to diverge. Nakatsuji21 has overcome this difficulty by introducing a scaled Hamiltonian in conjunction with the Arnoldi method. An outline of how in principle, the tight lower bound method presented here could be used in combination with the scaled Schrödinger equation and thus applied to atoms and molecules is presented in the Discussion.

ε0 ≥ E0 − H01

(2.11)

2 ⟨Ψ0|Ĥ |Ψ0⟩ ≡ a0ε0 2 + (1 − a0)E2̅

(2.12)

with the notation ∞

E1̅ =



(2.3)

1 [Ĥ − E0I ]̂ |Ψ0⟩ H01

which is orthogonal to the trial state |Ψ0⟩, and one readily finds that

a0 ≥

(2.5)

Multiplying eq 2.4 from the left by the true ground-state of the Hamiltonian and squaring gives the exact result

(2.6)



j=2

(2.7)

so that 2

(1 − a0) [⟨Ψ0|Ĥ |Ψ0⟩ − E1̅ ]2 a0

(E1̅ − E0)2 (E1̅ − E0)2 + H012

(2.15)

H012 E1̅ − E0

(2.16)

There is still one piece of the Temple formula which is missing and that is a lower bound estimate to the energy E̅ 1. As is evident from its definition (eq 2.13) E̅ 1 ≥ ε1. Although the first excited state energy is not known, a lower bound for it may be estimated from Weinstein’s formula. For any function ϕ such that ⟨ϕ|Ĥ |ϕ⟩ ≤ (ε1 + ε2)/2 one obtains the lower bound for the first excited state as

∑ |Ψ⟩⟨Ψ| j j

̂ ≡ |Ψ0⟩⟨Ψ0| + |Ψ⟩⟨Ψ| 1 1 + Q

(2.13)

(1 − a0) [⟨Ψ0|Ĥ |Ψ0⟩ − E1̅ ]2 + (1 − a0)(E2̅ − E1̅ 2) a0

ε0 ≥ ET ≡ E0 −

The identity operator in the Hilbert space of the Hamiltonian may be written as I ̂ = |Ψ0⟩⟨Ψ0| + |Ψ⟩⟨Ψ| 1 1 +

(1 − a0)

Inserting this expression into the lower bound given in eq 2.9 gives the Temple lower bound estimate

2 H012⟨φ0|Ψ⟩ 1

⟨φ0|Ψ0⟩2

, E2̅ =

The inequality in the second line is based on the fact that (E̅ 2 − E̅ 12) is by definition a variance and therefore positive. Since we do not know to estimate this variance reliably, we set it to zero to obtain the inequality in the second line which is equivalent to the lower bound to the square of the overlap matrix element

(2.4)

⟨Ψ|1 Ĥ |Ψ0⟩ = H01



∑ j = 1 ⟨φj|Ψ0⟩2 εj 2

(2.14)

To derive Temple’s lower bound, we introduce the first Arnoldi type20 normalized state

2

(1 − a0)

H012 =

2

[E0 − ε0]2 =

∑ j = 1 ⟨φj|Ψ0⟩2 εj

With some straightforward manipulation of eqs 2.11 and 2.12 one finds that

bounds the true ground-state energy ε0 from above. The same trial function may also be used to compute the standard deviation of the estimate for the ground-state energy

|Ψ⟩ 1 =

⟨Ψ0|Ĥ |Ψ0⟩ = E0 ≡ a0ε0 + (1 − a0)E1̅

(2.2)

⟨Ψ0|Ĥ |Ψ0⟩ − E0 2

(2.10)

For this purpose, we note that the first two moments of the Hamiltonian may be written down as

We also assume that a (normalized) real state |Ψ0⟩ is chosen as a “reasonable” approximation to the ground-state. The energy

H01 =

(2.9)

a0 ≡ ⟨φ0|Ψ0⟩2

(2.1)

E0 = ⟨Ψ0|Ĥ |Ψ0⟩

⟨φ0|Ψ0⟩2

This lower bound expression is not yet practical since one needs to obtain an estimate for the overlap squared

II. IMPROVING TEMPLE’S LOWER BOUND TO THE GROUND-STATE OVERLAP II.a. Temple’s Lower Bound Formula. We assume that the Hamiltonian operator has a complete set of eigenvalues and (normalized) real eigenfunctions Ĥ |φn⟩ = εn|φn⟩

1 − ⟨φ0|Ψ0⟩2

ε1 ≥ ⟨ϕ|Ĥ |ϕ⟩ −

2

|⟨φ0|Ψ⟩| = 1 − ⟨φ0|Ψ0⟩ − ⟨φ0|Q̂ |φ0⟩ ≤ 1 − ⟨φ0|Ψ0⟩ 1

2

⟨ϕ|Ĥ |ϕ⟩ − ⟨ϕ|Ĥ |ϕ⟩2

(2.17)

Temple’s lower bound is then a closed formula based on knowledge of the first two moments of the Hamiltonian using reasonable guesses for the ground and first excited states. However, it is not very accurate due to three reasons. One is ignoring the term ⟨φ0|Q̂ |φ0⟩ in eq 2.8. The second is ignoring the variance as in the first line of eq 2.14. The third is the

(2.8)

The ”heart” of Temple’s lower bound is expressed in the upper bound on the right-hand side, which is based on the fact that the projection operator Q̂ is necessarily non-negative. Inserting this inequality into the identity 2.6, noting that due to the 4080

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087

Article

Journal of Chemical Theory and Computation

to be smaller than the variance H012 associated with the original wave function |Ψ0⟩. In addition, due to the enlarged basis set we may expect a better estimate for the first excited state energy and its variance. Finally, the eigenvalue λ0 ≤ E0 so that we may expect that the lower bound to the overlap squared ⟨φ0|Φ0⟩lb2 obtained from the Temple formula (eq 2.15) will be significantly larger than the lower bound ⟨φ0|Ψ0⟩lb2 obtained for the original wave function using only |Ψ 0 ⟩ and |Ψ 1 ⟩. One will thus usually find that ⟨φ0|Φ0⟩lb2⟨Φ0|Ψ0⟩2 is larger than ⟨φ0|Ψ0⟩lb2. The use of this improved lower bound to the overlap is essential for obtaining a tight lower bound for the ground-state energy.

resulting rather inaccurate lower bound to the overlap matrix element given in eq 2.15. II.b. Improving the Estimate for the Square of the Ground-State Overlap. An improved lower bound to the overlap squared may be derived by considering an orthonormal finite basis set |Ψ0⟩, |Ψ1⟩, ..., |ΨN⟩ which is used to obtain the (N + 1) × (N + 1) matrix representation of the Hamiltonian. The Hamiltonian matrix may be diagonalized to give the N + 1 normalized eigenfunctions |Φ0⟩, |Φ1⟩, ..., |ΦN⟩. Clearly, this basis set will provide better estimates (denoted λ(N+1)j) of the true eigenvalues and eigenfunctions of the Hamiltonian. One may thus expect that the overlap squared ⟨φ0|Φ0⟩2 will be closer to unity than the original overlap squared ⟨φ0|Ψ0⟩2. The identity operator in the N + 1-dimensional space may be written as IN̂ + 1 = |Φ0⟩⟨Φ0| + PN̂

III. TIGHT LOWER BOUNDS As in the previous section, we assume the existence of a trial real normalized state |Ψ0⟩ with an associated mean energy E0 which gives a “reasonable” estimate of the ground-state energy. Following the Arnoldi procedure, we construct a finite set of orthonormal states. For example, the first added state is

(2.18)

so that ⟨φ0|Ψ0⟩2 = ⟨φ0|[|Φ0⟩⟨Φ0| + PN̂ ]|Ψ0⟩2

(2.19)

|Ψ⟩ 1 =

The matrix element ⟨Ψ0|Φ0⟩ is known from the diagonalization of the Hamiltonian matrix. Consider then the form (a + b)2. If both a and b are positive or both negative then (a + b)2 ≥ a2. Projection operators are positive operators. Our original wave function is assumed to be a reasonable approximation to the true ground-state wave function |Ψ0⟩ = |φ0⟩ + |Δφ0⟩

|Ψn⟩ = Π nj =−01

(2.20)

1 det[Hn xn − HÎ n xn]|Ψ0⟩ Hj , j + 1

(3.2)

where we used the notation Hnxn for the matrix representation of the Hamiltonian in the |Ψj⟩ basis set and Inx n is the identity operator in the n dimensional space. For example

⟨φ0|Φ0⟩⟨Φ0|Ψ0⟩ = ⟨φ0|Φ0⟩⟨Φ0|φ0⟩ + ⟨φ0|Φ0⟩⟨Φ0|Δφ0⟩ (2.21)

has two terms. The first one on the right-hand side is positive since it is a diagonal element of the positive projection operator |Φ0⟩⟨Φ0|. The second element is assumed to be small, this is the essence of the statement that the initial wave function is a “good” guess. Assuming that it really is small implies that it cannot change the overall positive sign. The same argument will hold for the matrix element ⟨φ0|P̂ N|Ψ0⟩. In other words when the initial wave function is a “good” guess, it is safe to assume that both terms contributing to the product in eq 2.19 are positive and so use the inequality (a + b)2 ≥ a2. This then implies the important lower bound:

H4x4

ij E0 H01 0 0 yz jj zz jj z jj H01 E1 H12 0 zzz j zz j = jj z jj 0 H12 E2 H23 zzz jj zz jj zz j 0 z 0 H E 23 3 { k

(3.3)

The Hamiltonian matrix in the Arnoldi basis is tridiagonal. The diagonal elements of the matrix are defined as Ei = ⟨Ψ|i Ĥ |Ψ⟩ i , i = 0, 1, ...

(3.4)

with associated variances 2

2 σi 2 = ⟨Ψ|i Ĥ |Ψ⟩ i − Ei , i = 0, 1, ...

(2.22)

(3.5)

The off diagonal coupling matrix elements are readily shown to be

How can we know that the initial wave function is a good guess? One may use the Temple lower bound for the overlap squared (eq 2.15). If it gives a value of 0.9 or thereabouts, then it would be reasonable to assume that indeed it is a good guess. Can one derive strict conditions which would ensure that the correction cannot change the sign? This is much more difficult, and I am not aware of them. Suffice it to say that, typically, initial guesses are very good, with accuracy of at worst a few percent. This is not sufficient for chemical accuracy, but it should be sufficient to ensure that the improved lower bound for the overlap squared is correct. Due to the diagonalization we may expect the variance associated with the improved ground-state wave function

H01 = ⟨Ψ0|Ĥ |Ψ⟩ 1 = σ0 H12 = ⟨Ψ|1 Ĥ |Ψ2⟩ =

(3.6)

σ12 − H012

(3.7)

and more generally Hj , j + 1 = ⟨Ψ|j Ĥ |Ψj + 1⟩ =

σj 2 − Hj − 1, j 2 , j = 0, 1, ...

(3.8)

As an aside, we note that in the previous paper19 there is an error in eqs 16 and 17 where the term σ12 + σ02 should be σ12 − σ02. Consider then the overlap squared of the exact ground-state wave function with |Ψk⟩:

2

σΦ20 = ⟨Φ0|Ĥ |Φ0⟩ − ⟨Φ0|Ĥ |Φ0⟩2

(3.1)

and the standard deviation H01 is as given in eq 2.3. More generally the nth state constructed in this process takes the form

in the sense that the correction term |Δφ0⟩ is small. The product

⟨φ0|Ψ0⟩2 ≥ ⟨φ0|Φ0⟩2 ⟨Φ0|Ψ0⟩2

1 [Ĥ − E0I ]̂ |Ψ0⟩ H01

(2.23) 4081

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087

Article

Journal of Chemical Theory and Computation ⟨φ0|Ψk⟩2 = Πkj =−01

1 Hj , j + 1

2

Equation 3.14 is the central new result of this section. It generalizes and improves upon the Temple lower bound expression in many ways. One notes that the ratio Bn may be

(det[Hk xk − ε0Ik xk ])2 ⟨φ0|Ψ0⟩2 (3.9)

An

considered as an average energy and An as a normalization constant which is larger than unity. Comparing eq 3.14 with the “standard” Temple lower bound as given in eq 2.9, one notes that the first term in the sqrt expression is reduced by the

Unitarity implies that n−1 2 2 ⟨φ0|Ψ⟩ 1 ≤ 1 − ⟨φ0|Ψ0⟩ −

∑ ⟨φ0|Ψ⟩j 2 (3.10)

j=2

factor An. Second, the term

Diagonalizing the Hamiltonian matrix with the orthogonal matrix Onxn, using the notation

OnTxnHn xnOn xn = λn xn

ij λn0 0 jj jj jj 0 λn1 = jjjj jj 0 0 jj jj j 0 0 k

0 yz zz z 0 0 zzzz zz ... 0 zzzz zz 0 λn , n − 1zz {

(3.11)

Using eq 2.6, the inequality of eq 3.10 may be rewritten as H012(1 − ⟨φ0|Ψ0⟩2 )

⟨φ0|Ψ0⟩2 Ä É 2Ñ n−1 Å ÅÅ 2 Ñ ÅÅ k − 1 (λkj − ε0) ÑÑÑ (λk 0 − ε0) 2 Ñ − H01 ∑ ÅÅÅΠ j = 1 Ñ Å Hj , j + 12 ÑÑÑÑ H012 k=2 Å ÅÇ Ö H012(1 − ⟨φ0|Ψ0⟩2 )

⟨φ0|Ψ0⟩2 ÄÅ É 2Ñ ÅÅ Ñ ÅÅ k − 1 (λkj − λn0) ÑÑÑ ÑÑ(λk 0 − ε0)2 − ∑ ÅÅÅΠ j = 1 2 ÑÑ Å H ÑÑÖ j , j + 1 k=2 Å ÅÇ

1 d2 1 Ĥ = − + x4 2 2 dx 2

(3.13)

where in the second line we used the fact that the ground-state eigenvalue λn0 ≥ ε0 for any dimension n. This gives us a quadratic inequality for the ground-state energy ε0 which is readily solved: H012(1 − ⟨φ0|Ψ0⟩2 ) A n⟨φ0|Ψ0⟩2



Table 1. First Few Symmetric Eigenvalues for the Quartic Oscillator, Adapted from Reference 22a

A n2 (3.14)

(3.15)

∑ αj

iΓy Ψ0(x) = jjj zzz kπ {

∑ λj0αj

1/4

(3.17)

j=2 n−1

Cn =

+

2

∑ λj0 αj j=2

0.53018105 3.72784897 8.13091301 13.26423554 18.96150051

As a trial wave function, we use the Gaussian

n−1

E02

εm

0 2 4 6 8

(3.16)

j=2

Bn = E0 +

m

a m denotes the number of nodes in the wavefunction; εm is the energy of the mth state.

n−1

An = 1 +

(4.1)

The lowest symmetric eigenvalues have been reported in ref 22 and are given in Table 1.

A nCn − Bn 2

where we used the notation ÄÅ É 2Ñ ÅÅ Ñ ÅÅ k − 1 (λkj − λn0) ÑÑÑ ÑÑ, k = 2, ..., n − 1 αk = ÅÅÅΠ j = 1 2 ÑÑ ÅÅ H ÑÑÖ j , j + 1 ÅÇ

⟨φ0|Ψ0⟩2

IV. TIGHT BOUNDS FOR THE GROUND-STATE OF A QUARTIC OSCILLATOR IV.a. Preliminaries. The quartic Hamiltonian we use is

n−1

Bn − An

(1 − ⟨φ0|Ψ0⟩2 )

may be estimated as described in the previous section. To exemplify the improved quality of the present lower bound expression, we consider in detail in the next section the specific case of a quartic oscillator. It should be stressed that in principle the lower bound estimate described in this paper will always converge to the correct answer. This is in a sense a trivial observation. As the size of the Hamiltonian matrix increases, the estimate for the ground-state wave function converges to the true ground-state wave function, the variance vanishes, and already from Temple’s formula, one obtains the exact ground-state energy. When using the original Temple formula, this convergence is typically very slow, and the size of the basis set needed to obtain useful lower bounds becomes too large. The methodology presented in this paper is useful since the convergence is much faster than that obtained using only the Temple estimate and as we shall see in the next section for the quartic oscillator can become comparable to the convergence of the upper bound estimate.

(3.12)

ε0 ≥

is a variance and therefore

thus increases the lower bound. Finally, the ratio

0

det[Hn xn − ε0In xn] = det[λn xn − ε0In xn] = Π nj =−01(λnj − ε0)



A n2

positive, and this further reduces the term under the sqrt and

we have that

(E0 − ε0)2 ≤

A nCn − Bn 2

exp( −Γx 2/2)

(4.2)

where the width parameter Γ will be determined further below. The first excited trial state has the form

(3.18) 4082

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087

Journal of Chemical Theory and Computation Ä ÉÑ ÑÑ 1 ÅÅÅÅ 1 2 2 4 Ψ1(x) = ÅÅ (Γ − Γ x + x ) − E0(Γ)ÑÑÑΨ0(x) ÑÖ H01(Γ) ÅÇ 2

Article

The lowest excited state energy is found to be (4.3)

10Γ 9 + 3Γ 6 − 330Γ 3 + 1224 E1(Γ) = ⟨Ψ|1 Ĥ |Ψ⟩ = 1 8Γ 2(Γ 6 − 6Γ 3 + 12)

with 2Γ 3 + 3 2 E0(Γ) = ⟨Ψ0|Ĥ |Ψ0⟩ = , H012(Γ) = ⟨Ψ0|Ĥ |Ψ⟩ 1 2 8Γ Γ 6 − 6Γ 3 + 12 = (4.4) 8Γ 4

(4.5)

The associated variance is

2

2 σ12(Γ) = ⟨Ψ|1 Ĥ |Ψ⟩ 1 − E1(Γ)

=

7Γ18 − 150Γ15 + 1806Γ12 − 13824Γ 9 + 66132Γ 6 − 168912Γ 3 + 181278 8Γ 4(Γ 6 − 6Γ 3 + 12)2

(4.6)

Φ21(x). Now, the largest lower bound is found at Γ ≃ 2.10 for which E̅ 1,W ≃ 2.758. This then will be the value of the estimate of the first excited state energy which will be used when obtaining the lower bound from the two-dimensional basis set computation (n = 2). Using E̅ 1,W = 2.758, we may obtain the ”standard” Temple lower bound (eq 2.16) plotted in Figure 2 (blue dotted line) as

At this point, we may obtain Weinstein’s estimate for Temple’s excited state energy as in eq 2.17. In Figure 1, we

Figure 1. Weinstein lower bound estimate for the lowest excited state energy. The dotted (red) and dashed (blue) lines respectively show the energies E1(γ) and Weinstein’s lower bound E̅ 1,W(Γ) as functions of the width parameter Γ. Weinstein’s lower bound is valid only in the region where the energy E1(Γ) is lower than the average of the first and second excited state eigenvalues, shown as the black solid line. The largest valid lower bound in this initial computation is thus found at Γ = 1.983 for which E̅ 1,W = 2.614977. For further explanation see the text.

Figure 2. Upper and lower bound estimates of the ground-state energy of the quartic oscillator. The (red) dashed-dotted and (orange) long dashed lines are the upper bounds E0(Γ) and λ20(Γ) respectively. The dotted (blue) line is the Temple lower bound obtained from eq 2.16, the dashed (green) line is the Temple lower bound obtained after diagonalization of the two-dimensional representation of the Hamiltonian (eq 4.7). The solid (black) line is the exact ground-state energy. Note how the upper bound improves significantly after diagonalization of the two-dimensional Hamiltonian matrix while the Temple lower bound much less so.

plot E1(Γ) and the resulting Weinstein lower bound E̅ 1,W = E1(Γ) − σ1(Γ) as functions of the width parameter Γ. Since Weinstein’s formula provides a lower bound to the first excited state energy provided that the estimate E1(Γ) ≤ (ε1 + ε2)/2, one notes from the figure that only values of Γ ≥ ∼ 1.7 may be considered. In this region, the largest lower bound is found at Γ ≃ 1.983 for which E̅ 1,W ≃ 2.614977. We note though that to obtain this lower bound it was necessary to consider the first four moments of the Hamiltonian. Putting it differently, it was necessary to create the first Arnoldi excited state wave function Ψ1(x) as in eq 4.3. One may then obtain a better lower bound on the lowest excited state energy by constructing the two-dimensional matrix of the Hamiltonian using the functions Ψ0(x) and Ψ1(x) and diagonalizing it. This will given the improved eigenvalues λ20 and λ21 and their associated eigenfunctions Φ20(x) and Φ21(x) (the first index denotes the number of basis functions used to obtain the diagonalized function.). One may then obtain an improved estimate of the Weinstein lower bound which now takes the form ε1 ≥ λ21 − σλ21 where σλ212 is the variance obtained with the excited state wave function

a function of the width parameter Γ. The best lower bound to the ground-state energy based on Temple’s expression (eq 2.16) is found in the vicinity of Γ = 1.616 and its value is E0,LBT = 0.5078. As may be discerned from eq 4.4, the best upper bound for E0(Γ) (red dashed-dotted line in Figure 2 is found at Γ = 31/3 and its value is ∼0.5408. The distance of this upper bound from the exact result is ∼0.0107 which is a factor of ∼2 smaller than the distance (∼0.0224) of the exact eigenvalue from the best Temple lower bound. After diagonalizing the two-dimensional representation of the Hamiltonian one obtains improved upper (λ20) and lower bounds (σλ202 is the variance of the Hamiltonian using the diagonalized ground-state wave function Φ20(x)): 4083

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087

Article

Journal of Chemical Theory and Computation λ 20,LBT = λ 20 −

σλ20 2 λ1̅ − λ 20

including the third moment of the Hamiltonian which is used to obtain E1(Γ). 3. Diagonalize the two-dimensional Hamiltonian and use Weinstein’s lower bound formula (eq 2.17) to estimate the best lower bound to the lowest excited state energy. 4. Use the ground-state eigenvalue and eq 3.14 to obtain improved upper and lower bounds. This calls for up to four moments of the Hamiltonian since the variance of the upper state Φ21 depends on the fourth moment. The next round of operations is to create the threedimensional representation of the Hamiltonian. The steps are 1. Create the second Arnoldi wave function Ψ2(x). 2. Construct the 3-dimensional matrix representation of the Hamiltonian using the Arnoldi wave functions Ψ0(x), Ψ1(x), and Ψ2(x). 3. Diagonalize the Hamiltonian and use Weinstein’s lower bound expression to improve the lower bound for the first excited state. 4. Use the ground-state eigenvalue and eq 3.14 to obtain improved upper and lower bounds. This calls for six moments since the variance of the second excited state depends on the sixth moment. The results of this second round are denoted as n = 3, since this involves the three-dimensional Arnoldi representation of the Hamiltonian. They are shown in Figure 4. The dotted lines

(4.7)

where λ̅1 is the analog of E̅ 1 and is also bounded from below by the first excited state energy ε1. These upper and lower bounds are also plotted in Figure 2. The upper (λ20, long dashed orange line) and lower Temple bound (λ20,LBT, dashed green line) are shown as a function of the width parameter Γ. One notes the improvement obtained from diagonalization, but at the same time, one also sees that the quality of the upper bound is superior to that of the lower bound. The best upper bound found is λmin 20 = 0.53213854 at Γ = 1.40647, while the best lower bound λmax 20,LBT = 0.51947684 is found at Γ = 1.51302. Both upper and lower bounds are improved, but the upper bound much more than the lower bound. The distance of the best upper bound from the exact value (0.00196) is now a factor of ∼5 smaller than the distance of the exact eigenvalue from the optimal Temple lower bound (0.0107). It is this poor convergence of the lower bound which needs improvement. To demonstrate the efficacy of the method developed in this paper, we plot in Figure 3 the dependence of the lower bound

Figure 3. Two-dimensional representation based upper and lower bounds for the ground-state energy as a function of the Gaussian width parameter Γ. The (orange) long dashed line is the upper bound λ20(Γ) (also shown in Figure 2). The dashed (green) line is the Temple lower bound obtained after diagonalization of the twodimensional representation of the Hamiltonian (eq 4.7, shown also in Figure 2). The (blue) dotted line is the lower bound derived in this paper (eq 3.14) as applied to the two-dimensional representation of the Hamiltonian. Note that this last result is competitive in accuracy with the upper bound.

Figure 4. Comparison of two and three-dimensional based upper and lower bounds for the ground-state energy as a function of the Gaussian width parameter Γ. The dotted (blue) lines show the upper and lower bounds based on the two-dimensional representation of the Hamiltonian, the dashed dotted (red) lines show the same for the three-dimensional representation. The lower bounds are obtained from eq 3.14.

obtained from eq 3.14 for the two-dimensional case and compare it to Temple’s lower bound obtained for the same (eq 4.7) and plotted also as the green dashed line in Figure 2. The best lower bound (0.52623479) is found at Γ = 1.9264. Its distance from the exact eigenvalue (0.00395) is now only a factor of 2 larger than the distance of the upper bound from the exact eigenvalue. The present method has improved the lower bound estimate by a factor of ∼2.5. IV.b. Tight Bounds. The process of obtaining tight bounds whether upper or lower is to systematically increase the dimensionality of the basis set to be used. In our case, this means starting with the Gaussian wave function (eq 4.2). The steps for the first round as described above were as follows: 1. Use the Gaussian to bound the ground-state energy from above and find the lowest upper bound as a function of the width parameter Γ. 2. Obtain the two-dimensional matrix elements as functions of Γ. For this purpose one needs up to and

depict the upper and lower bounds for the two-dimensional case, the dashed dotted lines for the three-dimensional case. The upper bounds show the eigenvalues λ20 and λ30. The lower bounds are the results obtained from eq 3.14 using the values of E̅ 1 = 2.758 and 2.973 for the cases n = 2 and 3, respectively. This figure clearly shows how the improved lower bounds converge much more rapidly to the eigenvalue and that this convergence is comparable to the same for the upper bounds. This process is then systematically repeated for increasing dimensionality of the basis set. The best Weinstein lower bounds to the first excited state energy are given in Table 2 for the dimensions n = 2−5. The optimal lower and upper bounds to the ground-state for a dimensionality of up to n = 5 are given in Table 3. From these results, it is evident that the improved lower bound, as expressed in eq 3.14 becomes very accurate as the dimensionality increases. The distance of the lower bounds 4084

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087

Article

Journal of Chemical Theory and Computation Table 2. Weinstein Lower Bounds to the First Excited State Energy as a Function of the Dimensionality of the State Space Useda n

E1,W(n)

2 3 4 5

2.758 2.973 3.388 3.576

a Note that by n = 5 one is reasonably close to the exact first excited state energy (∼3.728).

Figure 5. Average of the lower and upper bounds is plotted as a function of the width parameter Γ. The solid blue line is the average for the n = 3-dimensional case, the upper and lower dashed dotted lines are the upper and lower bounds, and the solid black line is the exact ground-state energy.

from the exact energy is very similar to that of the upper bounds. The tightness of the lower bound implies that by averaging the upper and lower bounds one obtains a much better estimate of the ground-state energy. In Figure 5 we plot the average (solid blue line) of the lower and upper bounds as a function of Γ and compare the mean with the upper and lower bounds for the n = 3-dimensional case. The average is much closer to the exact result (solid black line) for a large range of the parameter space, reflecting the fact that the upper and lower bounds are comparable in quality. The mean values of the best upper and lower bounds as a function of the dimensionality n are tabulated in Table 3 together with their distance from the exact answer. As expected, due to the tightness of both upper and lower bounds, the gap of the average value is reduced by a factor of ∼5−10. It took only five iterations to reduce the error in the estimate of the eigenvalue by 6 orders of magnitude.

Figure 6. Temple and improved lower bounds to the overlap squared of the approximate ground-state wave function with the exact groundstate wave function. The two lower lines show the Temple lower bound to the overlap squared (eq 2.15) for n = 2 (blue dotted line) and n = 3 (blue long dashed line) dimensions. The two upper lines show the lower bound to the overlap squared as obtained from eq 2.22. The lower (green) long dashed line is for n = 2, the upper (green) solid line is for n = 3.

V. DISCUSSION A general theory has been presented for obtaining systematic improvement of a lower bound to the ground-state energy, based on the use of the Arnoldi method of creating a basis set. Two essential ingredients are needed. One is the inclusion of higher-order overlaps, as discussed in section III. The other is obtaining an improved estimate for the overlap squared of the exact ground-state wave function with the zeroth order trial function as discussed in section 2. The Temple lower bound for the overlap squared can only be improved by obtaining a higher lower bound for the lowest excited state energy. In contrast the improved overlap method (eq 2.22) gives a lower bound to the overlap squared (a0,P) which increases systematically with increasing dimensionality and thus contributes significantly to the improved lower bound estimate of the ground-state energy. This may be seen more clearly by inspection of Figure 6 where the Temple lower bounds to the overlap squared (eq 2.15) are plotted as a function of Γ for n = 2 (blue dotted line) and n = 3 (blue long dashed line)

dimensions. The improvement of Temple’s lower bound to the overlap squared for n = 3 as compared to n = 2 is a reflection of the increased estimate of the lower bound to the excited state energy. The same figure also shows the lower bounds to the overlaps squared, based on eq 2.22 for n = 2 (green long dashed upper line) and n = 3 (green solid line). The improvement as compared with Temple’s lower bound is evident and does not only reflect the increase in the estimate for the excited state energy. Although the lower bound presented in this paper is competitive with the upper bound, one should note that computation of the upper bound always demands one order

Table 3. Tight Lower Bounds for the Quartic Oscillatora n

Γ*UB

λn0(Γ*UB)

Γ*LB

ELB(Γ*LB)

⟨E0⟩

ΔUB

ΔLB

Δ⟨E0⟩

2 3 4 5

1.40647 2.13384 2.33434 2.53905

0.53213854 0.53068300 0.53023261 0.53018746

1.9264 2.2267 2.4144 2.6005

0.52634791 0.52944139 0.53012022 0.53017226

0.52924322 0.53006220 0.53017642 0.53017986

0.00195749 0.00050195 0.00005156 0.00000641

0.00383314 0.00073966 0.00006083 0.00000874

0.00093782 0.00011886 0.00000464 0.00000120

a n denotes the dimensionality of the Hamiltonian matrix, ΓUB * and ΓLB * are the width parameter values at which the best upper and lower bounds are found, respectively. λn0(Γ*UB) and ELB(Γ*LB) denote the values of the best upper and lower bounds to the ground-state energy. ⟨E0⟩ is the average of the lower and upper bound energies. ΔUB and ΔLB are the distance of the upper and lower bounds respectively from the exact ground-state eigenvalue, Δ⟨E0⟩ is the distance of the mean energy from the exact ground-state eigenvalue.

4085

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087

Article

Journal of Chemical Theory and Computation

obtain an improved estimate for ⟨φ0|Ψ0⟩2 using eq 2.22. This would then provide through eq 2.9 an improved lower bound estimate for the ground-state energy. But this is only half of the new methodology; one would still want to obtain an expression for the overlaps ⟨φ0|Ψk⟩2 so as to include also these higher-order contributions. Here we note that when using the scaled Hamiltonian procedure, all the Arnoldi like states created have finite first and second moments of the Hamiltonian. Inspection of the final lower bound formulas (eqs 3.14−3.18) shows that all that is needed to implement the lower bound are first and second moments of the Hamiltonian and eigenvalues of the Hamiltonian matrix. In other words, although formally the Arnoldi basis set without scaling would lead to divergences, in practice these are not necessarily a problem. The nth moment of the Hamiltonian taken with an exact eigenstate does not diverge. The lower bound is based on overlaps of approximate wave functions with the exact ground-state wave function so that here too the higher moments appear as higher powers of the ground-state energy and do not diverge. Implementation of this procedure to an atom such as He will be the topic of future work. In summary, a procedure which can lead to tight lower bounds to the ground-state energy has been presented and tested. The results are promising and suggest that, in the future, one will be able to obtain both lower and upper tight bounds to atomic and molecular energies.

less in the computation. That is, when the lower bound calls for 2n moments of the Hamiltonian, the upper bound needs only 2n − 1 moments. However, especially as the basis set is increased, this difference is not essential. A by product of these results is that a simple average of the upper and lower bound estimate of the eigenvalue will significantly improve the accuracy of the estimate. The critical reader might question how general the results presented for the quartic oscillator are. In this context, one should note that eq 2.6 is exact. To derive Temple’s lower bound one needs three ingredients. One is the overlap of the first Arnoldi state with the exact ground-state. To obtain Temple’s result, one uses the simple estimate as given in eq 2.8. The results of Section III show that the higher-order contributions need not be ignored. Using the Arnoldi basis, they may all be expressed in terms of the overlap of the approximate ground-state wave function with the exact ground-state function. In other words, the first source of lack of tightness in Temple’s expression has been removed. The second source of inaccuracy is the quality of the lower bound to the overlap matrix element squared. The improvement given in eq 2.22 can only increase the accuracy of the lower bound estimate. It has very little to do with the specifics of the system studied but rather with the quality of the functions used. If the initial guess is of high quality, then eq 2.22 should significantly improve the lower bound to the overlap squared. The only drawback of Temple’s lower bound which has not been improved in this paper is the estimate for the lowest excited state energy. This is not much of an impedance since the lower bound is not very sensitive to the precise magnitude of the estimate. Suppose that E̅ 1 = ε1 − δε1. Then the effect on the lower bound relative to the exact excited state energy is ∼|δε1|/(ε1 − ϵ0). This is typically a small number unless the first excited state is almost degenerate with the ground state and so does not affect the lower bound too much. The reasons underlying the tightness of the lower bound for the quartic oscillator are not specific to the potential but reflect the improved theory. In addition, one may further use the present theory to improve the lower bound estimate to the first excited state by using a rough estimate to the energy of the second excited state. This should be useful especially when the first two states are very close to each other in energy. The tight lower bound as presented especially in section IV calls for the computation of moments of the Hamiltonian which are three or higher. For atoms with Coulomb potentials, these moments diverge. This problem has been neatly solved by Nakatsuji and co-workers by introducing their scaled Hamiltonian.21 However, due to the scaling, eq 3.9 relating the overlap of the ground-state with the higher Arnoldi states ⟨φ0|Ψk⟩2 to the overlap with the initial state ⟨φ0|Ψ0⟩2 is no longer straightforward. It is though possible to implement the tight lower bound method considered in this paper even for Coulombic systems. First we note that the lower bound for the overlap squared as given in eq 2.22 is readily implemented. Specifically, let us suppose that one carries out a scaled Hamiltonian computation which provides good wave functions for the first N states of the Hamiltonian, denote them as |Ψj⟩. Since they are not exact, one may use them to construct an NxN matrix representation of the Hamiltonian. This may then be diagonalized to obtain an improved basis set |Φj⟩ such that ⟨Φ0|Ψ0⟩ is known. Furthermore, one may use the Temple relation (eq 2.15) to obtain a lower bound to the overlap of the exact ground-state wave function from the state |Φ0⟩ and so



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Eli Pollak: 0000-0002-5947-4935 Notes

The author declares no competing financial interest.



ACKNOWLEDGMENTS This work was generously supported by a grant of the YedaSela Foundation of the Weizmann Institute of Science.



REFERENCES

(1) Nakashima, H.; Nakatsuji, T. How Accurately Does the Free Complement Wave Function of a Helium Atom Satisfy the Schrödinger Equation? Phys. Rev. Lett. 2008, 101, 240406. (2) Temple, G. The Theory of Rayleigh’s Principle as Appliedto Continuous Systems. Proc. R. Soc. London, Ser. A 1928, A119, 276. (3) Delves, L. On the Temple Lower Bound for Eigenvalues. J. Phys. A: Gen. Phys. 1972, 5, 1123−1130. (4) Hill, R. N. Tight Lower Bounds to Eigenvalues of the Schrödinger Equation. J. Math. Phys. 1980, 21, 2182−2192. (5) Marmorino, M. G.; Gupta, P. Surpassing the Temple Lower Bound. J. Math. Chem. 2004, 35, 189−197. (6) Donchev, A. G.; Kalachev, S. A.; Kolesnikov, N. N.; Tarasov, V. I. The Upper and Lower Bounds of Energy for Nuclear and Coulomb Few-Body Systems. Phys. Part. Nucl. Lett. 2007, 4, 39−45. (7) Weinstein, D. H. Modified Ritz Method. Proc. Natl. Acad. Sci. U. S. A. 1934, 20, 529−532. (8) Stevenson, A. F. On the Lower Bounds of Weinstein and Romberg in Quantum Mechanics. Phys. Rev. 1938, 53, 199. (9) Schollwöck, U. The Density-matrix Renormalization Group in the Age of Matrix Product States. Ann. Phys. 2011, 326, 96−192. (10) Baumgratz, T.; Plenio, M. B. Lower Bounds for ground-states of Condensed Matter Systems. New J. Phys. 2012, 14, 023027. 4086

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087

Article

Journal of Chemical Theory and Computation (11) Ren, J.; Yi, Y.; Shuai, Z. Inner Space Perturbation Theory in Matrix Product States: Replacing Expensive Iterative Diagonalization. J. Chem. Theory Comput. 2016, 12, 4871−4878. (12) Cohen, M.; Feldmann, T. Lower Bounds to Eigenvalues. Can. J. Phys. 1969, 47, 1877−1879. (13) Scrinzi, A. Lower Bounds to the Binding Energies of tdμ. Phys. Rev. A: At., Mol., Opt. Phys. 1992, 45, 7787−7791. (14) Lüchow, A.; Kleindienst, H. Accurate Upper and Lower Bounds to the 2S States of the Lithium Atom. Int. J. Quantum Chem. 1994, 51, 211−224. (15) Marmorino, M. G. Equivalence of Two Lower Bound Methods. J. Math. Chem. 2002, 31, 197−203. (16) Marmorino, M. G.; Almayouf, A.; Krause, T.; Le, D. Optimization of the Temple Lower Bound. J. Math. Chem. 2012, 50, 833−842. (17) Toth, S.; Szabados, A. Energy Error Bars in Direct Configuration Interaction Iteration Sequence. J. Chem. Phys. 2015, 143, 084112. (18) Marmorino, M. G.; Black, V. Lower Bounds to the GroundState Expectation Value of Non-negative Operators. J. Math. Chem. 2016, 54, 1973−1985. (19) Pollak, E. An Improved Lower Bound to the ground-state Energy. J. Chem. Theory Comput. 2019, 15, 1498−1502. (20) Arnoldi, W. E. The Principle of Minimized Iterations in the Solution of the Matrix Eigenvalue Problem. Q. Appl. Math. 1951, 9, 17−29. (21) Nakatsuji, H. Scaled Schrödinger Equation and the Exact Wave Function. Phys. Rev. Lett. 2004, 93, 030403. (22) Reid, C. E. Energy Eigenvalues and Elements for the Quartic Oscillator. J. Mol. Spectrosc. 1970, 36, 183−191.

4087

DOI: 10.1021/acs.jctc.9b00344 J. Chem. Theory Comput. 2019, 15, 4079−4087