TRIQS/DFTTools: A TRIQS application for ab initio calculations of correlated materials

We present the TRIQS/DFTTools package, an application based on the TRIQS library that connects this toolbox to realistic materials calculations based on density functional theory (DFT). In particular, TRIQS/DFTTools together with TRIQS allows an efficient implementation of DFT plus dynamical mean-field theory (DMFT) calculations. It supplies tools and methods to construct Wannier functions and to perform the DMFT self-consistency cycle in this basis set. Post-processing tools, such as band-structure plotting or the calculation of transport properties are also implemented. The package comes with a fully charge self-consistent interface to the Wien2k band structure code, as well as a generic interface that allows to use TRIQS/DFTTools together with a large variety of DFT codes. It is distributed under the GNU General Public License (GPLv3).


Introduction and Motivation
When describing the physical and also chemical properties of crystalline materials, there is a standard method that is used with great success for a large variety of systems: density functional theory (DFT). This powerful theory states that the electron density of a system uniquely determines all ground-state properties. However, in order to find the electron density in practice one usually resorts to the Kohn-Sham scheme involving an approximate exchange-correlation potential, such as the local-density or the generalizedgradient approximations. The solution of the self-consistent Kohn-Sham equations yields the Kohn-Sham orbitals (also called Bloch bands) |ψ kν and the corresponding Kohn-Sham energies ε kν . The index ν labels the electronic band, and k is the wave vector within the first Brillouin zone (BZ). Below, we will always imply the Kohn-Sham scheme when talking about DFT.
The solutions of the Kohn-Sham equations form bands of allowed states in momentum space. These states are filled up to the Fermi level by the electrons according to the Pauli principle, and one can use this simple picture for example to explain the electronic band structure of materials such as elementary copper or aluminium. Following this principle one can try to classify all existing materials into metals and insulators. A system is a metal if there are an odd number of electrons per unit cell in the valence bands, since this leads to a partially filled band that cuts the Fermi energy and thus produces a Fermi surface. On the other hand, an even number of electrons per unit cell can lead to completely filled bands with a finite excitation gap to the conduction bands, i.e. insulating behaviour. 2 However, there are certain compounds where this classification into metals and insulators fails dramatically. This happens in particular in systems with partially filled dand f -shells. There, DFT predicts metallic behaviour due to the open-shell setting, but in experiments many of these materials actually show insulating behaviour. This cannot be explained by band theory and the Pauli principle alone, and a different mechanism has to be invoked. The insulating behaviour in such systems, among which are the Mott insulators [1], comes by solely due to correlations between the electrons. Moreover, many interesting properties, such as transport, photoemission, etc. are temperature-dependent, which cannot be treated with a ground-state theory such as DFT alone.
In order to describe correlated materials within ab initio schemes one has to go beyond a standard DFT treatment. One step in this direction is the so called DFT+U approximation [2] which can induce the opening of a gap through a static shift of the Kohn-Sham eigenvalues. However, this technique cannot properly describe the physics of correlated metals, whose low-energy properties are dominated by dynamical many-body effects. In the mid-nineties, a very successful method was introduced to tackle this problem, namely the combination of DFT with the dynamical mean-field theory (DMFT) [3,4,5]. The main idea is to add local quantum correlations in the framework of DMFT to the nonlocal quantum Hamiltonian obtained by DFT. This method has been since applied to a large variety of problems, including high-temperature superconductors, organic insulators, iron pnictides, oxide perovskites, and other strongly-correlated materials.
Since the first proposal of DFT+DMFT, there has been a tremendous amount of new developments of numerical tools, including (fully charge self-consistent) interfaces between DFT codes and the DMFT [6,7,8,9,10,11,12], and exact algorithms for solving the DMFT equations [13,14,15,16,17,18]. As a result, many more classes of materials can be studied within this framework on an ab initio basis, with much higher accuracy than just a decade ago. Importantly, these technical advances allow studies of realistic materials at experimentally interesting ambient and even lower temperatures (see, e.g., Ref. [19] for Sr 2 RuO 4 ), in contrast to the older techniques that allowed studies down to temperatures of 1000 K only (see, e.g., Ref. [20] for the same material as above). The drawback of these new methods, however, is that they require a substantial amount of coding effort to keep up with the state-of-the-art techniques in the field.
In this paper we present version 1.3 of the TRIQS/DFTTools application, which is based on the TRIQS library [21]. The purpose of this application is to provide a complete set of tools to perform fully charge self-consistent ab initio calculations for correlated materials, and it enables researchers in the field of correlated materials to do these calculations without dedicating time to code their own implementation. The TRIQS/DFTTools package provides a fully charge self-consistent interface to the Wien2k [22] band structure program, as well as an interface that is usable with many other DFT or tight-binding codes, however, without the possibility of full charge self-consistency. Physical properties that can be calculated range from spectral functions and (projected) densities of states (DOS), charges and magnetic moments, to more involved quantities like optical conductivity and thermopower. The package can handle spin-polarised cases, as well as spin-orbit coupled systems [23]. It is released under the Free Software GPLv3 license.
The paper is organised as follows. We start with a brief general discussion of the structure of TRIQS/DFTTools in Sec. 2. In Sec. 3 we show how the interface between DFT and DMFT is set up, in particular for the Wien2k DFT package. We will also introduce the projective Wannier functions that are used as the basis set for the DMFT calculations. We continue in Sec. 4 with a simple example of a DFT+DMFT calculation, without charge self-consistency. In Sec. 5 the concept of full charge self-consistency is introduced, followed by a discussion of post-processing tools in Sec. 6. Information on how to obtain and install the code is given in Sec. 7, and we finally conclude in Sec. 8.

