Quantification of electron correlation effects – Quantum Information Theory versus Method of Increments

Quantification of electron correlation effects – Quantum Information Theory versus Method of Increments

Christian Stemmle christian.stemmle@fu-berlin.de    Beate Paulus Institut für Chemie und Biochemie - Takustr. 3, 14195 Berlin, Freie Universität Berlin, Germany
Abstract

Understanding electron correlation is crucial for developing new concepts in electronic structure theory, especially for strongly correlated electrons. We compare and apply two different approaches to quantify correlation contributions of orbitals: Quantum Information Theory (QIT) based on a Density Matrix Renormalization Group (DMRG) calculation and the Method of Increments (MoI). Although both approaches define very different correlation measures, we show that they exhibit very similar patterns when being applied to a polyacetelene model system. These results suggest one may deduce from one to the other, allowing the MoI to leverage from QIT results by screening correlation contributions with a cheap (“sloppy”) DMRG with a reduced number of block states.

\pdfstringdefDisableCommands

I Introduction

The objective of electronic structure theory is finding sufficient accurate solutions to the Schrödinger equation for chemical or physical applications. Unfortunately the correlated motion of this many-body problem turns out hard to be solved. An approach yielding numerically exact solutions is long know with the Full Configuration Interaction (FCI) method Löwdin (1957), usually based on a Hartree-Fock (HF) Hartree (1928); Fock (1930); Slater (1930) calculation which is missing electron correlation effects. However FCI is unfeasible for large systems due to factorial scaling with the number of electrons involved. The main challenge of electronic structure theory thus remains as finding approximate solutions, yielding results within a defined error margin.

Various approaches provide different trade-offs between accuracy and computational scaling. Density Functional Theory (DFT) Hohenberg and Kohn (1964) established itself for being capable of dealing with most systems of interest at reasonable computational cost, but has problems describing systems involving strong correlations and lacks systematic improvability. Additionally DFT requires some a priori knowledge of data, to validate whether a certain functional is suitable for the system in question. In other words, although DFT calculations may achieve high precision (i.e. have small statistical error), they may lack accuracy (i.e. have large systematic error), making them questionable for predictions.

Wave function based methods on the other hand are systematic truncations to the FCI problem, allowing them to treat any system with any desired accuracy and precision at the cost of higher computational effort. It is thus highly desirable to find superior trade-offs, making large system accessible for wave function based methods. Coupled-Cluster theory Bartlett (1981) evolved to be the standard method for single-reference systems, i.e. systems where the major part of the electron-electron interactions are accounted for by HF and thus most electron correlation effects are sufficiently describe by including dynamical correlations, by creating excitations based on the HF solution. Systems inadequately described by such single-reference approaches are called strongly correlated and require multi-reference treatments. A common approach for describing strong correlations are multi-configurational methods, like Multi Configuration Self Consistent Field (MCSCF) and its variant Complete Active Space Self Consistent Field (CASSCF) Roos et al. (1980); Szalay et al. (2012), including selected configurations of same or similar contributions to the wave function, while the HF method only considers a single configuration. Dynamic correlation can then be added in a second step by generating excitations based on many reference configurations, hence the name multi-reference methods.

The previously mentioned active space in CASSCF is a selection of a small subset of orbitals, which are deemed most relevant for static correlation. Hence this method still suffers from the same unfortunate scaling as the FCI problem, only less severe since usually the orbital space is much smaller. For cases where larger active spaces are required the Density Matrix Renormalization Group (DMRG) White (1992) is a suitable alternative. In DMRG the diagonalization of a large Hamiltonian is sidestepped by chunking it into many smaller matrices and keeping only the most important contributions. Another approach primarily developed for periodic and extended systems is the Method of Increments (MoI) Stoll (1992a, b), where the correlation energy is expanded in terms of groups of occupied orbitals (called centers). Both, DMRG and MoI, can benefit from localized orbitals to reduce long range correlations, which leads to faster convergence and reduced computational cost.

Another interesting aspect common to DMRG and MoI is the possibility to analyze and quantify the origin of correlation effects. While in the MoI the obtained increments can be directly interpreted, we can apply Quantum Information Theory (QIT) Horodecki et al. (2009); Szalay et al. (2015) to extract this data from the DMRG wave function (or in general any correlated wave function). Although there are technical differences between both, as will be explained in section II, we will show in section IV that we indeed find very similar trends. The model and computational details are presented in section III. Conclusion follows in section V.

Ii Theory

ii.1 Method of Increments

