Efficient description of strongly correlated electrons with mean-field cost

0.54 MB·

by  Katharina Boguslawski

Efficient description of strongly correlated electrons with mean-field cost
Katharina Boguslawski, Pawel(cid:32) Tecmer, and Paul W. Ayers∗
Department of Chemistry and Chemical Biology,
McMaster University, Hamilton,
1280 Main Street West, L8S 4M1, Canada
Patrick Bultinck
Department of Inorganic and Physical Chemistry,
Ghent University, Krijgslaan 281 (S3), 9000 Gent, Belgium
4
Stijn De Baerdemacker and Dimitri Van Neck†
1
Center for Molecular Modelling, Ghent University,
0
Technologiepark 903, 9052 Gent, Belgium
2

g (Dated: August 5, 2014)
u
A Wepresentanefficientapproachtotheelectroncorrelationproblemthatiswell-suitedforstrongly
interacting many-body systems, but requires only mean-field-like computational cost. The perfor-

2 mance of our approach is illustrated for the one-dimensional Hubbard model for different ring
lengths, and for the non-relativistic quantum chemical Hamiltonian exploring the symmetric disso-

] ciation of the H hydrogen chain.
l 50
e
-
r
t The accurate description of the electron–electron in- tion [34], i.e.
s
teraction at the quantum-mechanical level is a key prob-
at. lemincondensedmatterphysicsandquantumchemistry. |Ψ (cid:105)=exp(cid:32)(cid:88)P (cid:88)K caa† a† a a (cid:33)|Φ (cid:105), (1)
m Since most of the quantum many-body problems are ex- AP1roG i a↑ a↓ i↓ i↑ 0
traordinarily difficult to solve exactly, different approxi- i=1a=P+1
-
d mation schemes emerged [5–9], among which the density
where a† and a (σ =↓,↑) are the fermionic cre-
n matrix renormalization group (DMRG) algorithm [10– pσ pσ
o 12] gained a lot of popularity in both condensed mat- ation and annihilation operators, and |Φ0(cid:105) is some
c independent-particle wavefunction (usually the Hartree–
ter physics [11] and quantum chemistry [13–20] over the
[ Fockdeterminant). Indicesiandacorrespondtovirtual
lastdecade. SincetheDMRGalgorithmoptimizesama-
and occupied sites (orbitals) with respect to |Φ (cid:105), P and
3 trix product state wavefunction, it is optimally suited 0
v for one-dimensional systems; though DMRG studies on K denote the number of electron pairs (P = N/2 with
9 N being the total number of electrons) and orbitals, re-
higher-dimensional and compact systems have been re-
1 spectively, and {ca} are the geminal coefficients. This
ported [14, 17, 19, 21]. Yet, novel theoretical approaches i
0 wavefunction ansatz is size-extensive and has mean-field
8 aredesirablethatcanaccuratelydescribestrongcorrela- scaling, O(cid:0)P2(K −P)2(cid:1) for the projected Schr¨odinger
. tioneffectsbetweenelectronswherethedimensionofthe
1 equation approach [31].
Hilbert space exceeds the present-day limit of DMRG or
0 To ensure size-consistency, however, it is necessary
general tensor-network approaches [22] allowing approx-
4 to optimize the one-electron basis functions [31], where
1 imately 100 sites or 60 (spatial) orbitals, respectively.
all non-redundant orbital rotations span the occupied–
:
v Another promising approach, suitable for larger occupied, occupied–virtual, and virtual–virtual blocks
i strongly-correlated electronic systems, uses geminals with respect to the reference Slater determinant |Φ (cid:105).
X 0
(two-electron basis functions) as building blocks for the We have implemented a quadratically convergent algo-
r
a wavefunction [23–30]. One of the simplest practical rithm: weminimizetheenergywithrespecttothechoice
geminal approaches is the antisymmetric product of 1- of the one-particle basis functions, subject to the con-
reference-orbital geminals (AP1roG) [31–33]. Unique straint that the projected Schr¨odinger equations for the
among geminal methods, AP1roG can be rewritten as geminal coefficients hold. Specifically, we use a Newton–
a fully general pair-coupled-cluster doubles wavefunc- Raphson optimizer and a diagonal approximation of the
orbital Hessian to obtain the rotated set of orbital ex-
pansion coefficients. Our algorithm is analogous to the
orbital-optimized coupled cluster approach [35–37]. Due
to the four-index transformation of the electron repul-
sion integrals, the computational scaling deteriorates to
∗ [email protected] O(cid:0)K5(cid:1). The orbital-optimized AP1roG (OO-AP1roG)
† [email protected] approach was implemented in a developer version of the
2
00..0011
100.0
e [t]e [t] 00..0000
sitsit −−00..0011 95.0
r r
pepe −−00..0022 90.0
) )
GG
oo −−00..0033 κ 85.0
1r1r %
EEAPAP −−00..0044 80.0
−− −−00..0055
actact 75.0
EEExEx −−00..0066 NNssiitteess== 160 NNssiitteess== 160
( (∆∆ −−00..0077 NNssiitteess== 5104 70.0 NNssiitteess== 5104
−−00..0088 Nsites= 122 65.0 Nsites= 122
00..0011 00..11 11 22 44 88 1100 2200 3300 5500 110000 220000 0.001 0.01 0.1 1 2 4 8 10 20 30 50 100 200
UU [[tt]] U [t]
FIG. 1. Deviation of the OO-AP1roG total energies from FIG. 2. Percentage of the correlation energy %κ for different
exact values (blue dashed line) for different strengths of the strengths of the repulsive on-site interaction in the half-filled
repulsiveon-siteinteractionforthe1-DHubbardmodel(with 1-D Hubbard model (with periodic boundary conditions) for
periodic boundary conditions) for N = 6,10,14,50,122. N =6,10,14,50,122capturedbyOO-AP1roG.Theexact
sites sites
TheexactvaluesforsmallU (U <0.001t)forN =50,122 valuesforsmallU (U <0.001t)forN =50,122couldnot
sites sites
could not be converged. be converged.
Horton program package [38]. (see Figure I in the Supplementary Information). In
(a) The half-filled one-dimensional Hubbard the limit U → 0, the geminal coefficient matrix cor-
Hamiltonian. First, we consider the 1-D Hubbard rectlyapproachesthezeromatrixindicatingthatasingle
model Hamiltonian with periodic boundary conditions, Slater determinant is sufficient to describe the quantum
state exactly. For increasing U, {ca} becomes diagonal-
HˆHub =−t(cid:88)(cid:16)a†(j+1)σajσ+a†jσa(j+1)σ(cid:17)+U(cid:88)nj↑nj↓, dominant and adopts a diagonal striucture in the limit of
j,σ j U →∞. Thus,inthelimitofinfinite(repulsive)interac-
(2) tion, OO-AP1roG optimizes a perfect-pairing (seniority-
wherethefirsttermrepresentsnearest-neighborhopping zero) wavefunction [40, 41],
and the second term is the repulsive on-site interaction.
The operators a†jσ and ajσ are again the fermionic cre- (cid:89) [(a†i,↑+a†i+1,↑)(a†i,↓+a†i+1,↓)−(a†i,↑−a†i+1,↑)(a†i,↓−a†i+1,↓)]|0(cid:105)
ation and annihilation operators on a lattice with sites
i=1,3,..
j = 1,...,N , and n = a† a is the local number (3)
sites jσ jσ jσ
operator. To conclude, OO-AP1roG has mean-field-like scaling,
Figure 1 shows the differences in total energies ob- but can recover about 71% of the correlation energy in
tained for OO-AP1roG with respect to reference data the weak interaction regime, about 80% for intermediate
obtainedfromthesolutionoftheLieb-Wuequations[39] interactionstrengths,andapproximately93%inthecase
(Nsites = 6,10,14,50,122). OO-AP1roG can reproduce of strong on-site interaction for all chain lengths studied
the exact total energies in the limit of zero and infinite (a numerical comparison is presented in Table I of the
(repulsive) on-site interaction. The largest deviations Supplementary Information).
from the exact solution (up to 0.075t per site) are found Figure 3 shows the single-orbital entropy for different
fortheintermediateregionoftheon-siteinteraction,that lengths of the 1-D lattice as a function of the repulsive
is, for 2t < U < 50t. Figure 2 shows the percentage of on-site interaction U. The single-orbital entropy is the
the correlation energy captured by OO-AP1roG calcu- analogue of the one-site entropy, but determined in the
lated as %κ= EOO-AP1roG−EHF ·100. In the limit of zero naturalorbitalbasis: itiscalculatedasthevonNeumann
Eexact−EHF
andinfiniteU,theOO-AP1roGmodelbecomesexact;for entropy from a single-orbital density matrix (a many-
U = 0 the wavefunction can be exactly described by a particle reduced density matrix of one orbital). It mea-
singleSlaterdeterminantandthusthecorrelationenergy surestheentanglementofoneorbitalwiththeremaining
approaches zero, while for U → ∞, the quantum state (N −1) ones [18]. In particular, since the optimized
sites
can be represented by the perfect pairing wavefunction. orbitals are localized on two neighboring sites, the von
For growing (repulsive) U, the percentage of the correla- Neumann entropy describes the correlation of a pair of
tion energy covered by OO-AP1roG increases gradually. sites and the other part of the system. In the following,
For small values of U, the geminal coefficient matrix we will refer to the single-orbital entropy as the pair en-
{ca} is sparse and thus far from perfect pairing, which tanglementE inaccordancewiththelocalentanglement
i p
is represented by a diagonal geminal coefficient matrix determined for the on-site basis [44].
3
RHF
0.70 −18.00 AP1roG
OO−AP1roG
DMRG
0.60 CCSD(T)
0.30 e] −20.00 MP2
e PBE
0.50 0.20 Hartr −22.00
Ep 0.40 0.10 gy [
0.30 0.00 ner −24.00
0 1 2 3 4 5 E
0.20 Nsites= 6 −26.00
N = 10
sites
N = 14
sites
0.10 NNsites== 51022 −28.00
lns(2ite)s 1 1.5 2 2.5 3 3.5 4
0.00 H−H distance [Å]
0 2 4 6 8 10 20 50 100 200
U [t]
FIG. 4. Symmetric dissociation of H chain using the STO-
50
6G basis set [49] obtained from different methods. The
FIG. 3. Pair entanglement E (single-orbital entropy) for
p DMRG reference data are taken from Ref. 45.
different strengths of the repulsive on-site interaction in the
half-filled1-DHubbardmodel(withperiodicboundarycondi-
tions)forN =6,10,14,50,122calculatedbyOO-AP1roG.
sites
TABLEI.Spectroscopicconstants: equilibriumbonddistance
(R ), potential energy depth (D ) and harmonic vibrational
e e
frequency(ω )forthegroundstateofH (STO-6G).Differ-
The pair entanglement takes its minimum value at e 50
ences with respect to the DMRG reference data are listed in
U = 0t where the wavefunction can be exactly repre-
parenthesis.
sented by a single Slater determinant. It is easy to
verify that all orbital pairs are uncorrelated in a one- Method Re [˚A] De [eV] ωe [cm−1]
RHF 0.940(−0.030) 199.0(+109.3) 25089(+2268)
determinant wavefunction and thus Ep =0. For increas- AP1roG 0.941(−0.029) 198.2(+108.5) 23013(+2252)
ing on-site interaction, the pair entanglement smoothly MP2 0.955(−0.015) 144.1(+54.4) 24568(+1747)
PBE 0.971(+0.001) 146.6(+56.9) 23662(+841)
accumulates (see Figure 3) and reaches its maximum OO-AP1roG 0.966(−0.004) 82.2(−7.5) 23013(+192)
value of ln2 in the large U limit (for U → ∞, the DMRG[45] 0.970 89.7 22821
single-orbital density matrix has the diagonal elements
{0,0,0.5,0.5}). Note that OO-AP1roG yields similar
pair entanglement profiles for all chain lengths studied reference, the DMRG potential energy curve determined
and correctly reproduces the small and large U limits. inRef.45wasused,whichcanbeconsideredastheexact-
(b)Symmetric dissociation of the H molecule. diagonalization limit. None of the standard quantum
50
Thenon-relativisticquantumchemicalHamiltonianinits chemical methods, like MP2, CCSD(T) or DFT using
second quantized form reads the PBE exchange–correlation functional, yield qualita-
tively correct energy curves for the symmetric stretch-
Hˆ = (cid:88)hpqa†pσaqσ+12 (cid:88) (cid:104)pq|rs(cid:105)a†pσa†qτasτarσ+Hnuc, idnegptohf dtheeteHrm50incehdafinro.mInDpFaTrtiacnudlaMr,Pth2eisptootoendteiaelpe,nwehrgilye
pq,σ pqrs,στ
(4) CCSD(T) does not converge for interatomic distances
where the first term comprises the kinetic energy and larger than 2.0 ˚A. Note that the lack of size-consistency
nuclear–electron attraction, the second term is the in AP1roG is cured by orbital optimization in the OO-
electron-electron interaction, and the third term repre- AP1roG approach. The latter yields a potential energy
sents the nuclear–nuclear repulsion energy, respectively. curve that is closest to the DMRG reference data along
In Eq. (4), indices p, q, r and s run over all one-particle the whole dissociation pathway and leads to a proper
basis functions, while σ and τ denote the electron spin dissociation limit of H50. Moreover, the OO-AP1roG
({↑,↓}). The Hamiltonian as defined in Eq. (4) was method gives spectroscopic constants (presented in Ta-
used for the study of the symmetric stretching of the ble I) that are in excellent agreement with DMRG ref-
H hydrogen chain, which is a commonly-used molecu- erencedata,outperformingstandardquantum-chemistry
50
lar model for strongly-correlated systems and which re- approaches.
mains a challenging problem for conventional quantum- Wavefunctions constructed as antisymmetric products
chemistry methods [45–48]. of nonorthogonal geminals, like the AP1roG wavefunc-
In Figure 4, the performance of AP1roG and OO- tion scrutinized here, provide an alternative approach to
AP1roG is compared to restricted Hartree–Fock (RHF), electronic structure, with mean-field scaling. Because
second-order Møller-Plesset (MP2) perturbation theory, it uses electron pairs as a building block, AP1roG is a
coupled cluster theory with singles, doubles and pertur- suitable way to describe strong correlations dominated
bative triples (CCSD(T)), and density functional theory by electron pairing. However, in order to ensure size-
using the PBE [50] exchange–correlation functional. As consistency, the single-particle (orbital) basis used to
4
construct the electron pairs must be optimized. Our S.D.B. and D.V.N acknowledge financial support from
results show that orbital-optimized AP1roG is a robust FWO-Flanders and the Research Council of Ghent Uni-
method all the way from the weakly-correlated to the versity. S.D.BisanFWOpostdoctoralfellow. Wethank
strongly-correlated limit, in both molecules and periodic Ward Poelmans, Brecht Verstichel and Paul Johnson for
systems. results for the Hubbard model obtained by solving the
Acknowledgments. K.B. acknowledges the finan- Lieb-Wu equations. We had many helpful discussions
cial support from the Swiss National Science Founda- on geminals and exactly-solvable models with Peter Li-
tion (P2EZP2 148650). P.T. and P.W.A gratefully ac- macher and Paul Johnson. B.V., H.V.A., W.P., S.W.
knowledge financial support from the Natural Sciences and D.V.N. are Members of the QCMM alliance Ghent-
and Engineering Research Council of Canada. P.B. with Brussels.
[1] E. Dagotto, Rev. Mod. Phys. 66, 763 (1994). [26] J. V. Ortiz, B. Weiner, and Y. Ohrn, Int. J. Quantum
[2] K. Capelle and V. L. Campo, Phys. Report. 528, 91 Chem. S15, 113 (1981).
(2013). [27] P. R. Surjan, in Correlation and Localization (Springer,
[3] P.-O.Lo¨wdin,inAdv.Chem.Phys.,Vol.I(Wiley&Sons, 1999) pp. 63–88.
Inc,1958)Chap.Reviewofdifferentapproachesanddis- [28] W. Kutzelnigg, Chem. Phys. 401, 119 (2012).
cussion of some current ideas, pp. 209–321. [29] P.R.Surja´n,A.Szabados,P.Jeszenszki, andT.Zoboki,
[4] R. J. Bartlett and M. Musial(cid:32), Rev. Mod. Phys. 79, 291 J. Math. Chem. 50, 534 (2012).
(2007). [30] J. K. Ellis, R. L. Martin, and G. E. Scuseria, J. Chem.
[5] W.Foulkes,L.Mitas,R.Needs, andG.Rajagopal,Rev. Theory Comput. 9, 2857 (2013).
Mod. Phys. 73, 33 (2001). [31] P. A. Limacher, P. W. Ayers, P. A. Johnson, S. De
[6] D.I.Lyakh,M.Musia(cid:32)l,V.F.Lotrich, and.J.Bartlett, Baerdemacker, D.VanNeck, andP.Bultinck,J.Chem.
Chem. Rev. 112, 182 (2012). Theory Comput. 9, 1394 (2013).
[7] G. K.-L. Chan and S. Sharma, Annu. Rev. Phys. Chem. [32] P.Limacher, P.Ayers, P.Johnson, S.DeBaerdemacker,
62, 465 (2011). D.VanNeck, andP.Bultinck,Phys.Chem.Chem.Phys
[8] L. Pollet, Rep. Prog. Phys. 75, 094501 (2012). 16, 5061 (2014).
[9] R. Rodr´ıguez-Guzma´n, C. A. Jim´enez-Hoyos, R. Schut- [33] P.A.Limacher,T.D.Kim,P.W.Ayers,P.A.Johnson,
ski, andG.E.Scuseria,Phys.Rev.B87,235129(2013). S.DeBaerdemacker,D.VanNeck, andP.Bultinck,Mol.
[10] S. R. White, Phys. Rev. Lett. 69, 2863 (1992). Phys. , 853 (2014).
[11] U. Schollwo¨ck, Rev. Mod. Phys. 77, 259 (2005). [34] T. M. Henderson, J. Dukelsky, G. E. Scuseria, A. Sig-
[12] O. Legeza, R. M. Noack, J. So´lyom, and L. Tincani, noracci, andT.Duguet,arXivpreprintarXiv:1403.6818
in Computational Many-Particle Physics, Lect. Notes (2014).
Phys., Vol. 739, edited by H. Fehske, R. Schneider, and [35] T. Helgaker, P. Jørgensen, and J. Olsen, Molecular
A.Weiße(Springer,Berlin/Heidelerg,2008)pp.653–664. Electronic-Structure Theory (Wiley, 2000).
[13] K. H. Marti and M. Reiher, Z. Phys. Chem. 224, 583 [36] G.E.ScuseriaandH.F.SchaeferIII,Chem.Phys.Lett.
(2010). 142, 354 (1987).
[14] K.Boguslawski,K.H.Marti,O.Legeza, andM.Reiher, [37] U.Bozkaya,J.M.Turney,Y.Yamaguchi,H.F.Schaefer,
J. Chem. Theory Comput. 8, 1970 (2012). and C. D. Sherrill, J. Chem. Phys. 135, 104103 (2011).
[15] K. Boguslawski, P. Tecmer, O. Legeza, and M. Reiher, [38] Horton 1.2.0 2013, written by T. Verstraelen, S. Van-
J. Phys. Chem. Lett. 3, 3129 (2012). denbrande, M. Chan, F. H. Zadeh, C. Gonzalez, K. Bo-
[16] S. Wouters, P. A. Limacher, D. Van Neck, and P. W. guslawski, P. Tecmer, P. A. Limacher, A. Malek (see
Ayers, J. Chem. Phys. 136, 134110 (2012). http://theochem.github.com/horton/).
[17] Y. Kurashige, G. K.-L. Chan, and T. Yanai, Nature [39] E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445
Chem. 5, 660 (2013). (1968).
[18] K. Boguslawski, P. Tecmer, G. Barcza, O. Legeza, and [40] E. Neuscamman, Phys. Rev. Lett. 109, 203001 (2012).
M. Reiher, J. Chem. Theory Comput. 9, 2959 (2013). [41] L. Bytautas, T. M. Henderson, C. A. Jim´enez-Hoyos,
[19] P. Tecmer, K. Boguslawski, O. Legeza, and M. Reiher, J. K. Ellis, and G. E. Scuseria, J. Chem. Phys. 135,
Phys. Chem. Chem. Phys 16, 719 (2014). 044119 (2011).
[20] S. Wouters, W. Poelmans, P. W. Ayers, and [42] H. Bethe, Z. Phys. 71, 205 (1931).
D. Van Neck, Comput. Phys. Comm. XX, [43] E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445
DOI:10.1016/j.cpc.2014.01.019 (2014). (1968).
[21] M. C. Chung and I. Peschel, Phys. Rev. B 62, 4191 [44] S.-J. Gu, S.-S. Deng, Y.-Q. Li, and H.-Q. Lin, Phys.
(2000). Rev. Lett. 93, 086402 (2004).
[22] V. Murg, F. Verstraete, O. Legeza, and R. M. Noack, [45] J.Hachmann,W.Cardoen, andG.K.-L.Chan,J.Chem.
Phys. Rev. B 82, 205105 (2010). Phys. 125, 144101 (2006).
[23] A. C. Hurley, J. Lennard-Jones, and J. A. Pople, Proc. [46] T.TsuchimochiandG.E.Scuseria,J.Chem.Phys.131,
R. Soc. Lond. A 220, 446 (1953). 121102 (2009).
[24] A. J. Coleman, J. Math. Phys. 6, 1425 (1965). [47] L.Stella,C.Attaccalite,S.Sorella, andA.Rubio,Phys.
[25] D. M. Silver, J. Chem. Phys. 50, 5108 (1969). Rev. B 84, 245117 (2011).
5
[48] N. Lin, C. A. Marianetti, A. J. Millis, and D. R. Reich- [50] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.
man, Phys. Rev. Lett. 106, 096402 (2011). Lett. 77, 3865 (1996).
[49] W. J. Hehre, R. F. Stewart, and J. A. Pople, J. Chem.
Phys. 51, 2657 (1969).