General structure of TRIQS/DFTTools
The central part of TRIQS/DFTTools is a collection of Python modules that allow an efficient implementation of the DMFT self-consistency cycle. Since the input for this DMFT part has to be provided by DFT calculations, there needs to be another layer that connects the Python-based modules with the DFT output. Naturally, this interface depends on the DFT package at hand. Within TRIQS/DFTTools we provide an interface to the Wien2k band structure package, and a general interface that can be used in a more general setup. Note that only the Wien2k interface allows fully charge self-consistent calculations. Support for other codes is in progress.
The basic structure is schematically presented in Fig. 1. After performing a DFT calculation, the Kohn-Sham orbitals are used to construct localised Wannier orbitals, and all required information is converted into an hdf5 file. This file is used by the Python modules of TRIQS/DFTTools to perform the DMFT calculation. The fully charge self-consistent loop can be closed by taking the interacting density matrix and using it to recalculate the ground state density of the crystal in Wien2k. This leads to a new Kohn-Sham exchange-correlation potential, and thus to new orbitals. The full loop has to be repeated until convergence.
TRIQS/DFTTools also provides tools for analysing the data, such as spectral functions, densities of states, and transport properties. Depending on the method of solution of the Anderson impurity problem within the DMFT loop, the user may need to perform analytic continuation of the self-energy in order to use these tools 1 .

