Efficient implementation of the Gutzwiller variational method

Efficient implementation of the Gutzwiller variational method

Nicola Lanatà University of Gothenburg, SE-412 96 Gothenburg, Sweden    Hugo U. R. Strand University of Gothenburg, SE-412 96 Gothenburg, Sweden    Xi Dai Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China    Bo Hellsing University of Gothenburg, SE-412 96 Gothenburg, Sweden
July 9, 2019

We present a self-consistent numerical approach to solve the Gutzwiller variational problem for general multi-band models with arbitrary on-site interaction. The proposed method generalizes and improves the procedure derived by Deng et al., Phys. Rev. B. 79 075114 (2009), overcoming the restriction to density-density interaction without increasing the complexity of the computational algorithm. Our approach drastically reduces the problem of the high-dimensional Gutzwiller minimization by mapping it to a minimization only in the variational density matrix, in the spirit of the Levy and Lieb formulation of DFT. For fixed density the Gutzwiller renormalization matrix is determined as a fixpoint of a proper functional, whose evaluation only requires ground-state calculations of matrices defined in the Gutzwiller variational space. Furthermore, the proposed method is able to account for the symmetries of the variational function in a controlled way, reducing the number of variational parameters. After a detailed description of the method we present calculations for multi-band Hubbard models with full (rotationally invariant) Hund’s rule on-site interaction. Our analysis shows that the numerical algorithm is very efficient, stable and easy to implement. For these reasons this method is particularly suitable for first principle studies – e.g., in combination with DFT – of many complex real materials, where the full intra-atomic interaction is important to obtain correct results.

78.20.Bh, 71.10.Fd, 71.10.-w

I Introduction

In the 60th’s Martin Gutzwiller published a series of papers Gutzwiller (1963, 1964, 1965) where he introduced a variational method for studying ferromagnetism in transition metals. His brilliant idea was to variationally determine a projected wavefunction represented as


where the local operators improve the non-interacting wavefunction in accordance with the on-site interaction by modifying the weight of local electronic configurations.

In spite of its simplicity, the average values of any operator on can only be computed numerically for realistic lattice models, e.g., using variational Monte Carlo. Yokoyama and Shiba (1987); Capello et al. (2005) For this reason, Gutzwiller introduced an approximate scheme, known as the Gutzwiller approximation, to compute these average values analytically. Successively, the development of dynamical mean field theory (DMFT) Georges et al. (1996) has brought additional insights into the physical meaning of the Gutzwiller approximation. In fact, Metzner and Vollhardt showed that this approximation is exact in the limit of infinite coordination lattices, Metzner and Vollhardt (1987, 1988); Gebhard (1990) where the single-particle self-energy becomes purely local in space. Müller-Hartmann (1989)

Since its introduction the Gutzwiller wavefunction and approximation have proven to be very important tools to study strongly correlated systems. The understanding of many basic concepts, such as the Brinkman-Rice scenario for the Mott transition, Brinkman and Rice (1970) came originally from calculations based on the Gutzwiller method. From the computational point of view the list of interesting results that have been obtained by means of the Gutzwiller approximation Gutzwiller (1963, 1964, 1965) and its respective generalizations Dzierzawa et al. (1997); Gebhard (1990); Lanatà (2010); Schirò and Fabrizio (2010); Lanatà and Strand (2011) is impressively long, hence impossible to cite in an exhaustive way. Furthermore, the Gutzwiller approximation can be naturally combined with density functional theory (DFT), Hohenberg and Kohn (1964); Kohn and Sham (1965) applying, e.g., the local density approximation (LDA) Gunnarsson and Lundqvist (1976) for the exchange and correlation. LDA+Gutzwiller (LDA+G) Ho et al. (2008); Deng et al. (2009) has proven to be a powerful scheme for the study of real strongly-correlated metallic materials; Deng et al. (2009) giving a more accurate description than LDA+U, Anisimov et al. (1997) comparable to LDA+DMFT Kotliar et al. (2006) for ground state properties. Finally, the range of application of the Gutzwiller method has recently been extended to out-of-equilibrium calculations, such as electron transport across quantum dot systems Lanatà (2010) and quench dynamics in correlated electron systems. Schirò and Fabrizio (2010)

The Gutzwiller variational method requires a number of preliminary technical steps in order to make it really flexible and able to cope with systems of interest, as the number of variational parameters scales exponentially with the number of correlated orbitals involved in the calculation. This scaling is already problematic for transition metal systems with correlated orbitals if the required minimization is performed in a naive way. Based on the formalism introduced by Bünemann, Bünemann et al. (1998) Deng et al. Deng et al. (2009) recently derived a self-consistent numerical method that allows to efficiently perform calculations even for -orbital systems. The only limitation of this approach is the restriction to density-density type of interactions, which is actually due to the employed formalism and does not stem from the numerical approach in itself.

