# Explicit inclusion of electronic correlation effects in molecular dynamics

###### Abstract

We design a quantum molecular dynamics method for strongly correlated electron metals. The strong electronic correlation effects are treated within a real-space version of the Gutzwiller variational approximation (GA), which is suitable for the inhomogeneity inherent in the process of quantum molecular dynamics (MD) simulation. We also propose an efficient algorithm based on the second-moment approximation to the electronic density of states for the search of the optimal variation parameters, from which the effective renormalized interatomic MD potentials are fully determined. By considering a minimal one-correlated-orbital Anderson many-particle model based on tight-binding hopping integrals, this fast GA-MD method is benchmarked with that using exact diagonalization to solve the GA variational parameters. In the infinite damping limit, the efficiency and accuracy are illustrated. This novel method will open up an unprecedented opportunity enabling large-scale quantum MD simulations of strongly correlated electronic materials.

^{†}

^{†}preprint: submitted to Phys. Rev. B

## I Introduction

Electronic correlation effects in materials such as transition metal oxides, give rise to emergent phenomena including Mott insulator, magnetism, heavy fermion, and unconventional superconductivity. These phenomena defy the description of the density functional theory (DFT) within local density approximation (LDA), which has been successful in describing electronic and structural properties of good metals and several semiconductors. Other discrepancies show up for materials like elemental actinide solids. For instance, the equilibrium volume of -plutonium is experimentally 25% larger than the one given by the DFT-LDA approach, the greatest deviation known between experiment and theoretical value for this theory. The inadequacy of the DFT-LDA method for strongly correlated electron materials can be partly cured by including a direct treatment of quantum fluctuation effects by such quantum many-body approaches like the dynamical mean-field theory (DMFT). AGeorges96 (); GKotliar06 () Together with its success in describing key physical observables in many strongly correlated electron materials, however, the LDA+DMFT is computationally expensive and in practice limited to solid state systems with high crystalline symmetry, making it time consuming to describe the structural relaxation problems. The combination of LDA with the Gutzwiller variational method Gutzwiller1963 (); Gutzwiller1965 () has proved successful in providing an alternative but fast approach to the strongly correlated electron metals. Julien2006 () When combined with the MD, the computational efficiency of the GA method makes it ideal for the studies of such problems as material structure stability in strongly correlated electron metals.

The strategy seems to be straightforward in principle. However, practical MD simulations for low symmetry structures (such as defects, surfaces, clusters, and liquids) containing thousands of atoms (and accompanying electrons) present challenges. On the one hand, for an explicit treatment of strong electronic correlation effects, the ab initio method requires a definition of local correlated orbitals. On the other hand, interatomic forces derived from the correlated wavefunctions need to be calculated rapidly and repeatedly during the time evolution of the MD trajectory. The aim of this work is to present a generic framework of the GA-MD method, together with a path forward for improving the computational speed of the GA optimization procedure. We propose the construction of such parameterizations as presented in the tight-binding electronic structure method, and the use of the semi-empirical second moment approximation of electronic density of states for the calculation of local kinetic energy. The latter will significantly speed up the minimization procedure in the Gutzwiller variational method for the electronic structure, opening up the possibility of MD simulations to strongly correlated electronic materials.

The outline of the paper is as follows. In Sec. II, we give a detailed description of the density matrix formulation of the Gutzwiller approximation. It has the advantage of being applicable to crystals, as well as topologically and/or chemically disordered systems. In Sec. III, we derive an approximate but analytical solution to the optimization equations in the GA, where a high quality fitted solution on the whole range of physical interest is provided, and propose the second moment approach to the electronic density of states for the calculation of kinetic energy parameters. In Sec. IV, this efficient GA-MD method is demonstrated in a minimal Anderson model for heavy fermion systems based on tight-binding hopping integrals. A concluding summary is given in Sec. V.

## Ii Density matrix formulation of Gutzwiller method

### ii.1 Renormalization of hopping integrals

First, we review the Gutzwiller method briefly for which we closely follow Ref. Julien2006, but here we specifically include the topological disorder, i.e., during a typical MD process, no symmetry remains and all atoms are inequivalent. Among numerous theoretical approaches, the Gutzwiller method provides a transparent physical interpretation in terms of the atomic configurations of a given site. Originally, it was applied to the one-band Hubbard model Hamiltonian: Hubbard1963 ()

