Efficient and accurate treatment of electron correlations from first-principles

Efficient and accurate treatment of electron correlations from first-principles

Y. X. Yao Ames Laboratory–US DOE and Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA    J. Liu Ames Laboratory–US DOE and Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA    C. Liu Ames Laboratory–US DOE and Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA    W. C. Lu State Key Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun, Jilin 130021, China College of Physical Science and Laboratory of Fiber Materials and Modern Textile, Growing Base for State Key Laboratory, Qingdao University, Qingdao, Shandong 266071, China    C. Z. Wang Ames Laboratory–US DOE and Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA    K. M. Ho Ames Laboratory–US DOE and Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA
Abstract

We present an efficient ab initio method for calculating the electronic structure and total energy of strongly correlated electron systems. The method extends the traditional Gutzwiller approximation for one-particle operators to the evaluation of the expectation values of two particle operators in a full many-electron Hamiltonian. The method is free of adjustable Coulomb parameters, and has no double counting issues in the calculation of total energy, and has the correct atomic limit. We demonstrate that the method describes well the bonding and dissociation behaviors of the hydrogen and nitrogen clusters. We also show that the method can satisfactorily tackle great challenging problems faced by the density functional theory recently discussed in the literature. The computational workload of our method is similar to the Hartree-Fock approach while the results are comparable to high-level quantum chemistry calculations.

pacs:
31.10.+z, 31.15.E-,31.15.vn, 71.15.Nc

It is one of the outstanding challenges in physics, chemistry, and materials science to develop robust and efficient theoretical and computational methods to accurately calculate the electronic structure and total energy of materials containing strongly correlated electrons  not (). While accurate methods are available from quantum chemistry approaches (e.g., configuration interaction (CI)), these methods are too expensive for condensed matter systems. On the other hand, density functional theory (DFT) and related computational codes based on the Kohn-Sham approach Hohenberg and Kohn (1964); Kohn and Sham (1965) have been well developed, and are highly effective and successful for predicting the structures and properties of many materials, but they fail for systems with strongly correlated electrons. In the last three decades, there have been intensive efforts in developing new approaches to solve the outstanding problems in correlated electron systems  Anisimov et al. (1991, 1997); Georges et al. (1996); Savrasov et al. (2001); Kotliar et al. (2006); Zgid and Chan (2011); Lin et al. (2011); Knizia and Chan (2012); Ho et al. (2008); Yao et al. (2011, 2012); Deng et al. (2008, 2009); Wang et al. (2010); Lanatà et al. (2012); Schickling et al. (2012); Lanatà et al. (2013). Among these new developments, local-density approximation plus on-site Coulomb interaction parameter U (LDA+U) Anisimov et al. (1991, 1997), LDA+dynamical mean field theory Georges et al. (1996); Savrasov et al. (2001); Kotliar et al. (2006), and LDA+Gutzwiller Ho et al. (2008); Yao et al. (2011); Deng et al. (2008, 2009); Wang et al. (2010); Lanatà et al. (2012); Schickling et al. (2012); Lanatà et al. (2013) have emerged as the most popular methods for treating strongly-correlated electrons in solid-state systems. These methods handle electron correlations through the adjustable on-site Coulomb interaction parameters, while keeping the full description of the electronic structure through LDA. From the aspect of the theoretical predictive power, it is highly desirable to have a fully self-consistent ab initio theory that can treat correlated electron systems without adjustable parameters and with computational speed comparable to LDA or Hartree-Fock (HF) calculations Gebauer et al. (2013).

In this letter, we present an ab initio method for the electronic structure and total energy calculations of strongly correlated electron systems without adjustable Coulomb parameters. In our approach, the commonly-adopted Gutzwiller approximation (GA) for evaluating the one particle density matrix Gutzwiller (1965); Kotliar and Ruckenstein (1986); Bünemann et al. (1998); Bünemann and Gebhard (2007) is extended to treat the evaluation of the two-particle correlation matrix of the system. This approximation, which we call the correlation matrix renormalization (CMR) approximation Yao et al. (2014), allows the expectation value of a many-electron Hamiltonian with respect to Gutzwiller variational wave function (GWF) to be evaluated with reduced computational complexity. We show that the method describes well the bonding and dissociation behaviors of hydrogen and nitrogen clusters in comparison with the accurate and expensive quantum chemistry calculations. Furthermore, some of the most challenging problems faced by Kohn-Sham DFT-based calculations recently discussed in the literature  Cohen and Mori-Sánchez (2014); Cohen () can also be readily solved by our method. The method has no double counting issues in the calculation of total energy, and produces the correct atomic limit. The computational workload scales as N or better for systems of size N, similar to the HF calculationsYao et al. (2014).