Projective Wannier functions
In order to perform the DMFT calculations one has to choose an appropriate localised basis set. In this package, we provide a method to calculate projective Wannier functions 1 Note that TRIQS/DFTTools itself does not provide any method to perform the analytic continuation. |w α,σ km based on the Wien2k Kohn-Sham orbitals, as introduced in Ref. [24]. We denote orbitals in the correlated subspace C where many-body correlations are applied by Latin letters m, the index of the atom in the unit cell by α, and the spin projection by σ.
The basic idea of the projective Wannier function technique [25] is to expand a set of local orbitals |χ α,σ m over a restricted number of Bloch bands where ν is the band index and the sum is carried out over Bloch bands within an energy window W. Note that the functions |χ α,σ km are not orthonormal due to this truncation. The orbitals |χ α,σ km can be orthonormalised, giving a set of Wannier functions: km ′ are the overlap matrix elements. The functions |w α,σ km are, in fact, the Bloch transform of the corresponding real-space Wannier functions |w α,σ Rm , through |w α,σ km = R e ikR |w α,σ Rm . In practice it is much more convenient to work with projection operators, which transform quantities from the Bloch band basis to the Wannier orbital basis. We define these operators asP The matrix elements of these operators are calculated straight-forwardly using the above Wannier functions in the latter case, without weight around the oxygen positions. The Wannier functions are constructed using the dmftproj program, their real-space representations are generated with the help of the wplot program from the wien2wannier [10] package. Plots are produced using XCrysDen [26].
Wannier functions: These projectors can now be used to calculate, e.g., the local projected non-interacting Matsubara Green's function from the DFT Green's function, where we have dropped the spin index for better readability. Note that throughout the paper we assume a proper normalisation of the momentum sum over k in the first Brillouin zone, i.e. k 1 = 1. From this non-interacting Green's function one can calculate the density of states. As an example, we show in Fig. 2 the DOS for the prototypical material SrVO 3 . Three different Wannier constructions are displayed. First, a projection where only the vanadium t 2g bands are taken into account is shown in black. The projection that comprises vanadium t 2g as well as the oxygen p bands is shown in red, and the green line is the DOS for a projection using DFT bands up to 8 eV. The difference in the latter two projections is minor, and consists primarily in the small transfer of weight to large energies around 7 eV.
In Fig. 3 we show the real-space representation of the Wannier charge density. The left plot shows a t 2g -like Wannier function for a projection of t 2g bands only, which corresponds to the black line in Fig. 2. The right plot is the corresponding Wannier 6 function when all bands are included in the projection. The effect of increasing the projection window is obvious: it results in a better localisation of the resulting Wannier functions.
All necessary steps to construct the orthonormalised set of projection operators from Wien2k DFT calculations is contained in the program dmftproj that is shipped with TRIQS/DFTTools. In Wien2k the Bloch functions are expanded in the LAPW basis set as where N P W is the total number of plane waves in the basis set, N lo is the number of local orbitals, and N LO the corresponding number of auxiliary orbitals for semi-core states. For more details on the Wien2k basis set we refer the reader to the literature [22,27,24]. As initial localised orbitals |χ α,σ m in our projective Wannier function construction, we choose the solutions of the Schrödinger equation within the muffin-tin sphere u α,σ l (E 1l )Y l m at the corresponding linearisation energy E 1l . Inserting this ansatz into Eq. 1 and using the orthogonality relations of the solutions of the Schrödinger equation allows to calculate the matrix elements of the projection operators. After some algebra one arrives at the expressions with contributions from the LAPW and/or APW+lo orbitals as well as the contribution from the semi-core orbitals, dmftproj takes all the necessary data from Wien2k, evaluates Eqs. 8 and 10 to obtain the matrix elements of the projection operators, Eqs. 4 and 5. Further information needed for the DMFT calculation are i) the rotation matrix from local to global coordinate systems, ii) the Kohn-Sham Hamiltonian within the projection window, and iii) the symmetry operation matrices for the BZ integration. All this information is extracted by dmftproj from the Wien2k output. 7

Usage
As a prerequisite a self-consistent Wien2k DFT calculation has to be performed. At this point, the further steps for the interface are • writing of necessary data from the Wien2k calculation into files, which is done using the Wien2k command x lapw2 -almd, • running dmftproj to calculate the projective Wannier functions, • conversion of the text output files into the hdf5 file format using the Python modules provided.
After these steps, all necessary information for a DMFT calculation is stored in a single hdf5 file. For the details of these steps, and all input/output parameters we refer the reader to the extensive online documentation of TRIQS/DFTTools and the example below.

General interface for non-Wien2k users
The TRIQS/DFTTools package also contains an interface layer that can be used in many situations, albeit without the possibility of full charge self-consistency. In fact, as an input it takes the correlated subspace Hamiltonian H mn (k), where m and n are orbital indices, and the k-mesh spans the first BZ. Unlike the Wien2k converter, the calculation is done exclusively in Wannier orbital space. The converter takes this H mn (k), together with information such as the quantum numbers of the correlated orbitals and the required electron density, and produces a hdf5 file that can be used for the DMFT calculations. We note that this interface is completely independent of how the H mn (k) is produced. It can use, for instance, output from Wannier90 [28,29], which is available for many DFT codes, from NMTO calculations [30], or from any other tight-binding model representing the non-interacting electronic structure of the system.