However, in order to properly describe the physics of several strongly correlated materials, the full intra-atomic interaction is needed; not only its density-density component. The rotational invariant on-site Coulomb and exchange interaction is generally modeled in terms of the so called Kanamori parameters Kanamori (1963) commonly referred to with the symbols and . The strength of the two-electron spin-exchange interaction is determined by the parameter . In transition metal oxides with partly filled -shells the off-diagonal interactions – exchange coupling, spin-flip and pair hopping – are crucial. For example, the metallic property of SrVO can not be reproduced from theory without accounting for spin-exchange interaction. Werner et al. (2009) Several theoretical model studies points in the same direction. To mention a few, the transition from the paramagnetic to the ferromagnetic phase for multiband systems requires finite Bünemann et al. (1998) the spin-freezing transition predicted for multiband systems, which is expected to influence the Mott transition, Costi and Liebsch (2007) takes place when and is absent for Werner et al. (2008) Previous Gutzwiller model studies Bünemann et al. (1998) and the present work show that the critical required for the Mott insulator transition is substantially reduced when increasing the ratio from zero.

Motivated by the above examples we believe that an implementation of the Gutzwiller variational method valid for general on-site interactions and, at the same time, numerically efficient would constitute an important progress. Formal advancements pointing in this direction have lately been conceived by several authors, Bünemann and Weber (1997); Gulácsi et al. (1993); Bünemann et al. (1998); Attaccalite and Fabrizio (2003); Wang et al. (2006); Fabrizio (2007); Lanatà et al. (2008, 2009) although these progresses have been hampered by the lack of efficient numerical algorithms applicable to the most general case. In particular, Fabrizio and collaborators Attaccalite and Fabrizio (2003); Ferrero et al. (2005); Fabrizio (2007); Lanatà et al. (2008, 2009); Lanatà (2009) derived a mathematical formulation of the problem whose complexity is unaffected by the form of the on-site interactions and furthermore allows to easily incorporate symmetries into the variational function from the onset.

The main goal of this work consists of merging together the general formalism developed by Fabrizio and collaborators Attaccalite and Fabrizio (2003); Ferrero et al. (2005); Fabrizio (2007); Lanatà et al. (2008, 2009) mentioned above and the numerical procedure derived by Deng et al.Deng et al. (2009) overcoming the restriction to density-density interaction without increasing the complexity of the computational algorithm. Furthermore, the fully rotational invariant on-site interaction enables us to construct the variational wavefunction using the orbital rotational symmetry, which is instead broken if the off-diagonal terms are neglected. It will be shown that this extra symmetry can be used to reduce the number of variational parameters.

The outline of this paper is as follows. In Sec. II the Gutzwiller problem for a general tight binding Hamiltonian is introduced. In Sec. III the employed formulation of the Gutzwiller method Attaccalite and Fabrizio (2003); Fabrizio (2007); Lanatà et al. (2008, 2009); Lanatà (2009) is summarized. In particular, in Sec. III.3 it is shown that this formulation provides a natural extension of the formalism of Ref. Deng et al., 2009, having the same mathematical structure in the special case of density-density on-site interactions. In Sec. IV we discuss the implementation of symmetries of the wavefunctions. In Sec. V the numerical procedure to minimize the Gutzwiller energy is described in detail. In Sec. VI we briefly discuss how the proposed method can be adapted to a of LDA+G type of calculation. In Sec. VII we prove the reliability of the method presenting a comparison with other methods for the cases of two and five orbitals. Furthermore we discuss several technical details of the numerical procedure, such as convergence properties and computational speed. Finally, Sec. VIII is devoted to the conclusions.

Ii The Gutzwiller method

Let us consider the general tight binding Hamiltonian

where creates an electron in state (where labels both the spin and the orbital at site ) and are many-body Fock states expressed in the -basis. These states are defined by the occupation numbers where runs over integer numbers from 1 to , being the number of on-site single particle states,


Thus the number of Fock states is . The Hermitian matrix represents the local terms, interaction and crystal fields, in the -basis, i.e., the same basis in which was defined in Eq. (II). This basis will henceforth be denoted as the original basis.

The structure of the Gutzwiller variational function is given by Eq. (1), where is an uncorrelated variational wavefunction, that satisfies Wick’s theorem, and is a general operator acting on the local configurations at site


where the matrix , assumed to be real in this work, contains all the variational parameters needed to define the operator .