In the Method of Increments (MoI) Stoll (1992a, b); Paulus (2006); Voloshina and Paulus (2014); Fertitta et al. (2015) one partitions the correlation energy into contributions by assigning the occupied orbitals to different centers. These centers are correlated with the rest of the system (virtual orbitals) at different levels. The first level are the one-center or one-body increments

(1)

where is the reference energy of the uncorrelated system (Hartree-Fock), and the energy obtained by correlating center only, using a size extensive correlation method chosen by the user. Further levels are calculated by combining two, three or more centers and obtaining their total energies , and so on. Again, only the correlation energies

(2)
(3)

are considered. To avoid double accounting for correlations from lower levels, these have to be subtracted from the higher levels to obtain the corresponding increments

(4)
(5)

The correlation energy of the whole system can then be expanded in a series of increments

(6)

The number of increments at each level increases combinatorially making the method expensive at higher levels. But as higher level contributions converge to zero the expansion can be truncated, usually after the third level. Furthermore, using localized occupied orbitals allows neglecting contributions for pairs (or groups) of orbitals which are spatially separated, reducing the number of increments to be calculated further.

The individual correlation energies , , , etc. may also be interpreted directly as the correlation effects arising from each center, while increments like or measure the additional effect of correlating the combined group.

ii.2 Density Matrix Renormalization Group (DMRG) and Quantum Information Theory (QIT)

The Density Matrix Renormalization Group (DMRG) method was first invented by White (1992) and is suitable to treat strong correlations in large active spaces and is closely connected to Quantum Information Theory (QIT) Horodecki et al. (2009) which allows for quantification of correlation effects in terms of (groups of) orbitals. As both, DMRG and QIT, are described in detail in various reviews Legeza et al. (2008); Marti and Reiher (2010); Kurashige et al. (2013); Wouters et al. (2014); Szalay et al. (2015); Olivares-Amaya et al. (2015); Wouters et al. (2016), we restrict ourselves to a brief summary here.

In the Density Matrix Renormalization Group (DMRG) method the Full Configuration Interaction (FCI) wave function is approximated by trying to find its most important contribution. In essence the complete FCI Hamiltonian is not diagonalized at once, but optimized iteratively in subspaces of e.g. 2 orbitals. During this iterative scheme the most important contributions are kept and carried over to the next iteration. This procedure is facilitated by storing the wave function in a Matrix Product State (MPS), whose accuracy can be controlled by the dimensions connecting two matrices (called number of block states or virtual dimensions). Both, the iterative diagonalization of smaller subspaces and efficient storing of the wave function in MPS format, allow for treatment of large active spaces and is especially suitable for treating static (or strong) correlations.

For the MPS the orbitals are thought of to be arranged in a linear chain, of arbitrary order. The Configuration Interaction (CI) wave function is then given by

(7)

where indicates the orbitals position on the chain of length . The indices label the single-orbital basis states , which correspond to the spin occupations of a spatial orbital: , , and . The CI coefficients are given by the -order tensor and is capable of storing coefficients for any configuration with an electron count from 0 to . Thus its memory requirement grows exponentially as .

The Matrix Product State (MPS) facilitates these memory requirements by factorizing the tensor as a product of low order tensors and with controlled rank

(8)

where each matrix corresponds to one molecular orbital . The factorization can formally be obtained by succeeding application of singular value decompositions and is exact, i.e. the full tensor can be recovered without loss of information. As a result the size of the matrices is still growing exponentially towards the center of the chain Schollwöck (2005). For the MPS to reduce the memory requirements, one approximates the matrices by defining an upper limit to the matrix dimensions called number of block states or virtual bond dimensions . In practice the challenge of DMRG is now to find an appropriate value of and a suitable order of orbitals in the chain which will lead to an optimized set of matrices to describe the CI wave function.

QIT is closely connected to DMRG, as it allows for quantification of orbital correlation contributions (entanglement) based on the CI wave function. In turn the QIT results may be used to find better parameters for the DMRG calculation. One can easily calculate the -orbital reduced density matrices, by contracting the MPS over all but orbitals, e.g. the one-orbital reduced density matrix is given by

(9)
(10)

The correlation contribution of a single orbital are then given by the one-orbital von Neumann entropy Legeza and Sólyom (2003)

(11)

where are the eigenvalues of the one-orbital reduced density matrix . The above equation will give small values if all main configurations of the CI wave function have the same occupation in orbital , i.e. one of the will be close to while the other three are negligible. Largest values are obtained if the occupation in orbital for all main configurations is equally distributed over all (, , , ), i.e. when for all , then .

The sum of all one-orbital entropies gives a measure for the total correlation