(1) |

with

(2) |

and

(3) |

The Hamiltonian contains a kinetic part with a hopping integral from site to , and an interaction part with a local Coulomb repulsion for electrons on the same site. () is the creation (annihilation) operator of an electron at site with up or down spin . measures the number (0 or 1) of electron at site with spin . The Hamiltonian, Eq. (1), contains the key ingredients for correlated up and down spin electrons on a lattice: the competition between delocalization of electrons by hopping and their localization by the interaction. It is one of the most widely used models to study the electronic correlations in solids.

In the absence of the interaction , the ground state is characterized by the Slater determinant comprising the Hartree-like wave functions (HWF) of the uncorrelated electrons, . When is switched on, the weight of the doubly occupied sites will be reduced because of the cost of an additional energy per site. Accordingly, the trial Gutzwiller wave function (GWF) is built from the HWF ,

(4) |

The role of is to reduce the weight of the configurations with doubly occupied sites, where measures the number of double occupations and is a variational parameter. This method corrects the mean field (Hartree) approach, for which up and down spin electrons are independent, and, overestimates configurations with doubly occupied sites. Using the Rayleigh-Ritz principle, this parameter is determined by minimization of the energy in the Gutzwiller state , giving an upper bound to the true unknown ground state energy of . To enable a practical calculation, it is necessary to use the Gutzwiller approximation, which assumes that all configurations in the HWF have the same weight.

Nozieres Nozieres1986 () proposed an alternative way which shows that the Gutzwiller approach is equivalent to the renormalization of the density matrix in the GWF. It can be formalized as

(5) |

The density matrices and are projectors on the GWF and HWF, respectively. is an operator which is diagonal in the configuration basis; where is a diagonal operator acting on site ,

(6) |

Here, is an atomic configuration of the site , with probability in the GWF and in the HWF respectively, whereas is a configuration of the remaining sites of the lattice. Note that this prescription does not change the phase of the wave function as the eigenvalues of the operators are real. The correlations are local, and the configuration probabilities for different sites are independent.

The expectation value of the Hamiltonian is given by,

(7) |

The mean value of the on-site operators is exactly calculated with the double occupancy probability, . Therefore, is the new variational parameters replacing . Using Eqs. (5)-(6), the two-site operator contribution of the kinetic energy can be written as,

(8) |

where and are the only two configurations of spin at sites and that give a non-zero matrix element for the operator in the brackets. The summation is performed over the configurations of opposite spin . The probabilities in the HWF depend only on the number of electrons, whereas the in the GWF also depends on .

After some elementary algebra, one can show that the Gutzwiller mean value can be factored into,

(9) |

where these renormalization factors are local and can be expressed as

(10) |

In Eq. (9), is shorthand for the expectation value of over the HWF , that is, , and similarly for the average over the Gutzwiller state . We have also used as shorthand for , that is, the average number of electrons on the considered “orbital-spin” in the HWF. In the simple case when the state is homogeneous and paramagnetic, all quantities becoming site- and spin-independent.

In Eq. (9), the term contributing to the kinetic energy, , is renormalized by a factor of , which is less than one in the correlated state, and equal to one in the HWF. This factor can be interpreted as a direct measure of the correlation effect. Indeed Vollhardt Vollhardt1984 () has shown that where is the effective mass and is the bare mass of the electron. Thus a close to corresponds to a weakly correlated electron system and a smaller value reflects enhancement of the correlation effect. Equation (7) leads to the variational energy per site, and for the homogeneous and paramagnetic state, is given by

(11) |

which can be minimized numerically with respect to the variational parameter . In the above expression, the factor 2 accounts for the two-fold spin degeneracy and is the kinetic energy per site and per spin identical at all sites and spins for an homogeneous HWF,

(12) |

In the case of half filling (), minimization is analytical, and provides the optimal choice for double occupancy

(13) |

and

(14) |

If the Coulomb repulsion exceeds a critical value , , leading to an infinite quasiparticle mass with a Mott-Hubbard Metal-Insulator transition. This is also known as “the Brinkmann-Rice transition”, Brinkmann1970 () as these authors first applied the Gutzwiller approximation to the Metal-Insulator transition.