One-shot DFT+DMFT calculations
Listing 1 Sketch of the Python script for the DMFT self-consistency loop in one-shot DFT+DMFT calculations. Note that the solver imports, construction and execution must be appropriately added and modified. As mentioned in the introduction, TRIQS/DFTTools provides simple modules that allow the compact and efficient implementation of the DMFT self-consistency cycle. We will illustrate how to use the package by giving a one-to-one correspondence of the code statements to the necessary equations. In this section we give a simple Python script for so-called 'one-shot' DFT+DMFT calculations without full charge self-consistency. Note that the construction of the Wannier functions and the conversion procedures discussed in the previous section have already been performed. We assume that the hdf5 archive for the calculation, where everything is stored, is called YourDFTDMFTcalculation.h5. Furthermore, we assume that a suitable method to solve the quantum impurity problem is available, such as the TRIQS/CTHYB solver [18]. Note that TRIQS/DFTTools itself does not provide any impurity solver.
The DMFT calculation as given in Lst. 1 consists of the following steps. First, we initialise our calculation as done in lines 2 to 6 of Lst. 1. The SumkDFT module is loaded and initialised with the previously generated hdf5 archive YourDFTDMFTcalculation.h5. We initialise a solver module, which is used for the solution of the Anderson impurity model, for instance, the TRIQS/CTHYB package [18]. It is important that the user is familiar with the details of this solver initialisation, including set-up of the interaction matrix, definition of the interaction parameters U and J, etc. We request n loops = 15 self-consistency cycles to be done. This number has of course to be adapted to the problem at hand in order to reach convergence. After this initialisation step, we begin the self-consistency loop: 1) Impurity self-energy (line 11): The command SK.set_Sigma(Sigma_imp = [ S.Sigma_iw ]) takes an impurity self-energy Σ imp mm ′ (iω) in orbital space, and sets it as a member of the SK class. We assume that the object S.Sigma_iw is defined in the solver. At initialisation, the self-energy can be set to zero, or to its Hartree value to speed up convergence. 2) Chemical potential determination (line 13): The chemical potential is set according to the definition where N e is the required electronic density in the projection window. Please note that the operation [· · · ] −1 denotes matrix inversions throughout this paper. The self-energy in band space is defined by We use the projectors as introduced in Sec. 3, as well as a double-counting correction Σ DC [2]. 9 3) Local Green's function (line 15): The local Green's function is calculated as (13) where α is the index of the correlated atom in the unit cell. In Lst. 1 we have only one inequivalent correlated atomic shell, so we need to take only α = 0. 4) New hybridisation function (line 17): We now calculate the new Weiss field for the impurity problem using the Dyson equation 5) Solve the impurity problem (line 20): We call the solver method S.solve with the relevant parameters. Again, the user must be familiar with the usage of the solver of choice and its parameters, in order to produce meaningful and physically sound results. The solver produces the Green's function and the self-energy of the interacting impurity problem. Following the TRIQS/CTHYB convention, we assume for the example given here that the Green's function and the self-energy are stored in S.G_iw and S.Sigma_iw, respectively. 6) Double counting (lines 23-24): We calculate a new double counting (DC) correction from the impurity density matrix (line 23). As additional parameters one needs to provide a Hubbard interaction value U and Hund's exchange J, which have to be in accordance with the definition of the interaction matrix. TRIQS/DFTTools offers three different DC corrections: The flag use_dc_formula=0 corresponds to the fullylocalised limit (FLL) correction, use_dc_formula=1 is a variant of FLL that is used for Kanamori-type Hamiltonians [31], and use_dc_formula=2 is the around-mean-field (AMF) correction. This step completes the self-consistency loop.
Steps 1 to 6 are now iterated until convergence. After that, we can save the final solution into the hdf5 archive, as done in lines 27 to 32. We want to note here that it is sometimes useful to save all data produced by the user in a separate subgroup, e.g. dmft results, in order to have a clean structure in the hdf5 archive. For details how to handle the hdf5 archive we refer the reader to Ref. [21] and the online TRIQS documentation.
We want to stress here that this simple script of about 30 lines 2 of Python code is sufficient for a full DFT+DMFT calculation, which is possible only thanks to the modular structure of TRIQS/DFTTools and TRIQS.