We start with the full ab initio Hamiltonian for an interacting many-electron system in the second quantization form

(1)

Here , , , are the atomic site indices. , , , are orbital indices and is the spin index. are eigenstates of the local many-body Hamiltonian. The first term is the local on-site part which has been singled out for exact treatment, is the energy of the local many-electron configuration . The second and third terms describe the non-local one-body and two-body contributions. All interactions are included in this Hamiltonian without any adjustable parameters. When evaluating this Hamiltonian with the full CI wave function, one obtains an exact expression of the total energy which consists of non-local one-particle and two-particle density matrices in addition to the local on-site contributions. In our CMR approach, we evaluate the Hamiltonian in Eq. 1 with the GWF of the form

(2)

Where is a non-interacting electron wavefunction, i.e., a single Slater determinant. is the variational parameter which determines the occupation probability of the on-site configuration . The central part of the Gutzwiller approach is the suppression of the energetically unfavorable atomic configurations in the many-body wave function. Using the GWF of Eq. 2 and adopting the generally accepted GA for the expectation value of a one-particle operator Bünemann et al. (1998); Bünemann and Gebhard (2007), the total energy of the system in our CMR scheme can be expressed as

(3)

where is the occupation weight of the configuration .  can be evaluated following the standard GA rule for one-particle hopping operators, i.e., if and otherwise. To reach the expression Eq. 3, the validity of Wick’s theorem has been assumed. is the residual correlation energy due to the approximation involved in the CMR approach. In general, can be determined by comparing the total energies from the CMR with that from accurate CI or quantum Monte Carlo calculations for some exactly solvable structures. Since the dominant local onsite electron correlation effect has been taken into account by the GWF, the residual correlation energy due to the CMR approximation is expected to be small. In the test cases to be shown in this letter, we find that one way to include the effects of is to modify the renormalization -factor obtained from the GA.

Figure 1: (Color online) The total energy (upper panel) and double occupancy weight (lower panel) of H clusters as a function of bond length calculated from the CMR method, which agree very well with the results from the exact CI calculations. The CMR results are also much better than the HF results. The calculations are done using a minimal basis set.

We first demonstrate the CMR method by studying the dissociation behavior of the hydrogen molecules. The dissociation behavior of these hydrogen molecules has been the testing ground of methods for correlated electron calculations, because the electron correlation changes from the weak to strong regime as the hydrogen bond length increases. For these systems, the residual correlation energy is included by replacing the renormalization -factor obtained using the density-density type GA Bünemann et al. (1998) by the functional of , i.e., . The is determined by requiring that the total energy (without explicit term) and the probability of the local double occupancy for a H dimer obtained from the CMR calculation to be the same as the exact CI results. We first performed the test using the minimal basis set (one orbital for each H atom). In this case, the total energy and double occupancy probability from the CI calculation can be solved analytically and the has an analytical form in term of . Using the analytically determined from the reference system H, CMR calculations are performed for H-ring, H-ring, H-ring and H-cube structures. The results from our CMR calculations are presented in Fig. 1 in comparison with the full CI or the highly accurate multi-configurational self-consistent field (MCSCF) results. We found the bonding and dissociation behavior of the hydrogen clusters calculated from the CMR method agrees very well with the result from the high level quantum chemistry calculations. In contrast, the HF results show large systematic errors, especially at large separations where the electron correlation effect becomes prominent, as evidenced by the strong suppression of the energetically unfavorable local electron double occupancy weight.

We further tested the CMR method for the dissociation behavior of hydrogen clusters using a large basis set of 6-311G**, which contains 3 -orbitals plus 3 -orbitals. In this case, the needs to be determined numerically by fitting the CMR energies and the local configuration occupation probabilities of the H dimer to the exact CI results. Using such a numerically constructed functional , we have performed the CMR calculations for H-ring, H-ring, and H-cube with the same large basis set. In Fig. 2 we show that the CMR method yields again very good bonding and dissociation curves in close agreement with the MCSCF calculations. The inset of Fig. 2 shows the behavior of , which scales like at small and approaches as goes to .

Figure 2: (Color online) The total energy of hydrogen clusters as a function of bond length calculated from the CMR method, which agree well with the results from the high-level quantum chemistry CI or MCSCF calculations. The HF results are also shown for comparison. The calculations are done using large basis set. Insets: obtained by fitting the CMR total energy and local configuration weights with the exact full CI results.