(12)

of the wave function Legeza and Sólyom (2004); Szalay et al. (2017).

Analogues to the one-orbital von Neumann entropy, higher orders can be calculated from their corresponding -orbital reduced density matrix. For example the two-orbital von Neumann entropy is obtained from the two-orbital density matrix Legeza and Sólyom (2006)

(13)
(14)

and quantifies the correlation contributions of the combined two-orbital subsystem . To quantify correlations between and the mutual information Rissler et al. (2006) has to be calculated

(15)

ii.3 Comparison of Increments and Entropies

Both, MoI and QIT, allow for quantification of correlation contributions by single orbitals or groups of them. If we assign each orbital to its own center in the MoI, we can directly compare the increments (, ) with the QIT quantities (, ). To stress this assignment, we will switch terminology from here on, and call the one- and two-center increments, one- and two-orbital increments respectively.

Indeed a close similarity is obvious when comparing eq. 4 and eq. 15. The two equations suggest a close connection of the two-orbital increment with the mutual information , as well as of the one-orbital increments with the one-orbital entropy and the two-orbital correlation energy with the two-orbital entropy . Therefore, one might ask whether it is possible to infer from the QIT quantities to the increments and vice versa?

It should be pointed out however, that there are some main differences between entropies and increments. First of all, in the MoI one only accounts for correlation effects of the current increment with the rest of the system. For example when calculating the two-orbital increments the correlations effects of all other orbitals are not considered, i.e. each increment is only correlated with the virtual orbitals, but not the other increments on its level. On the other hand, the QIT entropies rely on a CI wave function which correlates all orbitals at once.

The second main difference is the definition of both quantifies. The increments are essentially energies obtained by applying the Hamiltonian to different truncations of the CI wave function, while for the entropies we first construct reduced density matrices (all from the same CI wave function), and then apply a logarithmic function (cf. eq. 11) to its eigenvalues, making it a non-linear mapping.

Thus there is no direct correspondence of increments and entropies, and quantitative differences are to be expected. We will show however, that qualitative agreement can indeed be observed.

Iii Model and Computational Details

As a basic molecular model system we choose conjugated trans-polyacetelenes, whose single and double bonds will provide varying degrees of correlation effects. MoI calculations have been successfully applied to such systems before, and allowed to discriminate between different bonding situations Yu et al. (1997). Additionally, linear systems are better suited for the linear MPS structure in DMRG calculations. For our study we choose the hexatriene (\ceC6H8), which is just large enough to allow for spatially well separated localized orbitals. Furthermore it provides two different kinds of double bonds, one at the center of the chain, the other two on each end of the chain. The used bond distances are for \ceC-C bonds, for \ceC=C bonds and for \ceC-H bonds, while all bond angles are set to .

Hartree-Fock orbitals are obtained using the cc-pVTZ basis set Kendall et al. (1992) and occupied orbital localized using Pipek-Mezey Pipek and Mezey (1989) localization. For the active space we exclude the 6 carbon core orbitals, and include all orbitals which can be constructed from the carbon and hydrogen shells, resulting in 16 occupied and 16 virtual orbitals. We thus use for both, DMRG and MoI, an active space with 32 electrons in 32 orbitals (CAS(32,32)). This is similar to previously reported CAS-MoI and DMRG calculations on \ceBe6 ring Fertitta et al. (2015).

For the DMRG calculations the Budapest DMRG program Legeza (2016) was used, while all other calculations were performed using Molpro Werner et al. (2015, 2012). For the DMRG and QIT results presented in table 1 and fig. 8 the Dynamic Block State Selection (DBSS) approachLegeza et al. (2003); Legeza and Sólyom (2004) was used, with an density matrix cutoff of and a maximum number of block states . Together with an optimized orbital ordering these parameters provide high accuracy for the correlation energy. For the QIT quantities however, it is sufficient to use provided an appropriate orbital ordering is used.

The QIT results will give us entanglement measures among occupied and virtual orbitals, as well as in between them in just one run. However, in the conventional MoI approach, we can only get the correlations for (groups of) occupied centers. In order to compare both measures (QIT and MoI) directly we define each occupied orbital to be its own center. Furthermore, to access the correlation measures for virtual orbitals, we can flip things around and expand the increments in virtual orbitals and correlate these centers with all occupied orbitals. For large virtual spaces, this will drastically reduce the active space for each individual increment calculation and in turn increase the number of centers. This approach has recently also been suggested by Eriksen et al. (2017) in an effort to improve parallelism of such calculations and make large orbital spaces more accessible.