In general, average values of operators with respect to must be computed numerically unless the lattice has infinite coordination number, in which case they can be evaluated analytically if the following equations – commonly named Gutzwiller constraints – are satisfied:


where is the local single-particle density-matrix operator with elements .

The variational problem to solve amounts to variationally determine both and by minimizing the average value of the Hamiltonian [Eq.(II)]


fulfilling Eqs. (5-6), where we have introduced


For a general tight binding model [Eq. (II)] this problem is complicated for two reasons: (i) and are not independent variables because of the Gutzwiller constraints [Eqs. (5-6)], (ii) the number of variational parameters scales exponentially with the number of orbitals.

Iii Reformulation of the Gutzwiller problem

In this section we briefly summarize the reformulation of the Gutzwiller problem derived in Refs. Lanatà et al., 2008Lanatà et al., 2009, and we show its formal analogy with the formulation of Bünemann and Weber Bünemann et al. (1998) in the special case of pure density-density local interaction.

iii.1 The mixed-basis representation

Let us introduce the so-called natural-basis Fabrizio (2007) operators , i.e., the operators such that


where are the eigenvalues of the local density matrix


Notice that the natural-basis operators are always well defined as is Hermitian, implying that there always exists a unitary transformation such that


Instead of expressing the Gutzwiller projector in terms of the original basis as in Eq. (4) we adopt the following mixed original-natural Lanatà et al. (2008) basis form


where, by assumption, are Fock states in the original -basis , while are Fock states in the natural basis, namely in terms of the -operators. In other words a generic state is identified by the occupation numbers – with – and has the explicit expression


For later convenience we adopt the convention that the order of the and the states is the same. For instance, if the second -vector in Eq. (12) is then the second -vector is .

iii.2 The -matrix

Let us introduce the uncorrelated occupation-probability matrix  Fabrizio (2007) with elements




We remind that are the elements of the diagonal density matrix of Eq. (9), and denote the occupation numbers of the natural states . We also introduce the matrix representation of the operators and


Notice that, if we respect the convention that the order of the and the states is the same, we have that


We now define the matrices and with elements  [Eq. (12)] and  [Eq. (II)]. Notice that is defined in the original-natural and is defined in the original-original basis. With the above definitions, the expectation value of any local observable can be calculated as




and the Gutzwiller constraints [Eqs. (5-6)] can be written as


The formalism is further simplified by defining the matrix


that was introduced in Ref. Lanatà et al., 2008. The expectation value of any local observable is given by


The Gutzwiller constraints take the form


and the variational energy [Eq. (7)] is, in the Gutzwiller approximation, given by Lanatà et al. (2008)




In conclusion, within the formalism summarized in this section, the variational energy is a functional of and , to be minimized fulfilling the Gutzwiller constraints [Eqs. (25-26)].

iii.3 Diagonal projector as a particular case

Let us assume that the coefficients of the matrix that define the projector in Eq. (4) are diagonal


and that the original basis coincides with the natural basis


From Eqs. (30-31) we have that


and, consequently, the following equation hold for all local operators


From Eq. (33) follows that the equations that characterize the Gutzwiller problem [Eqs. (25-29)] can be evaluated in terms of the variational parameters defined by Bünemann in Ref. Bünemann et al., 1998


A crucial observation in this work is that in the general case considered here, in which Eqs. (30-31) are not assumed, Eqs. (25-29) are expressed in terms of quadratic forms of the matrix elements of instead of , but in a formally identical way.

We conclude this section observing that the physical density matrix of the system


is not equal to the so called variational density matrix


when Eqs. (30-33) do not hold. In the general case the distinction between variational density matrix and physical density matrix needs to be taken into account.

Iv Symmetries of the variational function and matrix

In this section we discuss in detail how to build symmetries in the Gutzwiller variational function. The site label is dropped for simplicity. The procedure discussed here extends the method discussed in Ref. Lanatà et al., 2009 to general point symmetry groups. The problem amounts to define the form of the -matrix such that is invariant under the action of a matrix representation of a group in the many-body -local space. The transformation law


defines a representation Wigner (1959) of in the local Hilbert space generated by the Fock configurations , see Eq. (3).


In the original-original basis [Eq. (4)] the invariance condition


is equivalent to


Let us assume that the most general transformation that relates the original and the natural basis


commutes with , i.e., that


It can be shown, see appendix A.1, that Eq. (42) is verified for every group whose elements do not mix configurations that belong to different eigenspaces of the number operator . This assumption is obviously verified by every geometry group, but excludes, for instance, the particle-hole transformation. A first consequence of Eq. (42) is that the matrix has the same form in the original-natural and in the original-original representation. In fact