Away from half-filling, one has to minimize the variational energy of Eq. (7) numerically. Moreover if the system is inhomogeneous, which is the case for a MD simulation, all quantities (, Â ) may vary locally from one site to the other. Consequently, the general variational energy, a function of double occupancy probabilities on all sites, is

(15) |

Minimization must then be performed numerically for each site, i.e., derivation with respect to , leading to the local equation:

(16) |

Here is the local partial “effective” kinetic energy, i.e., the contribution of orbital-spin “” to kinetic energy, but calculated with an “effective” hopping, renormalized by ,

(17) |

### ii.2 Inequivalent sites: renormalization of levels

When sites are inequivalent, or if orbitals belong to different symmetries as in a multiorbital basis, it is necessary to add to the Hamiltonian an on-site energy term

(18) |

The Hubbard Hamiltonian is written as

(19) |

In this case, the starting HWF directly obtained from the non-interacting part of the Hamiltonian, is not automatically the optimal choice, i.e., having the lowest energy. For example, if we look for the ground state of Eq. (19) in the Hartree-Fock (HF) self-consistent field formalism, it is necessary to vary the orbital occupations. Practically, it can be achieved by replacing Eq. (19), by an effective Hamiltonian of independent particles with renormalized on-site energies :

(20) |

The HWF we are looking for is an approximate ground state of the true many-body Hamiltonian (19) and is the exact ground state of effective Hamiltonian (20). The additive constant accounts for double counting energy reference, so that the ground state energies are the same for both Hamiltonians:

(21) |

The optimal choice of parameters can be obtained by minimizing the ground state energy of with respect to . With the help of Hellmann-Feyman theorem, one can obtain the derivative of the kinetic energy

(22) |

On the other hand, differentiation of Eq. (21) in association with Eq. (22) and the mean field approximation recovers the well-known formula for the on-site energies

(23) |

where the constant is simply .

In the Gutzwiller approach, the same argument about the variation of orbital occupation, i.e., flexibility on the HWF , is true. It is necessary to find a way to vary this Slater determinant, from which the GWF is generated, so that the Gutzwiller ground-state energy is a minimum. One needs to find an equivalent of Eq. (23) in the Gutzwiller context. The average value of Eq. (19) on a GWF is given by:

(24) | |||||

Following the previous HWF self-consistent field approach, one can find an effective Hamiltonian of independent particles having as an exact ground state. This state generates the GWF which is an approximate ground state of the true Hamiltonian Eq. (19). In analogy with Eq. (21),

(25) |

leads to the expression:

(26) |

with effective but fixed renormalized hopping integrals and effective on-site energies , still to be determined. The Hellmann-Feynman theorem applied to provides again an expression similar to Eq. (22), but with effective hopping integrals. Taking into account the dependence of the ’s through (Eq. (10)) and differentiating Eqs. (24) and (25) with respect to the parameters , one obtains the equivalent expression to Eq. (23) in the Gutzwiller context:

(27) |

Here is the partial kinetic energy of orbital-spin

(28) |

with the -projected density of states (DOS) for a system described . Equation (25) leads to

(29) |

Except for a few very special conditions in one-band Hubbard model, the renormalization of correlated-orbital levels is not only important in the optimization of the total energy but also in giving a correct description of single-particle quasiparticle properties. JXZhu2012 () To solve the full problem of finding an approximate ground state to Eq. (19), one is faced with a self-consistency loop: First get the occupations from a HWF, and a set of ‘bare’ levels; then obtain a set of configuration parameters, the probabilities of double occupation, by minimizing Eq. (24) with respect to these probabilities, followed by the on-site level renormalization according to Eq. (27). The loop repeated until a convergence is achieved.

## Iii Approximated solutions of minimization equations

Due to the complicated expression of Eq. (10), it is non-trivial to solve Eq. (16). Graphically, the solution corresponds to the intersection of the function with a horizontal line (see Fig. 1).

This situation may seem insurmountable for application of the Gutzwiller method to MD, as analytical expressions are necessary to be able to derive forces. Fortunately, the function can be fitted with reasonable accuracy (see Fig. 1) by a logarithm function, giving an analytical approximate solution of Eq. (16). This choice was suggested by the shape of the true derivative of , keeping in mind the following physical constraints: The uncorrelated case () has to give the solution , for a given occupancy , and the probability of double occupancy is restricted in the range (otherwise there would be negative arguments in the square root of ), providing a rescaling of the logarithm argument. Finally, we adapt the coefficient in front of the logarithm, in such a way that the fitted function has the same slope as the true one in the uncorrelated limit . The final result reads:

(30) |

The physical constraints above fix uniquely all three coefficients

(31) | |||||

(32) | |||||

(33) |

The approximate value of double occupancy, the solution of minimization equation, within this approximation, is , where . The small remaining difference between this approximate and the true value can be corrected by a second order expansion around the approximate value leading to an accurate analytical expression

(34) |

Here , , and stand for the true , and its first and second order derivatives, respectively, calculated at the approximated value .

The relative error on this second order corrected value with respect to the exact solution is less than 1% over the whole range of values. This second order corrected local double occupancy is the one now used in the calculation of renormalization factor , Eq. (10). To check the validity of this approximation for the derivative, we also plot in Fig. 1 the comparison between true and approximate . Again we see the good accuracy, the small discrepancy being in the range of very small double occupancy, i.e., corresponding to high value of Coulomb repulsion , which are far from the values for realistic materials.

The other input for Eq. (28) necessary to perform tractable MD simulations is the partial energy. A common approximation is the well-known approach of second moments JFriedel:1964 () with the assumption of a rectangular electronic density of states of bandwidth and height for each spin. The resulting second moment of the rectangular band is . For a given MD snapshot, the second moment for a given atomic site “” can be constructed as a sum over atoms “” neighboring “” as (see Appendix for more details). Using the simple tight-binding theory, WAHarrison:1980 () the hopping integrals scale as a power law of the interatomic distance . With the number of electrons on each atom (assuming charge neutrality for all sites), the partial kinetic energy needed as input in Eq. (28) and is given by , similar to the result by Ackland. Ackland2003 () To compute and to account for the effect of Coulomb correlations, we use the hopping integrals, renormalized by -factors, , rather than the bare ones, : After minimization in the Gutzwiller method, the true interacting Hamiltonian , is replaced by an effective Hamiltonian of non-interacting quasiparticles, with renormalized hopping integrals (and potentially renormalized on-sites too, but it is useless in our case, as we assumed charge neutrality, so there are no average charge transfer between sites). To conclude, we see that inclusion of electronic correlations in MD simulations within the Gutzwiller method, just requires one more intermediate step, compared to usual one, for having the renormalization factors that reduce a bit the hopping integrals. Once they have been computed for a set of actual positions of atoms, the rest of the process is similar to the usual way developed in other semi-empirical approaches. MFinnis:1984 (); MSDaw:1983 ()

## Iv Model and results

### iv.1 Model

To illustrate the possibility of the method, we consider a minimal two-orbital model that mimics, e.g., heavy fermions or actinides systems, with one non-correlated band, called for convenience “d”, whereas the other, called “f”, bears a strong local Coulomb repulsion . For simplicity, each of these two orbitals has a spin- degree of freedom per site. This model is described by the following Hamiltonian (close in spirit to the Anderson lattice model) with the usual notation:

(35) | |||||

Here the -orbitals are coupled among themselves and with -orbitals, whereas the -orbitals are only coupled to their neighboring -orbitals. The power laws WAHarrison:1984 in distance from atom located at position to atom at for hopping integrals are -5 and -6 for dd-coupling and df-coupling, respectively:

(36) |

and

(37) |

where and are constants. After a Gutzwiller mean-field treatment of Hamiltonian Eq. (35), we obtain

(38) | |||||

whose parameters are obtained from the minimization procedure analogous to deriving Eq. (26) from Eq. (25), for a given set of atomic positions. When the converged Gutzwiller ground state has been obtained we calculate the forces on each atom. These attractive forces have a quantum origin, due to the hybridization through hopping integrals. To mimic the short range repulsion that accounts for the Pauli principle when atoms get too close, we add a phenomenological repulsive potential

(39) |

with a constant.

### iv.2 Calculation of forces

For a given set of atomic positions, the overall total energy of the system is the sum of the electronic approximate Gutzwiller ground state energy plus the short range repulsion potential,

(40) |

The - component of the force acting on atom , is the derivative of with respect to position component of this atom (same relations hold for - and -components):

(41) |

From the Hellmann-Feynman theorem, the first term, due to hybridization, can be split into elementary contributions:

(42) |

where the contribution of orbitals of site and of site ( or are either - or -orbitals) is related to the derivative of the hopping integral ,

(43) |

The factor 4 arises from the two-fold spin degeneracy and the Hermiticity. For the interacting case (), or are less than one for -orbitals but equal to one for -orbitals. For the non-interacting case (), the above formula is obtained by setting all . It can be shown that forces from the hybridization origin are always attractive.

The computed forces Eq. (41) are then inserted into New’s equation of motion (EOM) for each atom. The positions are advanced in time by a time step by numerically integrating the EOMs with the Verlet algorithm. The resulting new atomic positions are then taken as input into Eq. (38), and new atomic forces Eq. (41) are computed. The MD trajectory consists of the string of many time steps iterating back and forth through this two-step process.

### iv.3 Results

As a demonstration, consider a two-dimensional system that contains 16 atoms. In the calculation, we take , , and . Hereafter all energies are measured in units of . The bare level is chosen to be . The initial condition is 16 atoms forming a regularly spaced 4 x 4 array with a unit nearest-neighbor distance. To benchmark our method and see the efficiency of second moment plus accurate approximate solution for double occupancy, we also performed the calculation based on exact diagonalization. Since we are interested in finding only the equilibrium structure of the system, the velocity on each atom is set to zero before advancing by the numerical solution of the EOM. The resulting MD “trajectory” eventually finds a local minimum on the energy landscape as the atomic positions are converged and the residual forces are driven to the noise limit. We started with the non-interacting case (). Figure 2 shows initial and final positions of of atoms after 5 millions of MD time steps i.e., sufficient to converge and where residual forces can be considered as noise. The equilibrium structures obtained from the second moment and the exact diagonalization methods are quite similar with the root mean square deviation of only about 0.058 unit distance. The CPU time in the case of second moment method is a factor of 20 times faster than the exact diagonalization. The second moment approach scales linearly with the number of atoms whereas ED scales as . Therefore, if we had studied 10 times more particles, i.e., 160, the increase speed factor would have been around 2000. The combination of second moment for approximate electronic structure quantities and a fast Gutzwiller solver really opens the realm of possibilities for the simulations of large systems, where electronic correlations play an important role.

We then started from the converged atomic positions as shown in Fig. 2 and repeated the simulation with a local Coulomb repulsion . The comparison between and for the second moment approach is displayed in Fig. 3. As expected and in accordance to our previous experience with Pu-, Julien2006 () we see that the Gutzwiller -factors have the effect of reducing the hybridization and the resultant attractive forces, which leads to a slightly more expanded equilibrium structure as shown in the zoom-insets of Fig. 3. The same trend was also observed in the results (not shown here) obtained by exact diagonalization with versus . We note that the structure expansion is quite small. The reason lies in the fact that in the Anderson-like model, the attractive force is contributed not only from the - hybridization hopping but also significantly from the direct - hopping. In addition, the efficiency of the hybridization reduction in the present model is proportional to . It is in contrast to the one-orbital Hubbard model, where the effective hopping integrals are proportional to .

## V Summary and conclusion

We have derived for the first time a real-space version of the Gutzwiller method embedded into the MD simulations for strongly correlated electron systems. From the positions of a configuration of atoms, a Hamiltonian can be constructed in terms of hopping integrals, on-site energies and Coulomb repulsion terms. It is precisely these interaction terms that require treatment beyond mean-field HF theories, which allows for Gutzwiller method. This method is a variational method in the Rayleigh-Ritz sense, for which one minimizes the energy of the system via a set of local Gutzwiller variational parameter dependent on each site, thereby providing an approximate ground state energy for the system. This minimization can be computationally demanding, especially when all the sites are inequivalent. This has motivated the development of an accurate analytical solution for finding the optimized double occupancy on each site.

MD simulation is a repeated two step process. First, employing a Hellmann-Feynman theorem within the Gutzwiller ground state, we can calculate the quantum origin of the forces acting on the atoms. This ground state energy defines an interatomic potential, which explicitly accounts for correlation effects. Second, the atomic forces derived from this interatomic potential are input into the classical equations of motion and the atomic positions are evolving forward in time for one time step using a numerical integrator (e.g., the Verlet algorithm). The step one is then repeated, where the new atomic positions define a new Hamiltonian whose new ground state can be found by the Gutzwiller method and so on.