Full charge self-consistency
The one-shot DFT+DMFT calculation can be extended to a fully charge self-consistent calculation with only marginal additional effort. The basic concept is general, but the implementation will depend of course on the DFT code that you wish to use for the calculations. We are presenting the interface to the Wien2k DFT package, which is the first interface implemented in TRIQS/DFTTools. We describe this interface below. Work is in progress to include more DFT packages, in particular VASP [32,33].
First of all, we want to stress that the fully charge self-consistent calculation is controlled by Wien2k scripts. That is why a sound knowledge of Wien2k and how it is used is absolutely necessary for the considerations in this section.
Listing 2 Code snippet for the calculation of the corrected charge density due to correlations charge density with correlation effects. Therefore, the Python scripts that control the DMFT iteration have to be modified: • We do not start each DMFT iteration from scratch, but instead from the previous iteration. This is most conveniently done by storing all necessary data in the hdf5 archive.
• One DMFT iteration (n loops = 1) per full self-consistency cycle is usually sufficient. In cases where convergence is problematic, more iterations may be needed.
• In order to calculate the correction to the charge density, one has to calculate the density matrix from the correlated Green's function at the end of the DMFT step, as shown in Lst. 2. The first two lines set the chemical potential more accurately. The important step is line 3, where the correlated density matrix is calculated in Bloch band space where Σ νν ′ (k, iω n ) is given by Eq. 12. This non-diagonal density matrix is used instead of the Kohn-Sham density matrix to calculate the charge density in x lapw2 -qdmft, see Fig. 4. Line 4 in Lst. 2 saves properties to the hdf5 archive such that they are available in the next iteration.
Compared to one-shot calculations, the fully charge self-consistent calculations are in general not much more demanding. As a rough guide, the additional computational effort due to an increased number of self-consistency loops is normally around 50%. Of course, this extra effort depends largely on the problem at hand and whether the DMFT corrections to the charge density are relevant or not.
This charge self-consistent implementation also allows the calculation of total energies from DFT+DMFT. The detailed formulas and their derivation can be found in Ref. [34] and references therein.
For more options of the TRIQS/DFTTools modules, as well as a complete example script and tutorial for a charge self-consistent calculation, including total energy calculation, we refer the reader to the online documentation of TRIQS/DFTTools.

Post-processing
The main output of a DMFT calculation, be it one-shot or fully charge self-consistent, is the interacting Green's function and the self-energy. However, in many cases one is interested in other quantities that are more closely connected to observables of experiments. For that reason TRIQS/DFTTools provides methods to further process the data and calculate some physical properties • (orbitally-projected) density of states • transport properties (resistivity, thermopower, optical conductivity) For all these methods one needs to use a real-frequency self-energy. In case the solver provides only quantities on the Matsubara axis, one has to analytically continue the data, by methods like Pade approximants [21], the maximum entropy method [35], or others.
DOS and band structure. These two quantities are very similar in terms of equations, since they are both imaginary parts of retarded Green's functions. For the band structure, we use where the Green's function is defined equivalently to Eq. 16, but using a real-frequency self-energy. The trace is taken in matrix space of Bloch indices ν, ν ′ . Summing this quantity over the first BZ gives the local spectral function.
As an example, we show in Fig. 5 a comparison of the DFT and the DFT+DMFT band structures (left) as well as the local spectral function and DOS (right) of SrVO 3 . The projective Wannier functions were calculated in a large energy window of [−8.5, 7.5] eV. The impurity problem was solved using the TRIQS/CTHYB solver with interaction parameters U = 6.0 eV and J = 0.65 eV, and an analytic continuation was performed using stochastic maximum entropy [36].
The thin black lines are the DFT results. One can clearly identify the mass renormalisation of the t 2g bands, which is a bit smaller than 2 in our calculation. The oxygen p, as well as the vanadium e g states are only marginally affected. However, due to the hybridisation between vanadium t 2g and oxygen p states, the bands between −7 and −2 eV acquire a finite width.
13 Transport properties. These are evaluated in the Kubo linear-response formalism [5] neglecting vertex corrections [38,39], and at the moment their computation is based on the Wien2k DFT code. In addition to the self-energy, one also needs the velocities v α (k) which are proportional to the matrix elements of the momentum operator in direction α = {x, y, z}, calculated with the Wien2k optics package [40]. The conductivity and the Seebeck coefficient for the combination of directions αα ′ are defined as with β the inverse temperature. The kinetic coefficients K αα ′ n are given by Here N sp is the spin factor and f (ω) is the Fermi function. The transport distribution Γ αα ′ (ω 1 , ω 2 ) is defined as where V is the unit cell volume. In multi-band systems the velocities v α (k) and the spectral function A(k, ω) are Hermitian matrices in the band indices ν and ν ′ . The frequency-dependent optical conductivity is given by A related implementation, also using full dipole matrix elements but in the framework of maximally localised Wannier functions, has recently been published as the woptic package [41].
We illustrate the calculation of transport properties again for the example of SrVO 3 . We apply the same setup as for the calculation of the band structure given above. For 14 the k-summation, we use 4495 k-points in the irreducible BZ. At T = 290 K we obtain a DC resistivity of 47 µΩcm, which compares reasonably well with the experimental value of 70 µΩcm [42]. The Seebeck coefficient, −8 µV/K in our calculation, agrees remarkably well with the experimental room temperature value of about −11 µV/K [43]. The same holds for the optical conductivity (Fig. 6) obtained from DFT+DMFT (black solid line) when compared to experimental results (blue dotted line). To show the effect of the correlations, we also present results without self-energy (red dashed line), which corresponds to evaluating the optical conductivity directly from DFT. Sometimes, scattering processes not included in plain DFT+DMFT calculations (e.g., phonons, impurities, or non-local fluctuations) are important. In this case, the calculated resistivity constitutes a lower bound to the expected experimental value. Note that for all transport calculations the convergence with the number of k-points is essential and has to be checked carefully. Further details on the use of these tools, including all options and parameters, are given in the online documentation.