and from Eqs. (40) and (42) we have that


Notice that from the assumed invariance of respect to we have that, ,


Using Eq. (41) we can easily express in terms of as follows




so that, combining Eq. (45) and Eq. (46), we obtain that


Eq. (48) is equivalent, because of Eq. (42), to the following invariance relation for :


From the above considerations and Eq. (23) we can conclude that satisfies the same invariance relation of the coefficients of the Gutzwiller projector in the original-original basis, i.e., that


In other words, we have proven that has the same form of expressed in the original-original Fock representation [Eq. (4)].

The set of all the matrices that satisfy Eq. (50) is a linear space . Consequently, there exists a basis of matrices such that


In this work we assume that and are real. This means that is a linear space over the field of real numbers. Notice that this does not restrict the variational freedom as long as the local interaction is real. This excludes, for example, the spin-orbit coupling.

iv.1 Calculation of

In order to calculate it is convenient to apply a similarity transformation to


with the property to decompose in irreducible representations. Wigner (1959) More precisely, we need to calculate a unitary matrix such that is of the form


where (i) the representations are irreducible for all , (ii) if two representations are equivalent then they are also equal. In appendix A we derive a possible procedure to calculate explicitly the similarity transformation utilized in this section for a general geometry group .

Let us consider the linear space (over the real field) of the complex matrices that satisfy the following equation


It can be easily proven by means of the Schur lemma Wigner (1959) that the most general is of the form


where the blocks correspond to inequivalent representations of , and each block is of the form


being identity matrices of size , being the dimension of each one of the irreducible equivalent representations of repeated in the -th block, and being independent complex numbers. Eqs. (56) and (57) allow to define straightforwardly a basis of .

Let us define now the linear space of real matrices generated by the set of matrices obtained as


The linear space that we need, see Eq. (50), is obtained as


where is the linear space of all real matrices. It is convenient, finally, to orthonormalize the basis set of in order to have


iv.2 Independent variational parameters

All the relevant quantities that define the variational energy, i.e., (i) the Gutzwiller constraints [Eqs. (25-26)], (ii) the matrices [Eq. (29)] and (iii) the local-interaction energy [Eq. (27)] can be expressed in terms of quadratic forms in the coefficients defined in Eq. (52) as follows:


Notice that the tensors , and are fully determined by the symmetry of the wavefunction [Eq. (1)] and the number of orbitals. For this reason it is generally convenient to precalculate them before starting the numerical minimization of the variational energy. This point will be further discussed in Sec. VII.4.

iv.3 Simplified variational ansatz

The complexity of the numerical problem is considerably reduced if the mixing of different atomic configurations is neglected in the Gutzwiller projector. This amounts to assume that the matrix defined in Eq. (23) has the form


where are the orthogonal projectors onto the eigenspaces of the local atomic interaction . Although the variational parameters neglected in Eqs. (64-65) can play a crucial role in some case, Lanatà et al. (2008) this simplified ansatz merits to be mentioned for at least two reasons. (i) It still allows to solve exactly the problem in the atomic limit. (ii) The number of independent variational parameters is generally extremely lower in this approximation, allowing to perform calculations not feasible otherwise. Furthermore, once a variational result


is obtained assuming Eqs. (64-65), it can be used as a good starting point for the self-consistent search of the energy minimum with the more general variational space discussed before, see Eq. (52),


with the result to speed up the calculation. Notice, in fact, that , implying that


i.e., that .

V Numerical optimization of the variational energy

In this section we discuss in detail the self-consistent numerical strategy to minimize the energy [Eq. (27)] fulfilling the Gutzwiller constraints [Eqs. (25-26)].

Notice that the formulation of the Gutzwiller problem through Eqs. (25-27) is formally analog to the constrained formulation of DFT derived by Levy Levy (1979, 1982) and Lieb. Lieb (1983) In fact, the variational energy can be expressed as a functional of the variational density matrix , see Eq. (36),


where denotes the minimum over the set of variational parameters and , satisfying the Gutzwiller constraints [Eqs. (25-26)] at fixed . In this work Gutzwiller the problem is solved by calculating the density functional , and minimizing it respect to .

For clarity reasons we have structured the rest of this section as an exposition of the numerical procedure, omitting the mathematical proofs. The mathematical details can be found in the appendices.

v.1 Preliminary calculation

Let us consider the Gutzwiller renormalized non-local tight binding operator




For later convenience, we consider a general one-body Hamiltonian


where is a given local operator.

Let be the ground state of . It can be easily verified that




Eq (73) will be used in the following subsections, where two important inner parts of our numerical scheme are described in detail.

v.2 Slater determinant optimization step