A further approximation consists of avoiding exact diagonalization of the Hamiltonian from which, in principle, we can calculate the local densities of states (DOS) to obtain all necessary integrated quantities. Instead, we have proposed to use the second moment approximation, in which the true DOS is replaced by a rectangular approximation having the same second moment. Because the needed quantities to construct the variational ground state energy are basically integrated from DOS, they are less sensitive to the detailed structure of the DOS, thus validating the second moment approximation. The second moment of the energy is easily computed from a few Hamiltonian matrix elements, where we can set up the MD process without invoking exact diagonalization to solve iteratively the Gutzwiller minimization. Therefore, a very accurate approximate, but analytical solution, is available, making this process feasible. We concluded this first study with an application to a realistic case to show the potential of this approach, which to our knowledge, has never been developed and will open up new possibilities for simulations of correlated electron materials with molecular dynamics. We have applied the Gutzwiller Molecular Dynamics method to the one correlated orbital per site case as exemplified in the lattice Anderson model. The generalization of these ideas to multiple correlated orbitals case is not made directly because of higher degeneracy. That is, the number of variational parameters increases as the number of atomic configurations (namely as , with the degeneracy of the level) and therefore the number of local equations to be solved increases accordingly. We defer to future work on how to reduce the number of variational parameters, and as in the one correlated orbital per site model, to find an analytical approach to the problem.

###### Acknowledgements.

We thank S. Valone, W. A. Harrison, and J. M. Wills for useful discussions. J.-P.J. would like to thank the Los Alamos National Laboratory for the hospitality and financial support during his visits. This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at LANL under Contract No. DE-AC52-06NA25396, and was supported by the LANL ASC Program (J.D.K. & J.-X.Z.). This work was in part supported by the Center for Integrated Nanotechnologies, a U.S. DOE user facility. Some preliminary results have been reported on the 2014 CECAM Workshop on Gutzwiller Wave Functions and Related Methods.*

## Appendix A Second moment approximation for bond quantities

We follow closely the derivation by Pettifor and co-workers APSutton:1988 (); DGPettifor:1989 () for the bond order in a tight-binding model. To compute forces from the derivation of Hamiltonian, one needs average values like , with and being a short hand for site and spin-orbital states. Within these notations, and for the purpose of demonstration, we write the Hamiltonian in simple tight-binding form:

(44) |

where on-site energies can be identified as average value of the Hamiltonian on local state labelled by whereas hopping integrals couple states to . The bracket is the thermal average obtained for a general operator by

(45) |

where , , and are respectively the Hamiltonian, the chemical potential, the operator number of particles and the grand partition function. This average reduces to the ground state mean value at zero temperature. From exact diagonalization (as we do for the cluster example developed here), these quantities can easily be calculated from the weights of state (and similarly ) on eigenstate labelled by and of energy :

(46) |

where is the Fermi distribution.

For large systems, where the speed of calculation is a limiting factor, it might be desirable however to avoid this diagonalization, and to find an approximate but cheap way to get them. In a alternative route, they are obtained from the retarded Green function (with Zubarev notation Zubarev60 ()) , which is the Fourier transform of for the operators and . can be considered as an off-diagonal element of the Green function, and in the case of effective independent electrons (with only 1-body operators, as it is the case for effective Hamiltonians) it reduces to the usual resolvent:

(47) |

The diagonal element (i.e. same states and ) relates to the -projected density of states via:

(48) |

There are several useful relations and sum rules that fulfills the Green function (see Ref. Zubarev60, for demonstration):

(49) | |||||

and the so-called spectral theorem provides a direct way to compute the average values we are looking for:

(50) |

When there is no magnetic field, the following relation holds:

(51) |

The main idea in the present approximate calculation is the following: we express exactly the off-diagonal element of Green function as a linear combination of diagonal elements of Green functions of bonding and anti-bonding states,created by . Indeed, one can easily show

(52) |

from which the desired quantity is obtained thanks to relations (50), (48) and (51):

(53) |

