# Implementation of Dynamical Mean Field Theory

### Table of Contents

- Preliminaries
- Introduction
- Building the DMFT Hybridization Function - Routine lmfdmft
- Extracting the Impurity Self-Energy - CTQMC Software
- References

### Preliminaries

It is assumed that the user has a basic background on electronic structure methods. In particular, Density Functional Theory, GW methods and DMFT, as well as the LMTO basis set. For an overview of the theory, and in particular a specific account of the equations used in the implementation, we will refer to Paolo Pisanti’s PhD thesis [2], whose browsing is advised when reading this document.

The user should also have a basic understanding of the lmf code and familiarity with its usage. For this purpose we refer to lmf tutorial. In particular, this manual has to be considered complementary to a specific tutorial written by Lorenzo Sponza which can be found here and whose reading is encouraged.

### Introduction

A schematic flowchart clarifying the structure of the software and the main routines involved in the QSGW+DMFT loop is presented in the figure below. For a formal introduction on DMFT, its assumptions, main equations and key quantities into play we refer to chapter 3, section 3.1 of [2]. The details of our specific implementation of the QSGW+DMFT scheme are instead presented in chapter 4 of [2]. In particular in the next two chapters we will refer to the equations of section 4.2 and 4.3 to illustrate the several steps of the implementation and the different subroutines.

In chapter 1 we will introduce the structure, main input and output files of the **lmfdmft** routine (blue box in the image below). The object of this routine is the calculation, from first principles and using as inputs the results of a converged QSGW (or DFT) loop, of the so-called hybridization function , defined in DMFT to capture the coupling of the impurity with the surrounding bath in which is embedded. The key steps for such a purpose will be to define and calculate a projection operator to the correlated local subset, employ this operator to define local quantities such as local Green’s function and impurity levels, correct the eigenvalues with the embedded impurity self-energy injection, and tune the chemical potential accordingly.

In chapter 2 we will present the Continuous Time Quantum-Monte Carlo (CTQMC) impurity solver (red box in the image below), whose main operation is the calculation of the impurity self-energy . The software has been implemented by Kristjan Haule at Rutgers University and it was originally part of his Wien2k-DMFT code [1]. It has been isolated in such a way to work in connection with the LMTO suite in a QSGW+DMFT implementation. The CTQMC solver takes the hybridization function and the impurity level originating from the **lmfdmft** loop as main inputs, together with the effective interactions parameters . These parameters are tuned with respect to the material under study and they are typically specific to the DMFT implementation and the solver adopted.

The flowchart of the main routines of the present QSGW+DMFT implementation. The new interface (lmfdmft routine) has been singled out in a red box. This interface connects the main LMTO suite package to the CTQMC solver of DMFT. Image created by Lorenzo Sponza. The gold and green boxes appearing in the image represent the other programs necessary to complete the interface between the two packages. They will also be object of 2nd chapter.

### Building the DMFT Hybridization Function - Routine lmfdmft

This chapter is focused on the first block of the DMFT cycle, the calculation (and update) of the *hybridization function* from the inputs of a converged QSGW (or DFT) calculation. This operation is performed by the routine, part of the Full-Potential package.

A tutorial written by L. Sponza to run a full QSGW+DMFT calculation and in particular the executable can be found at this link.

#### Input Files

In this section an overview of the input files required and their syntax will be presented. The starting point of the DMFT loop are the results of a converged QSGW (or either a DFT) calculation, executed with the package.