For later convenience, in this section we solve the problem to calculate the state that realizes the minimum of the Gutzwiller variational energy at fixed


where is given by Eqs. (70-71). Note that the functional depends on only indirectly through , that is given by


see Eqs. (29) and (61).

It is convenient to account for the Gutzwiller constraints employing the Lagrange multipliers method. We introduce




Notice that has the same form of Eq. (72). Finally, we calculate the ground state of for such that the Gutzwiller constraints [Eqs. (25-26)] are satisfied. Once the Lagrange multipliers are known we compute Eq. (73).

In summary, the calculations described in this section associate the input variables and to the output matrix


see Fig. 1.

v.3 -matrix optimization step

In this section we derive the numerical procedure to minimize the Gutzwiller energy functional


where is given by Eqs. (70-71) and depends on through Eq. (77), keeping the Slater determinant fixed and respecting the Gutzwiller constraints


In order to solve this problem we adopt the following linearization procedure, that is based on appendix B. We consider the Hermitian matrix


where is defined by Eq. (80), and


Then, we calculate the ground state of for such that the Gutzwiller constraints [Eqs. (83-83)] are satisfied. The obtained vector is used to define a new through Eq. (77).

Notice that the matrix entirely encodes the dependency of the problem on in Eq. (84). In summary, the above calculations associate to the input variables and the output renormalizaton matrix , see Fig. 1.

v.4 Fixed point formulation

Figure 1: (Color online) Flow chart representing the numerical calculation of the functional .

A very important observation in our implementation is that the composition of the two optimization steps derived in Secs. V.2 and V.3 can be described as a functional that associates a given renormalization matrix to a new renormalization matrix


see Fig. 1. This operation lead to a reduction of the variational energy unless, by definition, solves the equation


In this case defines a stationary point of the energy functional. This observation amounts to formulate the minimization of the variational energy at fixed as a fixpoint problem, that can be solved in several ways. Heath (2002) A first possibility is to use as an input to obtain and iterate the procedure up to convergence. This procedure is commonly referred to as forward recursion method. Georges et al. (1996); Borghi et al. (2010) However, as it will be shown in Sec. VII.4, the application of the Newton method is generally much more efficient.

We underline that the size of is equal to the number of orbitals, and that the number of independent parameters that define it is often reduced by symmetry. For this reason the solution of Eq. (88) generally requires a few Newton steps to converge, see Sec. VII.4. The problem of the exponential scaling of the local many-body space affects the numerical algorithm exclusively through the solution for the ground state of , see Eq. (84). The size of the matrix is, in fact, equal to the dimension of the vector . Nevertheless, the calculation of the ground state of is not numerically problematic for two reasons. (i) The dimension of is reduced by symmetries, as will be shown in Sec. VII.4. (ii) The calculation of the ground state of does not require a full diagonalization. Less computationally demanding algorithms, such as the power method or the Lanczos method, can be employed. See Sec. VII.4 for further details.

We remark that the self-consistent numerical algorithm derived in this paper requires, as a starting point, only a initial “guess” for the variational density matrix and the matrix . It is not necessary to construct a good initial guess for the whole Gutzwiller wavefunction, i.e., for the matrix , while this would be necessary in order to perform a direct constrained minimization of the energy functional [Eq. (27)]. This implies that the stability of the algorithm is not affected by the exponential scaling of the number of parameters involved in the calculation. This point will be further discussed in Sec. VII.4.

Vi Application to LDA+Gutzwiller

In this section we briefly discuss how to combine the Gutzwiller scheme with a first principle calculation of the uncorrelated electron structure as an input, applying the DFT scheme with LDA. This combined scheme is named LDA+G. We also outline how the Gutzwiller solver has to be modified in order to account for the double counting. Anisimov et al. (1997) The double counting appears as a mean-field contribution in the exchange-correlation taken into account in the LDA calculation.

As a starting point we consider the Kohn-Sham reference system obtained within a converged LDA calculation


In order to be able to apply LDA+G the tight binding Hamiltonian [Eq. (89)] must first be expressed in terms of a proper localized basis set Deng et al. (2009) , where labels both spin and orbital . This can be done by means of the overlap matrix




In order to express in the same form of Eq. (II) we separate in a non-local part and in a local part (the crystal fields) as follows




and is the number of sites . The on-site electron interaction can be modeled by the Slater-Kanamori rotationally invariant atomic interaction Kanamori (1963)


We underline that the form [Eq. (101)] for is obtained by implicitly assuming that the single-particle basis is given by real orbitals, i.e., the cubic (crystal) harmonics. It can be proven that the condition ensures the rotational invariance in the orbital space. The additional condition