Then the second moment approximation is applied to the DOS –the imaginary part of Green function– calculated for bonding and anti-bonding states, respectively. The second moment approximation is based on the constraints that all necessary quantities are integrated quantities. Consequently, they are not sensitive to the fine details of the DOS, which will be replaced by rectangular DOS having the same second moment that the true ones. The rectangular -projected (centered on site energy with width and height ) DOS has its second moment given by: whereas a direct path-counting (see Ref. Cyrot-Lackmann73, ) gives . Identification between those two relations fixes uniquely the bandwidth , from which the approximated -projected DOS can be computed. This procedure has been widely used in semi-empirical molecular dynamics as in Ref. Ackland2003, , for approximate local DOS. The extension we suggest here is to used it also for the bonding and anti-bonding DOS, and , to obtain the bond quantities.

It should be noted that in the particular case of the second moment of the projected density of states on A (or B), the second moment defined by the mean value provides unusual terms as follows:

(54) |

Such terms do not appear in the traditional way of calculating the moments starting from a given orbital at a specific site , where the second moment is simply given by the number of paths starting from , returning after two jumps back to the same site: this later simpler case is obtained by setting in the last equation. Finally, we get the second moment on either or combination:

(55a) | |||

(55b) |

together with the center of band (resp. ):

(56a) | |||||

(56b) |

from which we can obtain the related bandwidth (similar expression holds for ). Finally using expression (53) for rectangular DOS and , we have the desired quantity:

(57) |

where is the chemical potential at zero temperature. All this demonstration can be straightforwardly be extended to multiband case adding orbital index and spin to labels and . This procedure presents the great advantage of being very rapid compared to exact diagonalization and one can check that the Green function related to approximate DOS also fulfills sum rule as given in Eq. (49).

## References

- (1) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- (2) For a review, see G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
- (3) M.C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
- (4) M.C. Gutzwiller, Phys. Rev. 137, A1726 (1965).
- (5) J.-P. Julien and J. Bouchet, Prog. Theor. Chem. Phys., B 15 , 509-534 (2006)
- (6) J. Hubbard, Proc. Roy. Soc. London, A 276, 238 (1963).
- (7) P. Nozières, Magnétisme et localisation dans les liquides de Fermi, Cours du Collège de France, Paris (1986).
- (8) D. Vollhardt, Rev. Mod. Phys., 56, 99 (1984).
- (9) W.F. Brinkmann and T.M. Rice, Phys. Rev. B, 2, 1324 (1970).
- (10) J.-X. Zhu, J.-P. Julien, Y. Dubi, and A. V. Balatsky, Phys. Rev. Lett. 108, 186401 (2012).
- (11) J. Friedel, Trans. Metall. Soc. AIME 230, 616 (1964).
- (12) W. A. Harrison, Electronic Structure and the Properties of Solids (W. H. Freeman and Co., San Francisco, 1980).
- (13) G. J. Ackland and S. K. Reed, Phys. Rev. B 67, 174108 (2003).
- (14) M. Finnis and J. F. Sinclair, Phil. Mag. A 50, 45 (1984).
- (15) M. S. Daw and M. I. Baskes, Phys. Rev. Lett. 50, 1285 (1983).
- (16) J. Bünemann, W. Weber and F. Gebhard, Phys. Rev. B 57, 6896 (1998).
- (17) J. Bünemann, F. Gebhard, and W. Weber, J. Phys.: Condens. Matter, 9, 7343 (1997).
- (18) J. Bünemann, and W. Weber, Phys. Rev. B 55, 4011 (1997).
- (19) T. Okabe, J. Phys. Soc. Jpn. 66, 2129 (1997).
- (20) H. Hasegawa, J. Phys. Soc. Jpn. 66, 1391 (1997).
- (21) J.P. Lu, Inter. J. Mod. Phys. B10, 3717 (1996).
- (22) S.Y. Savrasov and G. Kotliar, Phys. Rev. Lett. 84, 3670 (2000).
- (23) F. Gebhard, Phys. Rev. B 44, 992 (1991).
- (24) D.N. Zubarev, Soviet Phys. Uspeskhi 3, 320(1960).
- (25) J.P. Gaspard and F. Cyrot-Lackmann, J. Phys. C: Solid State Phys. 6, 3077 (1973).
- (26) A. P. Sutton, M. W. Finnis, D. G. Pettifor, and Y. Ohta, J. Phys. C: Solid State Phys. 21, 35 (1988).
- (27) D. G. Pettifor, Phys. Rev. Lett. 63, 2480 (1989).