Iv Results

Total Energy () Correlation Energy ()
HF -231.883732
DMRG(32,32) -231.971288 -0.087556
CCSD -232.891987 -1.008255
Table 1: Total energy and correlation energy for trans-hexatriene obtained by with different methods using a cc-pVTZ basis set.
Level Correlation Energy () Summed Correlation Energy ()
CAS(32,32)-MoI (occupied) 1 -0.053475 -0.053475
2 -0.035636 -0.089111
3 0.000988 -0.088122
4 0.000393 -0.087729
CAS(32,32)-MoI (virtual) 1 -0.022884 -0.022884
2 -0.063367 -0.086251
3 -0.003104 -0.089355
4 0.001568 -0.087787
CCSD-MoI (occupied) 1 -0.511070 -0.511070
2 -0.552483 -1.063553
3 0.060588 -1.002969
Table 2: Correlation energies obtained with various MoI variants for trans-hexatriene using a cc-pVTZ basis set. Occupied orbitals have been localized using the Pipek-Mezey method.

We start our discussion with the total energies calculated by different methods. The energies in table 1 provide references for the uncorrelated system (HF), the static correlation contribution (DMRG) and dynamic correlations (CCSD). In table 2 we present the MoI results at different levels of increments, again dealing with static (CAS(32,32)-MoI) and dynamic correlation (CCSD-MoI). The 3-orbital increments level yields in all cases agreement within , while with 4-orbital increments reduce the difference by one order of magnitude to . In case of the CAS-MoI, we also expanded the correlation energy in terms of virtual orbitals, yielding similar results when including the 3- and 4-orbital increments.

iv.1 Static Correlation

(a)
(b)
(c)
(d)
(e)
(f)
Figure 7: Isosurface plot () for the localized occupied molecular orbitals #12 to #16 and the canonical virtual molecular orbital #17.
Figure 8: Comparison of Increments and Orbital Entropies for quantifying orbitals correlations. Left and right column show results based on Methods of Increments (MoI) and Quantum Information Theory (QIT) respectively. The upper and middle row show 1-orbital and 2-orbital correlations respectively. The Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) are separated by black and white lines respectively. In the lower row the 2-orbital correlations among occupied orbitals is shown as a zoom in. MoI increments are based on CASCI(32,32)/cc-pVTZ calculations.