Getting started
Detailed information on the installation procedure can be found on the TRIQS/DFTTools website and current issues and updates are available on GitHub.

Obtaining TRIQS/DFTTools
The TRIQS/DFTTools source code is available publicly and can be obtained by cloning the repository on the GitHub website at https://github.com/TRIQS/dft tools. As the TRIQS project, including its applications, is continuously improving, we strongly recommend users to obtain TRIQS and its applications, including TRIQS/DFTTools, from GitHub. Bugfixes to possible issues are also applied to the GitHub source.
For users that wish to have a more stable version of the code, without the latest bugfixes and functionalities, we suggest downloading the latest tagged version from the GitHub releases page https://github.com/TRIQS/dft tools/releases.

Installation
Once the TRIQS library has been installed properly, the installation of TRIQS/DFTTools is straightforward. As for TRIQS we use the cmake tool to configure, build and test the application. With the TRIQS library installed in /path/to/TRIQS/install/dir, the TRIQS/DFTTools application can be installed by $ git clone https://github.com/TRIQS/dft_tools.git src $ mkdir build_dft_tools && cd build_dft_tools $ cmake -DTRIQS_PATH=/path/to/TRIQS/install/dir ../src $ make $ make test $ make install This will first download the code to the src directory, then build, test, and install the application in the same location as the TRIQS library.

Citation policy
We kindly request that the present paper is cited in any published work using the TRIQS/DFTTools package. Furthermore, we request that the original works to this package, Refs. [24] for the first implementation and [34] for full charge self-consistency, are cited accordingly, too. In addition, since this package is based on the TRIQS library, a citation to [21] is also requested, along with the appropriate citation to any solver used, e.g. TRIQS/CTHYB [18].

Contributing
TRIQS/DFTTools, as an application of TRIQS, is an open source project. We appreciate feedback and contributions from the user community to this package. Issues and bugs can be reported on the GitHub website (https://github.com/TRIQS/dft tools/issues). Before starting major contributions, please coordinate with the team of TRIQS/DFTTools developers.

Summary
In summary we presented the TRIQS/DFTTools package, an application based on the TRIQS library. This package provides the necessary tools for setting up DFT+DMFT calculations. We have shown simple examples demonstrating how the code can be used to calculate Wannier functions, spectral functions, as well as transport properties of the prototypical material SrVO 3 .