Very recently, Cohen, et al used some prototype systems to show the dramatic errors in the DFT-based calculations. These errors stem from the fact that the current approximations used in DFT calculations miss the energy derivative discontinuity with respect to the total electron numberPerdew et al. (1982); Cohen and Mori-Sánchez (2014); Cohen (). The prototype systems to reveal the failure of DFT are the stretched few-electron systems, e.g., one-electron systems like HZ and HZH and two-electron systems like HZ with Z being the proton with nucleus charge varying between 0 and 2. While the electron density from the exact calculations shows dramatic discontinuous changes in real space with a slight variation of near some critical points at large separations, all the DFT calculations predict an artificial continuous variation of the electron density Cohen and Mori-Sánchez (2014). Remarkably, our CMR method gives exact solutions for any single-electron systems, as it can be easily proved that the orbital renormalization factors are constantly one and the method reproduces the HF results, which are exact in the special class of one-electron systems. It is also remarkable that our CMR method yields the exact bonding and dissociation behaviors for both H and H (see Fig. 2), while all the available DFT calculations can hardly describe both cases equally well Cohen (). One can further show that because the CMR method reaches the correct atomic solutions at the large separation limit, the exact discontinuous electron transfer observed in the HZ system at large separations can be well reproduced. In Fig. 3 we compare the electron occupation and double occupancy weight of Z atom from the CMR, HF, DFT with the generalized gradient approximation (GGA) and CI. Although all the methods predict similar results near equilibrium bond length (0.75Å), the CMR method shows significant improvements over the mean field HF and GGA and follows closely the exact CI results with increasing separations, even at the chemically crucial bond breaking region (2Å) and beyond. The underlying physics for the large errors of the simple mean field approaches like the HF and GGA can be understood by noting that the mean field double occupancy weight evaluated using the CI orbital occupation, shown as the dotted line in the lower panel of Fig. 3, can severely deviate from the exact CI double occupancy weight—manifesting the multi-configuration nature of the exact solution which is beyond the single Slater determinant description.

Figure 3: (Color online) The electron occupation number (upper panel) and the local double occupancy weight  (lower panel) as a function of the nucleus charge on the Z atom calculated with the CMR, HF, GGA and CI methods for the HZ dimer at three seperations: near equilibrium (left), close to bond breaking (middle) and beyond (right). The dotted line in the lower panel shows the mean field double occupancy weight evaluated at the CI electron occupation of the Z atom.

Another challenging prototype system is the H cluster with varying electron filling Cohen (). The exact solution predicts a relatively big energy gap for the system at large separations and half-filling ; while all the DFT calculations fail to reproduce this result because of the incapability to treat the strong electron correlation effects. In Fig. 4 we show the total energy of the H cube from the CMR, HF, GGA and MCSCF calculations as a function of even number of electron filling, which keeps the system to have the closed shell ground state solution Cohen (). While all the four theories give similar total energies at the small bond length, the discrepancy between them becomes increasingly large with expanding the H cluster. Remarkably, the CMR energies agree with the highly accurate MCSCF results very well for all the bond separations and electron fillings, which proves that the key many-body correlation physics in this system has been perfectly captured by the CMR method. A better comparison between the four levels of theories is presented by the energy gap, defined as the second order finite difference , as shown in the insets of Fig. 4. Clearly, as the bond length increases, or the electron correlation effects become stronger, the simple mean field HF and GGA energy gap shows larger deviations from the exact gap, especially the gap at half filling. In contrast, the CMR calculations yield energy gaps in excellent agreement with the MCSCF results in all the cases.

Figure 4: (Color online) The total energy as a function of total number of electrons obtained from the CMR, HF, GGA and MCSCF calculations for the H cube with increasing bond length . Insets: the corresponding second order energy difference .

The CMR method is also successfully applied to systems with atoms containing multiple correlated orbitals, e.g., nitrogen clusters. For computational convenience, we describe the nitrogen atom with the minimum basis set and choose the and as the correlated orbitals. The same idea can be equally well carried over to the large basis calculations as shown previously for the hydrogen clusters. Two functionals, and , are introduced to modify the renormalization coefficients of and orbitals. The specific functional forms, following the procedure in the calculations of hydrogen clusters, are determined by matching the CMR total energy, , and local correlated Fock state occupation probabilities, , with the exact CI results of the N dimer. We apply the method to calculate binding energy curves of three nitrogen clusters of different geometries, i.e., the square, diamond and tetragonal shapes. In Fig. 5 we show the total energy as a function of bond length from the CMR, HF and MCSCF calculations. The good agreement between the CMR and MCSCF energies for all the structures demonstrates the good transferability of our method.

Figure 5: (Color online) The total energy of the nitrogen clusters as a function of bond length calculated from the CMR method, which agrees well with the results from the high-level CI or MCSCF calculations. The HF results are also shown for comparison.