Before investigating the increments and entropies we take a brief look at a selected set of the molecular orbitals fig. 7. For a complete list of the active space orbitals please refer to the Supplementary Information (Figs. S1 to S4). We will focus our discussion to the exemplary two-orbital correlations involving localized occupied molecular orbitals describing \ceC-H bonds (#12 and #13) and \ceC=C bonds (#14 to #16) together with the Lowest Unoccupied Molecular Orbital (LUMO) #17. The remaining occupied orbitals describe the \ceC-C bonds (#1 to #5) and further \ceC-H bonds (#6 to #11). Note that \ceC-H orbitals #12 and #13 as well as \ceC=C orbitals #15 and #16 are localized at the ends of the linear molecule, which results in examples for spatially close lying and well separated two-orbital correlations.

In fig. 8 we present the one- and two-orbital increments (left column) as well as one-orbital entropy and mutual information (right column). Note that the two-orbital increments are calculated by doing two separate incremental expansions: the first one in terms of occupied orbitals, the second one in terms of virtual orbitals. The results are however compiled in a single picture. This cut between occupied and virtual orbitals is indicated by black and white lines for the one- and two-orbital quantifies respectively. The off-diagonal occupied-virtual pairs for the two-orbital increments cannot be calculated by the MoI method, therefore these missing values are indicated in white.

Comparing both columns we find very similar patterns, but with varying amplitudes. It is evident from both the one-orbital increments and one-orbital entropy (upper row in fig. 8), that the orbitals with largest correlation contributions are #14 to #17, representing the bonds. Out of these four, the one-orbital increments identify orbitals #15 and #16, the bond on each end of the chain, to be the most important one, and the LUMO #17 to be the least important one. However, the 1-orbital entropies shows a reversed trend, and even includes the virtual orbital #21 as the second highest one. Both approaches agree though in the trend for the bond orbitals: orbitals #1 to #5 (\ceC-C) are smaller than orbitals #6 to #13 (\ceC-H), being much smaller than the orbitals.

The two-orbital correlations (middle row in fig. 8) exhibit similar patterns as well, again with varying amplitudes. Note however, that the largest contributions identified by the mutual information, namely orbital pairs #17#14, #17#15 and #17#16, are not included in the MoI as they include one occupied and one virtual orbital. Focusing on the occupied correlations only (bottom row in fig. 8), we see that both methods again agree in the general trend: most important correlation contributions arise from the three bonding orbitals, while the spatially separated pair #15#16 is negligible. All the degenerate pairs of bond orbitals (e.g. #12 and #13) show identical correlation contributions when paired with the centered orbital #14, as their distance is the same. In connection with the orbitals on the ends of the chain (#15, #16) an alternating structure of large and small correlations emerges from the varying spatial separation.

iv.2 Dynamic Correlations

The inclusion of dynamical correlations by increasing the virtual orbital space leads costly calculations for the DMRG based QIT results. The MoI however, can easily be based on the Coupled Cluster Singles Doubles (CCSD) method. The latter will be identical to Full Configuration Interaction (FCI) calculations for 1-orbital increments, as only two electrons are correlated in each individual calculation. Similar, for the 2-orbital increments with 4 electrons each, we can expect the CCSD error to be negligible. Furthermore, increments expanded in virtual orbitals will be the same as for the CAS-MoI results above, as the occupied orbital space remains unchanged. It will however add contributions of higher lying virtual orbitals.

CCSD-MoI results expanded in terms of occupied orbitals are shown in fig. 9. We observe a large energy shift for all increments, as the virtual orbital space increased drastically. This indicates dynamical correlations being similar important for all occupied orbitals. However, the effect is less pronounced for the orbitals (#14 to #16), resulting in similar values as for the \ceC-H bond orbitals (#6 to #12). Only the \ceC-C sigma bond orbitals (#1 to #5) remain distinguishable by their value.

The 2-orbitals increments exhibit a very similar pattern, but quite large differences in the magnitude. Although being the smallest 1-orbital increments, the orbitals #1 to #3 now have largest contributions in the 2-orbital increments. These 2-orbital increments correspond to the combinations of \ceC-C and bond, being localized on the same bonds (#1#14 representing the central \ceC-C bond; #2#15 and #3#16 the two outer most \ceC-C bonds; cf. Supplementary Information). We can therefore conclude, that smaller lower-level increments do not necessarily indicated negligible higher-level contributions, thus care must be taken when trying to deduce from one to the other. Furthermore we can observe that the 2-orbital increments for the two neighboring pairs is now one of the smallest contribution (but still with increased value compared to the CAS-MoI results). The colors for \ceC-C and \ceC-C are essentially reversed for static and dynamical correlation effects. Thus we can see that correlation effects between spatially overlapping and bond are mainly of dynamical nature, in contrast to the static correlations between neighboring, degenerate orbitals.

Figure 9: Occupied orbital increments including dynamical correlations based on CCSD/cc-pVTZ calculations.

V Summary and Discussion

We have performed calculation recovering static and dynamical correlations by applying the MoI and DMRG. By comparing the individual contribution (increments vs. entropies) we can see, that both show very similar patterns in terms of which contributions are most important. The actual values however, show large difference, as was expected due to the different nature of increments and entropies. Deducing from one set of results to the other thus seams an appropriate approach, care must be taken though when deciding on which cutoff parameter to use when neglecting contributions.

Comparing the computational effort for both methods we see that MoI is rather restricted. Each individual increment requires a separate calculation, and the number of increments increases combinatorially with the number of orbitals (centers) and level. Furthermore, MoI it is not capable of showing cross-correlations between occupied and virtual orbitals. This can be problematic when trying to determine a smaller appropriate active space. Here QIT has the advantage of providing all relevant information, just based on a single CI wave function. Additionally, QIT can be easily extended to consider contributions from combining 3, 4 or more groups since the wave function is readily available. As a further benefit the QIT entropies are insensitive to the quality of the DMRG calculation (number of block states). Indeed the one-orbital entropies and mutual information presented here (cf. fig. 8) can be obtained with same quality by using only block states (cf. Fig. S5 of the Supplementary Information). The time and memory limiting factor is then likely to be the construction and diagonalization of the reduced density matrices.

As a possible application we may first screen all static correlation effects by applying DMRG on a large active space. In a next step we can then try to select the increments for a MoI calculation, able to obtain a more accurate correlation energy and recover additional dynamic correlations, which DMRG is not suited for.

Acknowledgements.
We would like to thank Prof. Peter Fulde (Dresden) for discussion and stimulating this topic. Financial support from the International Max Planck Research School “Functional Interfaces in Physics and Chemistry” is gratefully acknowledged. The high performance computing facilities of the Freie Universität Berlin (ZEDAT) are acknowledged for computing time.

References

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
278651
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description