In order to maintain consistency with the tutorial, we will refer to *lmfinput/* as the routine containing the original input files and *it#_lmfrun/* as the routine containing the #-th iteration of **lmfdmft**.

```
$ ls -l lmfinput
ctrl.case
site.case
rst.case
sigm.case
indmfl.case
sig.inp
```

An experienced user will notice that the first 4 files are regular files belonging to the full-potential package, those are the results of the converged QSGW loop. For a specific discussion of these files, we refer to this tutorial for **lmf**. The last 2 files (*indmf1.case* and *sig.inp*) are instead specific to the **lmfdmft** program and are only used in DMFT cycles.

##### Input Files: A Brief Overview

The is the main input file of, for a full documentation read here. In order to process this input file in a QSGW+DMFT calculation, one line has to be added to the input file resulting from the QSGW loop, this can be done by the following command

```
echo 'DMFT PROJ=2 NLOHI=11,53 BETA=50 NOMEGA=1999 KNORM=0' >> ctrl.case}.
```

The token *DMFT_PROJ* refers to the kind of projection operation, n.2 corresponds to equation (4.6) of [2]. No other option is advised at state of the art. The token *DMFT_NHOLI* defines the energy window in which the DMFT projection operator is included. In particular the first number refers to the first band to be included, lower bound of the window and the second number to the last band, upper bound. This energy window is supposed to include the contributions from the correlated subset of orbitals chosen (say, the -orbitals of Cu), as well as the contributions from orbitals which these bands are hybridized with (e.g. typically in copper oxides, the O -levels). We refer to section 4.1 of [2] for a guideline on the origin of this energy window and how to use it. The choice of the correct energy window should be inferred by looking at the QSGW bandstructure of the solid and in particular the point. The token *DMFT_BETA* refers to the inverse temperature, in units of eV, *DMFT_NOMEGA* to the total number of points of the Matsubara frequency mesh on which the dynamical impurity quantities are defined. *DMFT_KNORM* is a token related to the kind of renormalization to apply to the projection operator, the value of 0 related to eq.(4.8) of [2] and the value of 1 to the eq. following eq.(4.9), but the second option is strongly discouraged and it will probably be eliminated.

The file *site.case* displays the atom type and coordinates of the single cell, no modification with respect to the file adopted in **lmf**. This file can also be incorporated in the *ctrl.case* file by listing it under the token *SITE*. In the case of magnetic calculations, a *site2.case* with the coordinates related to the super-cell is usually employed instead, with the switch _-vfile=2_to be added to the command line.

The file *rst.case* contains the charge density of the system.

The content of *sigm.case* file is the QSGW self-energy scaled by the exchange-correlation potential in such a way to eliminate this contribution when updating the Hamiltonian. For this reason when no *sigm.file* file is found the exchange-correlation potential will be read automatically and a DFT+DMFT loop will be carried on.

##### indfml File

The *indfml.case* is the key input file in **lmfdmft**, containing all the parameters required to run the DMFT loop starting from QSGW inputs. It was originally part of K. Haule’s Wien2k-DMFT code [1] and therefore it still contains features which are not meaningful in this implementation. It can be generated by the script **init_dmft.py** taking the *ctrl.case* file as input or it can be simply copied from other calculations and revised. An example:

```
#---- Example of indmfl file---
0.1 1.2 1 2 # hybridization Emin and Emax, measured from FS, renormalize for interstitials, projection type
1 0.0025 0.0025 600 -3.000000 1.000000 # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_max (in eV)
1 # number of correlated atoms
1 1 0 # iatom, nL, locrot
2 2 1 # L, qsplit, cix
#================ # Siginds and crystal-field transformations for correlated orbitals ================
1 5 5 # Number of independent kcix blocks, max dimension, max num-independent-components
1 5 5 # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'x^2-y^2' 'z^2' 'xz' 'yz' 'xy'
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 2 0 0 0
0 0 3 0 0
0 0 0 4 0
0 0 0 0 5
#---------------- # Transformation matrix follows -----------
0.70710679 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000
0.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.70710679 0.00000000 0.00000000 0.00000000 -0.70710679 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.00000000 0.70710679 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000 0.00000000
-0.70710679 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000
```

Let us review the content block by block.

```
0.1 1.2 1 2 \# hybridization Emin and Emax, measured from FS, renormalize for interstitials, projection type}
```

The first 2 entries indicate the minimum and maximum values for the energy window entering the projector. These values are calculated from the eigenvalues at the point in correspondence of the bands indicated in the *DMFT_NLOHI* token of the *ctrl.case* file. The 3rd value regards a normalization option for the interstitials which is not taken into account in this code and the 4th entry is equivalent to the token *DMFT_PROJ* of *ctrl.case*.

```
1 0.0025 0.0025 600 -3.000000 1.000000 \# matsubara, broadening-corr, broadening-noncorr, nomega, omega\_min, omega\_max (in eV)
```

The first entry is a switch for Matsubara frequencies (1 for true, 0 for false), but the impurity solver does not handle real frequencies so only the option 1 is advised. The 2nd and 3rd values are broadening parameters not taken into account. The same can be said about the 4th, 5th and 6th entry, controlling # of frequency points, minimum and maximum frequency.

```
1 # number of correlated atoms
1 1 0 # iatom, nL, locrot
2 2 1 # L, qsplit, cix
```

This block defines the local subset of correlated orbitals chosen. First line: the 1st entry indicates the total number of different correlated atoms (equivalent or non-equivalent ones). Second line: there will be a line like this for each of the correlated atoms specified at the line before (just 1 in this case). The index *iatom* points to the specific atom chosen with respect to the legend in the *site.case* file, *nL* to the number of different orbital characters for this atom (e.g. for both , Copper orbitals *nL* would be 2) and is a flag controlling local rotations (1 or 0). Third line: specifications for atom of the previous line. Gives the value of the angular momentum (orbital in this case), *qsplit* refers to the Harmonic orbitals basis in which the projector and hybridization function are written, *cix* indicates the correlated block which this correlated set of orbitals belongs to. If orbitals of different character are included in the same block they will all be listed in the same file for the hybridization function and impurity self-energy. In this case the user will have the possibility of handling matrix elements (referring to the indices of section 4.1 of [2]) of mixed character, say, Oxygen-2 and Copper-3 for what regards the local functions.

```
#================ # Siginds and crystal-field transformations for correlated orbitals ================
1 5 5 # Number of independent kcix blocks, max dimension, max num-independent-components
1 5 5 # cix-num, dimension, num-independent-components
```

This block gives indications regarding the correlated blocks included. The first line is a comment line. The entries at the second and third line define the same quantities except the second line refers to the maximum values for all the correlated blocks included and the second line for the specific block considered (in this case they are equivalent having selected just one block of orbitals). The 3 entries indicate the block index, its dimension (which corresponds to the value of ) and the number of independent components (in the case of degeneracy this number is smaller than the dimension of the matrix).

```
#---------------- # Independent components are --------------
'x^2-y^2' 'z^2' 'xz' 'yz' 'xy'
```

This block is just explanatory and not read in the code. The second line gives the -character order of orbital components of the index, with in the case of -shells.

```
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 2 0 0 0
0 0 3 0 0
0 0 0 4 0
0 0 0 0 5
```

This block refers to the matrix elements of the impurity self-energy contained in the *sig.inp* file. It indicates the matrix of components in the basis following the notation of the previous block. In this example just the diagonal matrix elements are included with no degeneracy. A typical case of degeneracy in 3 systems reduces the total number of components to 3 with the matrix looking like

```
#---------------- # Transformation matrix follows -----------
0.70710679 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000
0.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.70710679 0.00000000 0.00000000 0.00000000 -0.70710679 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.00000000 0.70710679 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000 0.00000000
-0.70710679 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000
```

This last block gives the transformation matrix from the harmonic basis chosen to complex spherical harmonics. It is another feature of Wien2k-DMFT code, not used in. A unitary transformation of this matrix (according to Kristjan’s Haule notation) gives the matrix appearing in the file which on the contrary is employed and will be described later. The reader will find a more extensive example of a indmfl file here in the case of super-cell and different blocks of correlated orbitals, the corresponding super-cell site file is here.

##### sig.inp file

The *sig.inp* file contains the impurity self-energy resulting from the impurity solver. Specifically it is obtained from the file *Sig.out* resulting from the *CTQMC* solver and applying an operation of broadening to it as explained here. It is read in **lmfdmft** to update the local Green’s function and the hybridization function after an operation of embedding of the local impurity self-energy. The structure of this file is the following: a first column listing the values of Matsubara frequencies included (as many values as the token *DMFT_NOMEGA* indicates), and the other columns listing the complex values of in the order and total number indicated in the block of the file (see previous section). An example of the first lines (referred to the convention of the previous section):

```
#---- Example of siginp file (after several iterations)---
# broad_sig.f90 : input parameters = Sig.out 230 l (55 20 230) k (1 2 3 2 3)
0.06283185 78.23564148 -0.04490779 78.62860107 -0.03800205 80.56986618 -0.45687344 78.62860107 -0.03800205 80.56986618 -0.45687344
0.18849556 78.22446442 -0.13346720 78.61967087 -0.11252947 80.72179031 -0.74727110 78.61967087 -0.11252947 80.72179031 -0.74727110
0.31415927 78.20281219 -0.21841954 78.60257339 -0.18401078 80.78148651 -0.94496335 78.60257339 -0.18401078 80.78148651 -0.94496335
0.43982297 78.17199707 -0.29778469 78.57825470 -0.25096695 80.76762009 -1.11261369 78.57825470 -0.25096695 80.76762009 -1.11261369
0.56548668 78.13374329 -0.37015086 78.54793167 -0.31232700 80.72345734 -1.26635915 78.54793167 -0.31232700 80.72345734 -1.26635915
```

The total number of column is 11. The first one lists the frequencies, the other 10 are the complex values (real and imaginary part) of the self-energies related to the 5 components of the -orbitals (as indicated by the *sigind* matrix). At the first iteration the self-energy is null, therefore the only non-zero column will be the first one listing the frequencies. The user will not need to prepare a blank file for the first iteration. When **lmfdmft** is executed a suitable file will automatically be created. The **lmfdmft** command will then have to be called twice in the first iteration, one to create the zero *siginp* file, and a second one for the actual execution of the first loop.

#### The lmfdmt cycle and its main subroutine: sudmft.f

This section will be dedicated to a focus on the **lmfdmft** executable, with a description of its key steps together with the corresponding subroutines. We will refer to sections 4.2 and 4.3 of [2] and in particular Figures 4.1 and 4.2 for reference. The main body of **lmfdmft** is the **sudmft.f** subroutine, this is a list of its main steps (with a reference to the subroutines executing the several operations):

Read the

*indmfl*file and store the parameters needed to build the projector. This operation is carried on by routine.Read the charge density from the

*restart*file, the QSGW self-energy from the*sigm*file (when no*sigm*file is found a LDA calculation is assumed) and diagonalize the QSGW Hamiltonian to return a set of quasi-particles eigenvalues and eigenvectors (see top block of Fig.4.2) in the FP-LMTO basis (introduced in sec. 4.1.1.1). This operation is carried on by routine. It is the first conceptual step of the DMFT loop, repeated for each iteration but with the same results. This is because the QSGW Hamiltonian and eigenvalues are frozen, whereas the other quantities are updated up to DMFT convergence.Read the impurity self-energy file

*siginp*and store the content in an array. This operation is carried on by**readsiginp**routine. At the first iteration*siginp*is null and the code will recognize not to take any double-counting into account.Construct the un-normalized DMFT projection operator (eq. (4.6)), by routine . Build the overlap matrices Olapm and use those to renormalize the projectors (eq. (4.8)), this operation follows the call of

**makeproj**in**sudmft**. The projector is then dumped to disk.Embed the impurity self-energy from the local into the full-space basis using the projector just generated (according to eq.(4.3b)) and remove the DC contribution (eq.(4.23)). This operation is carried on by

**embed_sigma**routine.Tune and update the chemical potential. This operation in quite elaborate, we refer to sec.4.4.3 for details. It requires the diagonalization of the Hamiltonian over all Matsubara frequencies after updating it with the embedded impurity self-energy (routines

**makehbar2**and**agonham**). The following step is to update the valence charge from the corresponding dynamical eigenvalues (routine**cmpvalcharg_matsub4**). Finally determine the correct correction to the chemical potential in such a way to return the correct electron number (routine**getVnew**).Define the local functions by means of the projection operator and the correction from the embedded impurity self-energy. The local Green’s function (eq. (4.24)) is generated by routine

**makegloc2**and the impurity levels (eq. (4.25)) by routine**makeeimp**. These functions are both written to disk.Build the Hybridization function (according to eq. (4.19)) and store it into disk. This operation is carried on by

**makedelta3**routine.

##### New feature: Total QSGW+DMFT charge density

The steps just summarized are the key features of routine, representing the first block of the DMFT cycle (the second block is the object of next chapter), they are therefore repeated for each iteration of the cycle. When DMFT convergence is reached the **lmfdmft** executable can be adopted with a specific flag to extract the total charge density, updated by means of the converged impurity self-energy. This procedure is carried on following the prescription of [1], and it is a recent implementation coded by L.Sponza. For a review of this operation, together with the routines of **lmfdmft** responsible for its execution, we refer to section 6.3 of this manual written by L. Sponza.

#### Output Files

In this section we present an overview of the main output files of the **lmfdmft** cycle. They will be the contents of the so-called *it#_lmf* folder, where the # stands for the given iteration. The iterations following the first one use as input the impurity self-energy resulting from the *CTQMC* solver. Here we present the simple case of the output files of the first iteration.

```
#---- Output of lmfdmft run (in chronological order)---
evec.case
wkp.case
moms.case
proj.case
delta.case
eimp1.case
gloc.case
log.case
log
```

The first three files are regular outputs of **lmf**, binary files containing respectively the eigenalues (*evec*), information regarding the -point mesh (*wkp*) and momentum (*moms*). They are temporary inputs used to produced the other outputs and they are not taken as inputs in the following steps of the DMFT cycle. The *proj.case* file is another binary output storing the DMFT projection operator. It is generated and stored to disk on the fly and then read to define the local functions. The *gloc.case* is instead an ASCII file, containing the local correlated Green’s function as a function of Matsubara frequencies. The file structure is the same as for the *sig.inp* file, with as many columns as the correlated orbitals selected in the *Sigind* block of the *indmfl* file. The *gloc.case* is not taken as input by other routines but is a useful reference for other operations (such as the analytic continuation). The files *delta.case* and *eimp1.case* are ASCII files representing the key output of **lmfdmft**. They are taken as input by the *CTQMC* impurity solver. What was said about the *gloc.case* file’s structure is equivalent for *delta.case*, storing the hybridization function in the correlated local basis. The *eimp1.case* file contains instead the impurity levels and the double-counting constant used. The file has a composite structure tailored to be parsed by the scripts of the interface in the *CTQMC* solver. An example:

```
#---- Example of eimp1.case file---
Edc=[ 25.500000, 25.500000, 25.500000, 25.500000, 25.500000, 25.500000, 25.500000, 25.500000, 25.500000, 25.500000] # Double counting
Eimp=[ 0.326716, 0.393662, 0.061104, 0.460181, 0.061165, 0.326705, 0.393612, 0.061085, 0.460140, 0.061127] # Eimp: P.(e+E.sbar)-sinp-mu+Edc
Eimp=[ 0.000000, 0.066946, -0.265612, 0.133464, -0.265551, -0.000011, 0.066895, -0.265631, 0.133423, -0.265589] # Eimp: shifted by Eimp[0]
Ed [ -25.173284, -25.106338, -25.438896, -25.039819, -25.438835, -25.173295, -25.106388, -25.438915, -25.039860, -25.438873] # Ed: Eimp-Edc for PARAMS
mu 25.173284 # mu = -Eimp[0] for PARAMS
```

The *log* file shows the several passages of the routine summarized in the previous section. Let us review a portion of its content (in the case of the first iteration for zero *sig.inp*).

```
#---- First block of log file---
iors : read restart file (binary, mesh density)
use from restart file: ef window, positions, pnu
ignore in restart file: *
```

from which we know that the *restart* file containing the charge density has been correctly read.

```
#---- Second block of log file---
... sudmft job=1: make projectors
Reading indmfl file ... 1 total (1 inequivalent) correlated blocks, among 1 sites
readindmfl: 5 nonzero (5 inequivalent) matrix elements
channels for 2nd spin block
m1 m2 cix chn(1) chn(2)
1 1 1 1 6
2 2 1 2 7
3 3 1 3 8
4 4 1 4 9
5 5 1 5 10
correlated channels:
chan equiv m1 m2 isp icix cix ib
1 1 1 1 1 1 1 1
2 2 2 2 1 1 1 1
3 3 3 3 1 1 1 1
4 4 4 4 1 1 1 1
5 5 5 5 1 1 1 1
6 6 1 1 2 1 1 1
7 7 2 2 2 1 1 1
8 8 3 3 2 1 1 1
9 9 4 4 2 1 1 1
10 10 5 5 2 1 1 1
indmfl file expects 1 DMFT block(s) (max dimension 5), 10 matrix elements
```

This section corresponds to the creation of projection operators after reading the *indfml* file. It gives direct infos about the correlated blocks of orbitals chosen and their character. In this case we notice how 10 channels are produced from a block of 5 in the *indmfl* file since the calculation is magnetic and the spin is 2.

```
#---- Third block of log file---
Reading DMFT sigma from file ...
DMFT sigma is zero ... no double counting
SUDMFT: projector #2 with k-integrated norm. 8 bands in window (4,11)
2000 Matsubara frequencies, interval (0.0628318,251.265) eV
User-specified double-counting in eV : --ldadc=25.5
```

First part: this being the first iteration, the file *sig.inp* is null. The code recognizes it and no DC is included. Second part: general information about the projector before these are renormalized, the Matsubara mesh, the energy window chosen and the DC calculated according to (4.29) of [2] by reading the values specified as inputs for Hubbard parameters and the nominal occupancy of the correlated orbitals .

```
#---- Fourth block of log file---
BNDFP: Write evals,evecs to file for 29 qp
bndfp: kpt 1 of 29, k= 0.00000 0.00000 0.00000
-3.8973 -3.8973 -3.8973 -0.5804 -0.1352 -0.1352 -0.1352 -0.0280 -0.0280
bndfp: kpt 1 of 29, k= 0.00000 0.00000 0.00000
-3.8973 -3.8973 -3.8973 -0.5804 -0.1352 -0.1352 -0.1352 -0.0280 -0.0280
```

Here a first portion of the eigenvalues generated by routine and stored in *evec.case* is reported. Being a magnetic calculation we notice the same -points contribution being listed twice (for majority and minority spin components). In this case the eigenvalues are identical since the magnetic moment is initially null.

```
#---- Fifth block of log file---
Renormalize projectors ...
Find chemical potential mu...
Seek mu for 14 electrons ... 13.842564 electrons at Ef0=-0.026874. D(Ef0)=57.196
getVnew: v1=0e0 N1=-1.574e-1 v2=-2.753e-3 N2=-5.867e-3 bracket=F est=-0.00286977
getVnew: v1=-2.753e-3 N1=-5.867e-3 v2=-2.87e-3 N2=1.382e-4 bracket=T est=-0.00286707
getVnew: v1=-2.87e-3 N1=1.382e-4 v2=-2.867e-3 N2=-4.814e-8 bracket=T est=-0.00286707
mu = -0.024006 = Ef0--0.002867. Deviation from neutrality = -9.59e-14
Electron charge: 14.000000 moment: 0.000000 spin resolved DOS: 25.553 25.553
```

After the projector are renormalized, the chemical potential finder is run. In this case 3 iterations are sufficient to converge to the correct value for the electron count, the correction to is stored.

```
#---- Sixth block of log file---
Make gloc, delta, eimp ...
Writing files delta and eimp1 ...
Check SC condition skipped (missing information about previous iteration)
gloc(N) recorded in gloc.ext file.
Exit 0 done making DMFT hybridization function
```

Last block of the *log* file. The local functions are generated and written to disk. Being this the first iteration no DMFT self-consistent condition can be met.

### Extracting the Impurity Self-Energy - CTQMC Software

This section will be dedicated to the second (and last) block of the DMFT cycle, the employment of the QSGW-based hybridization function to extract the impurity self-energy of DMFT by means of the impurity solver (see sec.31 of [2] for the details of the equations and main quantities in play). In practice this operation pivots on two main programs, the *CTQMC* (Continuous Time quantum Monte-Carlo) impurity solver and the interface connecting this program with the other block (i.e. the **lmfdmft** routine) completing the DMFT cycle. The *CTQMC*, mainly written in C++ code, is a program developed by K. Haule as part of his Wien2k-DMFT code [1] and which has been adapted in this implementation. The next section will mainly focus on the interface between the **lmfdmft** routine and the solver, the main inputs and scripts required, and a brief overview of the output files. This chapter is complementary to this tutorial by L. Sponza which provides detailed instructions on how to run the *CTQMC* software, with some more insight into the scripts and some tips on how to set the correct parameters.

#### Input Files and Scripts

In this section we will provide a summary of the overall input files required for both the interface and the *CTQMC* software, together with the main scripts needed. The folder containing the main input files and scripts is generally called *qmcinput*. In order to be consistent with the notation of the previous chapter we will refer to *it#_qmcrun* for the folders containing several iterations of *CTQMC*. It is useful to review the input files for both folders.

```
$ ls -l qmcinput
atom_d.py
broad_sig.f90
Trans.dat
PARAMS
$ ls -l it#_qmcrun
#output of atom_d.py
info_atom_d.dat
actqmc.cix
new.cix
#output of previous lmfdmft run
Delta_in -> ../it#_lmfrun/delta.fe
Delta.inp
Eimp_in -> ../it#_lmfrun/eimp1.fe
Eimp.inp
```

The files contained in *qmcinput* are frozen for all the iterations whereas the ones in *it#_qmcrun* are updated for each iteration.

##### The PARAMS File

The *PARAMS* file is the key input file of *CTQMC*. The correct choice for the parameters listed in this file can affect quite significantly the quality of the calculation. For this reason, the user is highly encouraged to read this tutorial on how to set the most important parameters. An example:

```
#---- Example of PARAMS file---
Ntau 1000
OffDiagonal real
Sig Sig.out
Naver 100000000
SampleGtau 1000
Gf Gf.out
Delta Delta.inp
cix actqmc.cix
Nmax 700 # Maximum perturbation order allowed
nom 150 # Number of Matsubara frequency points sampled
exe ctqmc # Name of the executable
tsample 50 # How often to record measurements
nomD 150 # Number of Matsubara frequency points sampled
Ed [ -25.173284, -25.106338, -25.438896, -25.039819, -25.438835, -25.173295, -25.106388, -25.438915, -25.039860, -25.438873] # Impurity levels updated by bash script
M 75000000.0 # Total number of Monte Carlo steps
Ncout 200000 # How often to print out info
PChangeOrder 0.9 # Ratio between trial steps: add-remove-a-kink / move-a-kink
CoulombF 'Ising' # Ising Coulomb interaction
mu 25.173284 # QMC chemical potential by bash script
warmup 100000 # Warmup number of QMC steps
GlobalFlip 200000 # How often to try a global flip
OCA_G False # No OCA diagrams being computed - for speed
sderiv 0.01 # Maximum derivative mismatch accepted for tail concatenation
aom 3 # Number of frequency points used to determin the value of sigma at nom
HB2 False # Should we compute self-energy with the Bullas trick?
U 5.0
J 0.8
nf0 6.0
beta 50.0
```

At the bottom of the file we notice the physical parameters, the Hubbard parameters , the inverse temperature and the occupancy of the correlated orbitals (set in accordance with the values previously selected). In addition, we notice the two lines copied from the *Eimp.inp* file, one labelled by *Ed* (4th line of the file), and one labelled by *M* (last line of the file). All the other parameters control the stochastic character of the Monte-Carlo calculation. For more details about these, and some tips on how to tune them, in order not to be redundant we refer again to L.Sponza’s tutorial.

##### Treating d-systems: the script atom\d.py and the file Trans.dat

The *atom_d.py* is a python script written by K. Haule. It is specifically set for -systems (not adopted otherwise) and its purpose is to initialize the atomic problem for the specific atom considered and to transform the Coulomb interaction via a rotation to the correct harmonics basis. The input file needed is called *Trans.dat*, it contains information about the local basis and in particular the transformation matrix from standard spherical Harmonics to a given user’s specified basis which is more convenient for the calculation. There are basically just two main cases of the file, the user can find an example at the following links, one for non-magnetic and one for spin-polarized calculations. The execution line for *atom_d.py* requires some specific flags indicating parameters such as angular momentum and the starting impurity levels (the correct instructions on how to run it can be found at this page ). The main output file of the script is file *actqmc.cix*. This file has information regarding the atomic basis and the matrix elements of and is used as an input by the *CTQMC* solver and has to be copied in the *it#_qmcrun* folder. Two log files are produced by the script, *info_atom_d.dat* and *new.cix*, neither is required by *CTQMC*.

##### The Delta.inp and Eimp.inp Input Files

These files are soft linked and renamed from *it#_lmf* folder (containing the results of the iteration preceding the *CTQMC* run). The *Delta.inp* is directly taken as an input by *CTQMC*, whereas the *Eimp.inp* is parsed and taken as input by 2 scripts composing the interface. In the first place the 3rd line, namely (see previous chapter for comparison):

```
Eimp=[ 0.000000, 0.066946, -0.265612, 0.133464, -0.265551,-0.000011, 0.066895, -0.265631, 0.133423, -0.265589]
```

is copied in the execution line to run the script, object of next section. The fourth and fifth lines,

```
Ed [ -25.173284, -25.106338, -25.438896, -25.039819, -25.438835,
-25.173295, -25.106388, -25.438915, -25.039860, -25.438873] # Ed
: Eimp-Edc for PARAMS
mu 25.173284 # mu = -Eimp[0] for PARAMS
```

are instead used in the file (introduced further on).

##### The broad\sig.f90 Script

This is the only script treating the output files. It takes the *Sig.out* file, output of *CTQMC* containing the impurity self-energy, and it applies a Gaussian broadening to it in order to lessen the statistical noise. Every channel of the self-energy (i.e. its orbital component) is convoluted with a Gaussian distribution with a frequency-dependent width. This script allows the user to specify the Gaussian width, the total number of frequency points of the broadened self-energy and the starting and ending frequency, depending on the noise and shape of the starting input function. The fortran routine **broad_sig.f90** has obviously to be compiled to obtain the *broad_sig.x* executable. For instructions on how to compile the program and run the executable see this page of the tutorial. The resulting broadened self-energy is recorded in file *Sig.out.brd* wheres files *broad.log* and *width.log* respectively list the parameters used in the broadening and the values of the Gaussian width as a function of the frequency.

#### Outputs of CTQMC

```
#Outputs of CTQMC run
Probability.dat
g_qmc.dat
g_hb0.dat
s_hb1.dat
g_hb1.dat
Gtau.dat
histogram.dat
Sig.out
Gf.out
ctqmc.log
statusfiles/
```

We can differentiate these files into 2 categories. In the first category are files which are useful to judge the quality of the calculation, and among these enter all the **.dat* files and the log file *ctqmc.log*. For a preliminary discussion on these files see this tutorial. A special remark must be undertaken for the so called *statusfiles*. For any *CTQMC* run, which is parallelized on a given number of cores, status files named *status.#* will be produced with information about the sampling. In order to improve the quality of the calculation, these files have to be copied to the folder containing the following step of *CTQMC* in such a way to be read by the program at the start. In the case a single run of *CTQMC* is not enough to produce a sensible impurity self-energy, which might be too affected by statistical noise, the program can be re-run with the same instructions and inputs and the status files just produced will be also read to improve the accuracy of the sampling. To the second group belong the actual output files of the calculation. The file *Gf.out* contains the impurity Green’s function and the file *Sig.out* the impurity self-energy, both files affected by the statistical noise of the Monte-Carlo sampling. In particular the Monte-Carlo sampling accurately accounts for the low frequency region of the self-energy, whereas it is known to be extremely noisy for the high-frequency regime. In order for the self-energy to approach its Hartree-Fock value in the high-energy limit, the high-frequency tails are analytically corrected and concatenated with the Monte-Carlo sampling according to boundary conditions on value and slope. The *CTQMC* program makes use of the Hubbard I approximation for the tails. To get rid of statistical noise the *Sig.out* undergoes the *broad_sig.x* broadening to return the *Sig.out.brd* file, which gets soft linked to the *sig.inp* file entering the following iteration in such a way to restart the cycle.

#### Bonus Track: A Fully Inclusive Interface

At this point the user should have all the means to run a full QSGW+DMFT loop from scratch using the QUESTAAL program, assisted by this manual and following L. Sponza’s tutorials step by step. If this is the case, an inclusive set of scripts will allow the user to run a full calculation in one go by setting all the required parameters in advance. These scripts automatize the several operations needed to connect the **lmfdmft** block to the *CTQMC* one in such a way as to build a more user-friendly interface. They have been written mainly by L.Sponza and they are organized as follows. A central script **Contlte.sh** reads a file called *LoopParams.sh* with all the overall parameters needed to set the DMFT calculation (provided that all of the input files discussed in this manual have already been set up in the respective input folders *qmcinput*, *lmfinput* ) and prepares the folders hosting the calculation of the two main programs. After that, it subsequently submits the corresponding scripts, which also have to be prepared in advance into the respective input folders, containing the execution lines for the two programs. A self-explanatory example of the main script **Contlte.sh** can be found at this link, together with an example of the *LoopParam.sh* file. An example of the two scripts, one for the **lmfdmft** run (which in this version runs on a single processor) and one for the *CTQMC* run (which runs in parallel on a given number of cores) is also provided.

### References

[1] K. Haule, C.H. Yee, and K. Kim. Dynamical mean field theory within the full-potential methods: Electronic structure of ceirin5, cecoin5, and cerhin5. Phys. Rev. B, 81:195107, 2010. [2] Paolo Pisanti. A novel QSGW+DMFT method for the study of strongly correlated materials. PhD thesis, King’s College London, 2016.