In summary, we have developed an efficient method for calculating the electronic structure and total energy of the systems with strong electron correlations. The method is based on the Gutzwiller type variational wavefunction and adopts a correlation matrix renormalization approximation in which both one-particle density and two-particle correlation matrices at mean field level are renormalized according to the local electron correlation effects. While the computation workload of this new approach is similar to that in HF calculations, the calculation results are much more accurate. The benchmark results for the bonding and dissociation behaviors of the hydrogen and nitrogen clusters show that our method well reproduces the results from the accurate and yet expensive quantum chemistry CI and MCSCF calculations. The CMR method is also demonstrated to be accurate for treating the electron correlation effects in some prototype systems where the current DFT and HF calculations fail. The extension of the method to crystalline solids is straightforward and promising. The work along this direction is underway.

Acknowledgements.
This work was supported by the U. S. Department of Energy, Basic Energy Sciences, Division of Materials Science and Engineering under the Contract No. DE-AC02-07CH11358 including the computer time support from the National Energy Research Supercomputing Center (NERSC) in Berkeley, CA.

References

  • (1) See, for example, DOE-Office of Science Report: “Basic Research Needs for Advanced Nuclear Energy Systems”, 2006.
  • Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  • Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  • Anisimov et al. (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
  • Anisimov et al. (1997) V. I. Anisimov, F. Aryasetiawan, and A. Lichtenstein, J. Phys. Condens. Matter 9, 767 (1997).
  • Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  • Savrasov et al. (2001) S. Y. Savrasov, G. Kotliar, and E. Abrahams, Nature 410, 793 (2001).
  • Kotliar et al. (2006) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
  • Zgid and Chan (2011) D. Zgid and G. K.-L. Chan, The Journal of Chemical Physics 134, 094115 (2011).
  • Lin et al. (2011) N. Lin, C. A. Marianetti, A. J. Millis, and D. R. Reichman, Phys. Rev. Lett. 106, 096402 (2011).
  • Knizia and Chan (2012) G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
  • Ho et al. (2008) K. M. Ho, J. Schmalian, and C. Z. Wang, Phys. Rev. B 77, 073101 (2008).
  • Yao et al. (2011) Y. X. Yao, C. Z. Wang, and K. M. Ho, Phys. Rev. B 83, 245139 (2011).
  • Yao et al. (2012) Y. X. Yao, C. Z. Wang, and K. M. Ho, International Journal of Quantum Chemistry 112, 240 (2012).
  • Deng et al. (2008) X. Deng, X. Dai, and Z. Fang, EPL (Europhysics Letters) 83, 37008 (2008).
  • Deng et al. (2009) X. Y. Deng, L. Wang, X. Dai, and Z. Fang, Phys. Rev. B 79, 075114 (2009).
  • Wang et al. (2010) G. Wang, Y. Qian, G. Xu, X. Dai, and Z. Fang, Phys. Rev. Lett. 104, 047002 (2010).
  • Lanatà et al. (2012) N. Lanatà, H. U. R. Strand, X. Dai, and B. Hellsing, Phys. Rev. B 85, 035133 (2012).
  • Schickling et al. (2012) T. Schickling, F. Gebhard, J. Bünemann, L. Boeri, O. K. Andersen, and W. Weber, Phys. Rev. Lett. 108, 036406 (2012).
  • Lanatà et al. (2013) N. Lanatà, Y.-X. Yao, C.-Z. Wang, K.-M. Ho, J. Schmalian, K. Haule, and G. Kotliar, Phys. Rev. Lett. 111, 196801 (2013).
  • Gebauer et al. (2013) R. Gebauer, M. H. Cohen, and R. Car, ArXiv e-prints (2013), eprint 1309.3929.
  • Gutzwiller (1965) M. C. Gutzwiller, Phys. Rev. 137, A1726 (1965).
  • Kotliar and Ruckenstein (1986) G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57, 1362 (1986).
  • Bünemann et al. (1998) J. Bünemann, W. Weber, and F. Gebhard, Phys. Rev. B 57, 6896 (1998).
  • Bünemann and Gebhard (2007) J. Bünemann and F. Gebhard, Phys. Rev. B 76, 193104 (2007).
  • Yao et al. (2014) Y. X. Yao, J. Liu, C. Z. Wang, and K. M. Ho, Phys. Rev. B 89, 045131 (2014).
  • Cohen and Mori-Sánchez (2014) A. J. Cohen and P. Mori-Sánchez, The Journal of Chemical Physics 140, (2014).
  • (28) A. J. Cohen, http://meetings.aps.org/link/BAPS.2014. MAR.M1.1.
  • Perdew et al. (1982) J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys. Rev. Lett. 49, 1691 (1982).
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 ...
53360
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