SciPost Phys. 10, 117 (2021)

DCore: Integrated DMFT software for correlated electrons
Hiroshi Shinaoka1? , Junya Otsuki2 , Mitsuaki Kawamura3 ,Nayuta Takemori2 and Kazuyoshi Yoshimi3
1 Department of Physics, Saitama University, Saitama 338-8570, Japan
2 Research Institute for Interdisciplinary Science, Okayama University,
Okayama 700-8530, Japan
3 Institute for Solid State Physics, University of Tokyo, Chiba 277-8581, Japan? shinaoka@mail.saitama-u.ac.jp

Abstract
We present a new open-source program, DCore, that implements dynamical mean-field
theory (DMFT). DCore features a user-friendly interface based on text and HDF5 files.
It allows DMFT calculations of tight-binding models to be performed on predefined lattices as well as ab initio models constructed by external density functional theory codes
through the Wannier90 package. Furthermore, DCore provides interfaces to many advanced quantum impurity solvers such as quantum Monte Carlo and exact diagonalization solvers. This paper details the structure and usage of DCore and shows some
applications.
Copyright H. Shinaoka et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 08-07-2020
Accepted 17-05-2021
Check forupdatesPublished 27-05-2021
doi:10.21468/SciPostPhys.10.5.117

Contents
1 Introduction

1

2 Methodology
2.1 One-shot DFT+DMFT calculations and tight-bonding models
2.2 Model
2.3 DMFT formalism
2.4 MPI parallelization

2
2
3
3
5

3 Structure of DCore and flowchart of simulation

6

4 Format of input and output files
4.1 Input file format
4.2 Output-file format

7
7
10

5 Installation

11

6 Examples
6.1 First example: Square-lattice Hubbard model

11
11

1

SciPost Phys. 10, 117 (2021)

6.2

6.3

6.1.1 Model construction
6.1.2 Self-consistent calculations
6.1.3 Dynamical quantities
6.1.4 CT-QMC solver
6.1.5 Wannier90 interface
6.1.6 Antiferromagnetic state
Solving multi-orbital model using QMC solver: t 2g model on a Bethe lattice
6.2.1 Model definition
6.2.2 DMFT calculation
DFT+DMFT example: SrVO3
6.3.1 Construction of Wannier functions
6.3.2 Model definition
6.3.3 Parameters for self-consistent calculations and post processing
6.3.4 Results

12
12
14
15
18
19
21
21
22
23
24
24
24
25

7 Conclusion

26

References

27

1

Introduction

Dynamical mean-field theory (DMFT) has become a standard theoretical tool for studying
strongly correlated electronic systems [1]. In a DMFT calculation, an original lattice model ismapped to an effective Anderson impurity problem whose bath degrees of freedom are selfconsistently determined. Although the DMFT formalism was originally proposed for models
such as Hubbard models, it can be combined with density functional theory (DFT) based on abinitio calculations to describe the electronic properties of strongly correlated materials [2, 3].
This composite framework, called DFT+DMFT, has been applied to various types of materials,
such as cuprates [4], Fe-based superconductors [5–7], and f -electron materials [8–10].
To make DFT+DMFT available to more users and to promote the development of a community of users and developers, there is a strong need for an open-source program package witha user-friendly interface. Some open-source software packages for performing DFT+DMFT
calculations have been developed. The TRIQS (Toolbox for Research on Interacting Quantum
Systems) project [11] aims to provide the necessary components to implement DFT+DMFT
codes, including a quantum Monte Carlo (QMC) impurity solver [12] and an interface to DFT
codes [13]. The user can implement their own DMFT code using the building blocks in Python.
The ALPS project provides a DMFT code for simple Hubbard models [14] and implementationsof several continuous-time QMC (CT-QMC) impurity solvers [15–17]. w2dynamics [18] provides a state-of-the-art implementation of a QMC algorithm and includes a Python program for
performing simple DFT+DMFT calculations. DMFTwDFT features a user-friendly interface andhas interfaces to various DFT codes and the CT-QMC impurity solver [19]. DFT+embedded
DMFT Functional (eDMFT) [20] features a full set of DFT+DMFT functionalities but it is notopen-source software.
DCore is an open-source program package that implements (DFT+)DMFT calculations for
multi-orbital systems. This package is built on top of existing TRIQS Python libraries and provides interfaces to external impurity solvers from the ALPS [15,16] and TRIQS projects [12,21]
as well as exact-diagonalization solvers based on pomerol [22] and HΦ [23]. DCore provides
2

SciPost Phys. 10, 117 (2021)

a flexible interface based on text and HDF5 files and does not require extensive expert knowledge. It also provides a set of command line tools to analyze a self-consistent solution andevaluate physical quantities such as momentum-resolved spectral functions.
The remainder of this paper is organized as follows. In Sec. 2, the methods used for
implementing DCore are described. In Sec. 3, the structure of DCore and a flowchart ofsimulations are shown. In Sec. 4, the formats of input and output files are explained. In Sec.
5, brief instructions on the installation of DCore are provided. In Sec. 6, several run examplesof DCore from simple Hubbard models to DFT+DMFT calculations for a correlated material
are demonstrated. In Sec. 7, a summary and the conclusions of this paper are given.

2

Methodology

In this section, a brief explanation of the methodology used for implementing DCore is provided. For more technical details of the DFT+DMFT method, we refer the reader to a review
article [2]. The current version of DCore implements only one-shot DFT+DMFT calculations
based on Wannier901 . DCore does not support charge self-consistent calculations [9], which
are supported by other programs such as eDMFT and DMFTwDFT. In a charge self-consistentcalculation, the DFT effective potential is updated using the density matrix obtained by DMFT
calculations. This is essential especially in performing lattice optimization. DCore specializes
in the analysis of multi-orbital Hubbard models.

2.1

One-shot DFT+DMFT calculations and tight-bonding models

In one-shot DFT+DMFT calculations, one first performs “band calculations” usually using
the local density approximation (LDA), where correlation effects in strongly correlated orbitals (e.g., partially filled d orbitals) are not described accurately. One projects the effective
one-body LDA Hamiltonian onto some tight-binding basis, which yields a multi-orbital tightbinding model. Then, one introduces local electron interactions for correlated orbitals into
the tight-binding model. The resultant multi-orbital Hubbard model is solved using DMFT. We
will discuss how to subtract the contributions of electron correlations at the LDA level fromDMFT results (double-counting corrections) in a later section. One may not introduce additional interactions at the DMFT level for weakly correlated orbital for which LDA is accurateenough (e.g., deep oxygen orbitals).

2.2

Model

DCore is designed to solve a tight-binding model with periodic boundary conditions in second quantization. We denote the primitive vectors of the lattice by a1 , a2 , a3 . Then, thecoordinates of a unit cell can be specified by an integer array of length 3, R = (R1 , R2 , R3 ), asP
3
i=1 R i a i . Each unit shell involves shells, which are sets of spin orbitals. We define a creation†
operator and an annihilation operator for the α-th spin orbital in the s-th shell as cRsα
and cRsα ,respectively. In the reciprocal space, we define the corresponding creation and annihilation
operators by
1 X ik·R †
†
cksα
=p
e cRsα ,N R
1 X −ik·R
cksα = pecRsα ,N R
1

DCore uses a Wannier90 converter in TRIQS/DFTTools [13].

3

(1)
(2)

SciPost Phys. 10, 117 (2021)

Figure 1: Illustration of the DMFT self-consistency cycle defined in Eqs. (6)–(10).
The most numerically expensive step is the solution of the impurity problem, i.e.,the calculation of the self-energy for a given G0−1 . 〈. . .〉k denotes an average over the
momentum space and µ is the chemical potential.
where momentum index k = (k1 , k2 , k3 ) and ki ∈ Z. N denotes the number of k points.
The Hamiltonian of a tight-binding model can be written as follows:

H=
=

N
X

X

X

k=1 s,s0 ∈all shells αβ
N
X

X

X
R,s∈crsh

X

R,R0 s,s0 ∈all shells αβαβ†Hss0 (k)cksα
cks0 β +αβ†Hss0 (R − R0 )cRsα
cR0 s0 β +Hint (R, s),

(3)

X

(4)

R,s∈crsh

Hint (R, s),where the first term in Eq. (3) denotes a non-interacting Hamiltonian. In the above equations,s and s0 run either over correlated (crsh) “shells” or over all “shells” (i.e., the correlated shells
and a non-interacting shell2 ). Hereafter, the number of the correlated shells is denoted byncor. α and β denote the spin-orbital index in each shell. The second term Hint (R, s) denotes
an interaction at the s-th correlated shell. This Hamiltonian is given as

Hint (R, s) =1 X††
Uαβγδ (s)cRsα
cRsβ
cRsδ cRsγ ,2 αβγδ

(5)where Uαβγδ (s) denotes a spin-full four-rank Coulomb tensor. Note that Hint is limited to
intra-shell interactions and local interactions in each unit cell.

2.3

DMFT formalism

DCore solves the following DMFT self-consistent equations, where the lattice Green’s function
G, the local Green’s function Gloc , the noninteracting impurity-model Green’s function G0 , andthe impurity/lattice self-energy Σimp /Σlattice are defined in the spin-orbital space:−1G(k, iωn ) = iωn + µ − H(k) − Σ(iωn ) + ΣDC,(6)
Gloc (iωn , s) = Ps 〈G(k)〉k (∀s ∈ crsh),G0−1 (iωn , s) = [Gloc (iωn , s)]−1 + Σ(iωn , s) (∀s ∈ crsh),imp
Σσσ0 (iωn , s) ← G0 (iωn , s), Uαβγδ (s) (∀s ∈ crsh; σ, σ0 ∈ {↑, ↓}),
imp
Σσσ0 (iωn , s) ← Σσσ0 (iωn , s).
2

(7)
(8)
(9)
(10)

The term “non-interacting shell” means that there are no additional electron interactions introduced at the
DMFT level for the shell.

4

SciPost Phys. 10, 117 (2021)

Here, ωn ≡ (2n+1)πT is a fermionic Matsubara frequency at temperature T , 〈· · ·〉k denotes anaverage over the momentum space, s indexes correlated shells, and Ps is the projector to thes-th correlated shell. Double-counting corrections are denoted by ΣDC in Eq. (6) [see Eq. (15)].
Equations (6)–(10) are illustrated in Fig. 1.
All Green’s functions and self-energies in Eqs. (6)–(10) are matrices at each Matsubara
frequency, and their inversions are regarded as inverse matrices. The full self-energy Σ(iωn )in Eq. (6) is defined for all spin and shell components, and has the spin-block structureΣ↑↑ (iωn ) Σ↑↓ (iωn )Σ(iωn ) =,(11)Σ↓↑ (iωn ) Σ↓↓ (iωn )where each spin-component has the internal shell structureΣσσ0 (iωn , s = 0) 0
···
0
.. 
..

.
0
.

Σσσ0 (iωn ) = 
,
.

..
Σσσ0 (iωn , s = ncor − 1) 0
0
···
0
0

(12)

where the last column and row correspond to the non-interacting shell. The block structuresof the Green’s functions are determined consistently with those of the self-energies in Eqs. (6)–(10).
The Green’s functions and the self-energies are assumed to be either spin-diagonal or dense
matrices in the spin space with both diagonal and off-diagonal components. The former issufficient for computing a paramagnetic state or a collinear magnetic state polarized along thez axis. The latter allows Σ(iωn ), G(iωn ), and H(k) to have spin-off-diagonal elements. This
must be used when H has spin-orbit coupling that mixes the spin-up/down sectors or when
one computes non-collinear magnetic states.
As illustrated in Fig. 2, by default, all correlated shells are assumed to be inequivalentshells that can have different self-energies. As described in the following sections, the user canassume that several correlated shells have the same self-energy. In this case, these correlatedshells define one inequivalent shell (see Case 2 in Fig. 2).
The chemical potential µ is either adjusted to obtain a desired number of electrons or
is fixed at a given value. In the former case, the value of µ is determined by using a bisection algorithm, which requires solving Eq. (6) multiple times at each self-consistent iteration. Note thatPthe number of electrons per unit cell is given by −Tr G(τ = β − 0+ ),∞where G(τ) ≡ β −1 n=−∞ e−iωn τ 〈G(k, iωn )〉k (β and Tr denote inverse temperature and matrix trace, respectively).
In Eq. (10), the new self-energy is computed by solving an impurity model for each inequivalent shell. The non-interacting and interacting parts of the impurity model are definedby G0 (iωn , s) and Uαβγδ (s), respectively. Some QMC impurity solvers accept the hybridization
function defined below as an input:∆(iωn , s) = iωn + µ − Ps 〈H(k)〉k − G0−1 (iωn , s).

(13)

After computing impurity self-energies Σimp (iωn ) for all inequivalent shells, we update Σ(iωn , s)
by using simple linear mixing as
Σnew (iωn , s) ← σmix Σimp (iωn , s) + (1 − σmix )Σold (iωn , s)

(14)

and go back to Eq. (6). Here, σmix (0 < σmix ≤ 1) is a mixing parameter.
DCore implements three different algorithms for computing the double-counting correction ΣDC :5

SciPost Phys. 10, 117 (2021)

Case 1 (default)
Corr. shell 0

Corr. shell 1

Corr. shell 2

Corr. shell 3

Inequivalent shell 0

Inequivalent shell 1

Inequivalent shell 2

Inequivalent shell 3

Corr. shell 2

Corr. shell 3

Inequivalent shell 1

Inequivalent shell 2

Case 2

Assumed the same self-energy

Corr. shell 0

Corr. shell 1

Inequivalent shell 0

Figure 2: Correlated shells and inequivalent shells• HF_DFT (default): Hartree-Fock (HF) contribution is subtracted from the DFT input.
• HF_imp: HF contribution computed from the local Green’s function Gloc [Eq. (7)] is
subtracted from the impurity self-energy.
• FLL: A formula called fully-localized limit [24] is used.
For the default algorithm (HF_DFT), the double-counting correction is computed as follows:X †csδ 〉0 ,(15)Uαγβδ (s) − Uαγδβ (s) 〈csγΣDC
αβ (s) =γδ

where 〈· · · 〉0 indicates the expectation value at the initial (Kohn-Sham) state with k-summation
included, and {α, β, γ, δ} index spin orbitals.

2.4

MPI parallelization

DCore uses MPI parallelization for performing numerically expensive parts of DMFT calculations. For instance, k-average of G(k, iωn ) in Eq. (6) is executed in parallel. Besides, some
impurity solvers, such as CT-QMC solvers, may use MPI parallelization. As described in thefollowing sections, all programs in DCore are launched as a single process. Some of theminternally invoke MPI processes. For some examples, please refer to Sec. 6.

3Structure of DCore and flowchart of simulation

Input parameter
Temperature,
numerical conditionetc.

: External code
: Small input file: Executable in this package: Data file
OpenMX

openmx2dcore

Input parameter
Correlated shell
U and J term
Input parameterk gridfor post processes
etc.
Loaded by dcore_pre, dcore,dcore_check, and dcore_post

Hopping integral(wannier90 format)
DFT codes(QE, ABINIT, SIESTA,FLEUR, VASP, WIEN2K, etc.)

dcore_pre

Impurity solvers
ALPS/CT-HYB
ALPS/CT-HYB-SEGMENTTRIQS/cthyb
TRIQS/Hubbard-I
pomerol, H©

DMFT configuration(HDF5 file)

Wannier90

dcore

dcore_check

DMFT result(HDF5 file)dcore_post

Figure 3: Structure of DCore.

6

SciPost Phys. 10, 117 (2021)

DCore contains a set of programs that perform DMFT calculations for models and materials.
The structure of the programs and data flow is summarized in Fig. 3.
DCore consists of four layers: (i) interface layer, (ii) DMFT self-consistent loop, (iii) convergence check and (iv) post-processing. These are respectively performed by the executables
dcore_pre, dcore, dcore_check and dcore_post. Input parameters are provided by a single text file, which is read by all three programs. Data generated by dcore_pre and dcore
are separately stored in an HDF5 file and passed to the next process.
(i) Interface layer: dcore_predcore_pre generates the HDF5 file necessary for the DMFT loop. Users specify parameters that define a model, such as the hopping parameters on a certain lattice, andinteractions. The hopping parameters are given either for preset models (e.g., squarelattice, Bethe lattice) or using the Wannier90 format.
(ii) DMFT self-consistent loop: dcore
dcore is the main program for the DMFT calculations. The effective impurity problem
is solved repeatedly to satisfy the self-consistency condition of the DMFT. For solving
the impurity problem, dcore calls an external program. As an external program, we canselect ALPS/CT-HYB [15], ALPS/CT-HYB-SEGMENT [16], or TRIQS/cthyb [12] as a CTQMC solver, TRIQS/Hubbard-I [21] as a Hubbard-I solver, or pomerol [22] or HΦ [23]as an exact diagonalization (ED) solver.
(iii) Convergence check: dcore_checkdcore_check reads the results of self-consistent calculations (e.g., the self-energy) fromthe output HDF5 file of dcore and outputs several figures that are useful for checking
the convergence of the results.
(iv) Post-processing: dcore_post
dcore_post computes some physical quantities from the converged solution of the
DMFT loop. Currently, the (projected) density of states A(ω) and correlated band structures (momentum-resolved single-particle excitation spectrum) A(k, ω) can be calculated. The analytic continuation of the self-energy from the Matsubara frequencies toreal frequencies is performed using the Padé approximation.

4

Format of input and output files

In this section, the format of input and output files for DCore is briefly explained.

4.1

Input file format

The input file of DCore consists of six parameter blocks, respectively named [model], [system],[impurity_solver], [control], [tool], and [mpi]. All the programs usually read thesame file, but that some blocks are, even if present, irrelevant to some of the programs (seeTable 1). Each block is described below.
1. [model] block
This block includes parameters for defining a model to be solved. The parameter types
are divided into four parts: (i) Basic parameters, (ii) Lattice parameters, (iii) Interaction parameters, and (iv) Local potential parameters. In the following, we describe theparameters defined in each part.
7

SciPost Phys. 10, 117 (2021)

Table 1: Blocks to be used for each program in DCore.
Block

model
systemimpurity_solver
control
tool
mpi

dcore_predcore

dcore_checkdcore_post

Yes
-

Yes
Yes
Yes
Yes
Yes

Yes
Yes
Yes
-

Yes
Yes
Yes
Yes
Yes

(i) Basic parameters
There are four basic parameters, namely seedname, nelec, norb, and spin_orbit.
seedname determines the name of the model HDF5 file. nelec and norb specify the
number of electrons per unit cell and the number of orbitals, respectively. spin_orbit
specifies whether the model has spin-orbit interactions.
(ii) Lattice parameters
In DCore, chain, square, cubic, and bethe lattices are prepared as predefined models, as shown in Fig. 4. We can select a lattice type using the lattice parameter andset values of transfer integrals using parameters t and t’. For treating more realistic
cases, such as DFT+DMFT calculations, hopping parameters in the Wannier90 format
can be imported by selecting wannier90. In this mode, we can set a number of correlated shells in a unit cell and a mapping from correlated shells to equivalent shells using
the ncor and corr_to_inequiv parameters, respectively. Figure 5 shows a schematic
diagram of the shell structure in wannier90 mode. For experts, the lattice data can becustom prepared using the external mode. In this mode, all necessary data should
be directly made in the dft_input group of the model HDF5 file. For details, see thereference manual of DFTTools [13]. The information of reciprocal lattice vectors andthe number of wave vectors can be also specified.
s (iii) Interaction parameters
The interaction part of the Hamiltonian is defined as Eq. (5). The interaction matrix Uαβγδ (i) is specified by the parameter interaction. In DCore, three types ofinteraction, namely (a) kanamori, (b) slater_f, and (c) slater_uj, are defined.
(a) If interaction = kanamori, the Kanamori-type interaction is used; i.e.,U(aσa )(bσ b )(cσc )(dσd ) = Va bcd δσa σc δσ b σd with Vaaaa = U, Vabab = U 0 , Vabba = J, Vaabb = J,
where a 6= b (a, b, c, d are orbital indices, and σa , σ b , σc , σd are spin indices). These(a) chain(b) square(c) cubic(d) betheDensity of states

t

¡2t

t'

Energy2tFigure 4: Predefined lattices in DCore. In the (a) chain , (b) square, and (c)cubic lattices, the solid (dotted) bonds indicate the transfer integrals t (t 0 ). For
the (d) bethe lattice, a semicircular density of states with energy ranges [−2t:2t] isdefined.

8

SciPost Phys. 10, 117 (2021)(b) spin_orbit=True(a) spin_orbit=False (default)[model]ncor = 2
norb = 3
corr_to_inequiv = 0, 0
spin_orbit = False

Corr. shell1(inequivalent shell 0)<latexit sha1_base64="69B/Uk05O8z1FK2i5guc56ZfHuc=">AAACx3icbVHLbtQwFPWEVxleLV2ysRhVYoGqBIFgwaISG9i1gplWmkSV49ykpn5E9s0wIysL/oAt/AJf1L+pk8miM+2VbB2fe659fG9eS+Ewjq9G0b37Dx4+2nk8fvL02fMXu3svZ840lsOUG2nsWc4cSKFhigIlnNUWmMolnOaXX7r86QKsE0b/wFUNmWKVFqXgDAM1S7+LSrHz3Ul8GPdBb4NkABMyxPH53uh/WhjeKNDIJXNunsQ1Zp5ZFFxCO04bBzXjl6yCeYCaKXCZ7+229CAwBS2NDUsj7dmbFZ4p51YqD0rF8MJt5zryrty8wfJT5oWuGwTN1w+VjaRoaPd3WggLHOUqAMatCF4pv2CWcQwdGh/cvGrBbJV5XHYO3cZ/fK7CuYAydLy37isLoFtvq7z1oWtvQ+s+dFvc0nGq4Rc3SjFd+NRC0c6TzKcIS1zXdpSfJG27qcxlA9vSnrtDu35+Szx46tVhvMn2MG+D2bvDJOCT95Ojz8Ogd8gr8pq8IQn5SI7IV3JMpoSTn+QP+Uv+Rd8iEy2i5VoajYaafbIR0e9rdovhyg==</latexit>
sha1_base64="93VFV+aD4nbcJv8zKsh00rWrEoU=">AAACx3icbVHLbtQwFPWEVxleLSzZWIwqsUBVgooAicVIbGBXBDOtNIkqx7lJTf2I7JthRlYW/AFb+AW+iL/BSWbRmfZKto7PPdc+vjevpXAYx/9G0a3bd+7e27s/fvDw0eMn+wdP5840lsOMG2nsWc4cSKFhhgIlnNUWmMolnOaXH7v86RKsE0Z/w3UNmWKVFqXgDAM1T7+KSrHz/Ul8FPdBr4NkAybT96SPk/OD0d+0MLxRoJFL5twiiWvMPLMouIR2nDYOasYvWQWLADVT4DLf223pYWAKWhoblkbas1crPFPOrVUelIrhhdvNdeRNuUWD5bvMC103CJoPD5WNpGho93daCAsc5ToAxq0IXim/YJZxDB0aH169aslslXlcdQ7d1n98rsK5gDJ0vLfuKwugW2+rvPWha69C6950W9zScarhBzdKMV341ELRLpLMpwgrHGo7yk+Stt1W5rKBXWnP3aAdnt8Rbzz16jDeZHeY18H89VES8JfjyfTDMGeyR56TF+QlSchbMiWfyAmZEU6+k1/kN/kTfY5MtIxWgzQabWqeka2Ifv4HgefiQw==</latexit>

55Corr. shell 1(inequivalent shell 1)

Uncorrelated
orbitals

Uncorrelated
orbitals

32
5

32

32

⌃

Index of Wannier functions

Corr. shell 0(inequivalent shell 0)

Corr. shell 0(inequivalent shell 0)

Index of Wannier functions

0

0

0

[model]ncor = 2
norb = 3, 3
spin_orbit = False[model]ncor = 2
norb = 3, 3
spin_orbit = True

[model]ncor = 2
norb = 3
spin_orbit = True
corr_to_inequiv = 0, 0

Corr. shell 0(inequivalent shell 0)

Corr. shell 0(inequivalent shell 0)

Corr. shell 1(inequivalent shell 1)↑↑

Corr. shell 1(inequivalent shell 0)⌃

Uncorrelated orbitals

Uncorrelated orbitals<latexit sha1_base64="69B/Uk05O8z1FK2i5guc56ZfHuc=">AAACx3icbVHLbtQwFPWEVxleLV2ysRhVYoGqBIFgwaISG9i1gplWmkSV49ykpn5E9s0wIysL/oAt/AJf1L+pk8miM+2VbB2fe659fG9eS+Ewjq9G0b37Dx4+2nk8fvL02fMXu3svZ840lsOUG2nsWc4cSKFhigIlnNUWmMolnOaXX7r86QKsE0b/wFUNmWKVFqXgDAM1S7+LSrHz3Ul8GPdBb4NkABMyxPH53uh/WhjeKNDIJXNunsQ1Zp5ZFFxCO04bBzXjl6yCeYCaKXCZ7+229CAwBS2NDUsj7dmbFZ4p51YqD0rF8MJt5zryrty8wfJT5oWuGwTN1w+VjaRoaPd3WggLHOUqAMatCF4pv2CWcQwdGh/cvGrBbJV5XHYO3cZ/fK7CuYAydLy37isLoFtvq7z1oWtvQ+s+dFvc0nGq4Rc3SjFd+NRC0c6TzKcIS1zXdpSfJG27qcxlA9vSnrtDu35+Szx46tVhvMn2MG+D2bvDJOCT95Ojz8Ogd8gr8pq8IQn5SI7IV3JMpoSTn+QP+Uv+Rd8iEy2i5VoajYaafbIR0e9rdovhyg==</latexit>
sha1_base64="93VFV+aD4nbcJv8zKsh00rWrEoU=">AAACx3icbVHLbtQwFPWEVxleLSzZWIwqsUBVgooAicVIbGBXBDOtNIkqx7lJTf2I7JthRlYW/AFb+AW+iL/BSWbRmfZKto7PPdc+vjevpXAYx/9G0a3bd+7e27s/fvDw0eMn+wdP5840lsOMG2nsWc4cSKFhhgIlnNUWmMolnOaXH7v86RKsE0Z/w3UNmWKVFqXgDAM1T7+KSrHz/Ul8FPdBr4NkAybT96SPk/OD0d+0MLxRoJFL5twiiWvMPLMouIR2nDYOasYvWQWLADVT4DLf223pYWAKWhoblkbas1crPFPOrVUelIrhhdvNdeRNuUWD5bvMC103CJoPD5WNpGho93daCAsc5ToAxq0IXim/YJZxDB0aH169aslslXlcdQ7d1n98rsK5gDJ0vLfuKwugW2+rvPWha69C6950W9zScarhBzdKMV341ELRLpLMpwgrHGo7yk+Stt1W5rKBXWnP3aAdnt8Rbzz16jDeZHeY18H89VES8JfjyfTDMGeyR56TF+QlSchbMiWfyAmZEU6+k1/kN/kTfY5MtIxWgzQabWqeka2Ifv4HgefiQw==</latexit>

Corr. shell 0(inequivalent shell 0)

Corr. shell 0(inequivalent shell 0)
Corr. shell 1(inequivalent shell 1)↓

Uncorrelated orbitals

Corr. shell 1(inequivalent shell 0)
Uncorrelated orbitals↓

Figure 5: Shell structure of Wannier90 file with (a) spin_orbit=True and (b)spin_orbit=False. As an example, we consider cases with atoms in a unit cell(ncor=2). For the right panels of (a) and (b), the two atoms are treated equivalently (i.e., the self-energy is assumed to be identical).
parameters at each inequivalent shell are specified by the parameter kanamori asinteraction = kanamori
kanamori = [(U_1, U’_1, J_1), (U_2, U’_2, J_2), ... ].
(b) If interaction = slater_f, the interaction matrix is constructed by the effective
Slater integrals F0 , F2 , F4 , F6 [25]. These Slater integrals and the angular momentum lat each inequivalent shell are specified by the parameter slater_f asinteraction = slater_fslater_f = [(angular_momentum, F_0, F_2, F_4, F_6), ... ]It is noted that F0 , F2 , F4 , F6 must all be specified.
(c) If interaction = slater_uj, the Slater-type interaction is used. When U, J, andthe angular momentum l are set, the effective Slater integral is calculated internally in
DCore according to the formula defined in Table 2. The U, J, and l at each inequivalentshell are specified by the parameter slater_uj asinteraction = slater_uj
slater_uj = [(angular_momentum1, U1, J1), ... ]
Table

2:

Formulaoftheeffective

Slaterintegralsforinteraction = slater_uj [25, 26].
l
1
2
3

F0
U
U
U

F2
5J
14
1.0+0.63 J
6435
1001 J
286+195× 451
675 +250× 2025

F40.63F2451
675 F2F610012025 F2(iv) Local potential parameters
An arbitrary local potential can be implemented using parameters local_potential_
9SciPost Phys. 10, 117 (2021)

matrix and local_potential_factor. local_potential_matrix describes, in
the Python dictionary format, a set of inequivalent shell indices ish and a filename thatdefines the local potential matrix. The parameter local_potential_factor defines
a prefactor to the potential matrix. For details, see DCore’s online manual.
2. [system] block
This block includes thermodynamic parameters and some technical parameters such as
the inverse temperature and the number of Matsubara frequencies. It is also possible
to specify whether the chemical potential is fixed; if it is fixed, its value can be included. There is also the option with_dc, used to consider double-counting correction[see Eq. (15)]. Other parameters are described in the online document [27].
3. [impurity_solver] block
This block specifies the impurity solver to be used and the necessary parameters for running the solver program. In addition, the solver-dependent parameters can be specified
in this block. For details, see the online document [28].
4. [control] block
This block includes parameters that control the self-consistency loop of DMFT. For example, calculation condition parameters, such as the maximum steps of DMFT loops
and the convergence criteria, the input file names of the initial self-energies, and some
calculation options, such as the restart option, are specified as described in the online
document [29].
5. [tool] block
This block includes parameters for post calculations of real-frequency quantities, dcore_post. A real-frequency mesh and k-path are configured here. For more details, see theonline document [30].
6. [mpi] block
This block includes parameters that are read by dcore and dcore_post. The usage ofthe [mpi] block is explained in Sec. 6.2 and in the online document [31].

4.2

Output-file format

Output files are generated by each program. The output files for each program are described
below.
1. dcore_predcore_pre generates seedname.h5, where seedname is defined in a model block in
the input file. It has two groups, namely dft_input and DCore. See the DFTTools website [13] for details on the data structure in the dft_input group. In the DCore group,

the values of interaction matrix Uαβγδ (s) for each equivalent shell , where α, β, γ, δ deinote the spin-orbital indices at each inequivalent shell, and local potential Vs,o1,o2, wheres denotes the spin and o1, o2 denote orbitals, are output.
2. dcoredcore generates seedname.out.h5. All data are stored in the dmft_out group. In thisgroup, the input parameters, the total number of iteration steps, and the physical quantities, such as the local self-energy in the imaginary frequency domain and the chemical
potential, at each iteration step are output. The data structure of the self-energy is described in the online document [32]. For instance, the self-energy at the first iteration
for the 0-th inequivalent shell will be stored below /dmft_out/Sigma_iw/ite1/sh0
10

SciPost Phys. 10, 117 (2021)

in seedname.out.h5. An external impurity solver is invoked as a subprocess from themain process of dcore. For each iteration and each inequivalent shell, dcore generates
temporary input and output files in the working directory, e.g., work/imp_shell0_ite1
for the first iteration for the 0-th inequivalent shell. The user can check these temporary
files to diagnose problems.
3. dcore_checkdcore_check generates several pairs of a text file including numerical data and a figure file in PNG format in check directory. sigma.dat includes the local self-energyΣimp (iωn , s) at the final step. The corresponding figure file sigma_ave.png show theaverage self-energy of the last seven iterations, where the average is taken byP
ΣAve (iωn ) =s∈inequivalent shells

s
PNorb

α,β

Σαβ (iωn , s)ss∈inequivalent shells Norb

P

.

(16)

The maximum frequency is specified with the parameter omega_check in the [tool]
block. All other files have a prefix iter_, which indicates that the history during theiteration is stored. The text file iter_mu.dat contains the chemical potential as a function of the iteration number. iter_sigma-ish0.dat, iter_occup-ish0.dat, anditer_spin-ish0.dat stores the renormalization factor, the occupation number, andthe spin moment, respectively. Here, the number in -ish0 indicates the index for inequivalent shells. All data in the file iter_*.dat are plotted in iter_*.png.
4. dcore_post

dcore_post generates three text files, namely seedname_dos.dat, which includes thetotal spectral function, seedname_akw.dat, which includes the single-particle excitation spectrum A(k, w), and seedname_momdist.dat, which includes the momentum
distribution function. A script (seedname_akw.gp) for displaying A(k, w) is also generated at the same time. All files and figures are output in the post directory.

5

Installation

DCore is a collection of pure Python programs depending on TRIQS [11] and TRIQS/DFTTools[13]. DCore supports TRIQS 3.0.x and Python3. Once these prerequisites are installed, theinstallation of DCore is quick and efficient:$ pip3 install dcore
A pure-Python package will be automatically downloaded, being installed into the system
site-packages directory. If you do not have root privileges, you can install DCore into the user
site-packages directory:$ pip3 install dcore --user
You can find the locations of the installed files as$ pip3 show -f dcore
Please make sure that the PATH environment variable includes the directory containing theDCore executable files such as dcore, dcore_pre.
The detailed installation instructions can be found online [33]. In addition, at runtime,
one of the external impurity solvers shown below must be available.
11

SciPost Phys. 10, 117 (2021)• ALPS/CT-HYB [15]• ALPS/CT-HYB-SEGMENT [16]• TRIQS/cthyb [12]• TRIQS/Hubbard-I [21]• pomerol [22] (as Hubbard-I solver and ED solver)
We provide an up-to-date list of supported impurity solvers and short instructions on installation and usage online [34].

6

ExamplesThe examples below show how to define a model, how to compute a self-consistent solution,and how to compute physical quantities.

6.1

First example: Square-lattice Hubbard model

As the first example, we consider a single-band Hubbard model on a square lattice. Without
the shell index, the Hamiltonian in Eq. (3) is reduced toXX†H=
εk ckσckσ + Uni↑ ni↓ ,
(17)
i

kσ

†
where niσ = ciσciσ is the local number operator. The energy dispersion εk is given by

εk = 2t(cos k x + cos k y ) + 4t 0 sin k x sin k y .

(18)

Note that the nearest-neighbor hopping parameter t should be negative to minimize εk at theΓ point. This is an effective model of cuprate superconductors and is one of the most studied
models in strongly correlated electron systems.
We explain below the steps for obtaining the DMFT solution with the input file dmft_
square_pomerol.ini, which is available in the online tutorial [35].
6.1.1

Model construction

We first set up the model using dcore_pre. Parameters for defining a model are provided inthe [model] block of the input file. The square-lattice Hubbard model can be constructed asshown below.

[model]seedname = squarelattice = square
norb = 1
nelec = 1.0
t = -1.0
kanamori = [(4.0, 0.0, 0.0)]
nk = 8The square lattice is predefined in DCore and can be invoked simply by setting lattice =square. For more complicated models that are not predefined in DCore, this parameter isreplaced with lattice = wannier90 and the lattice data are described in a separate filewith the Wannier90 format. Sec. 6.1.5 describes its details. We consider half-filling, which is
12

SciPost Phys. 10, 117 (2021)

specified by nelec = 1.0. The value of transfer integral t is specified by t. In this case, t is
set as −1, and the absolute value of t is taken as the unit of energy. The interaction parameter
U is input by the first component of kanamori. The second and third components correspondto U 0 and J, respectively, but are meaningless for single-band models. The system size isdetermined by nk = 8, which means that the k-average in Eq. (6) is evaluated by summingup 8 × 8 k-points.
With this input file, dcore_pre is executed as$ dcore_pre dmft_square_pomerol.ini
This program generates an HDF5 file named square.h5, which includes all the informationon the model Hamiltonian, such as H(k).
6.1.2

Self-consistent calculations

We now proceed to the main program dcore. For executing the program, we need three additional blocks in the parameter file, namely [system], [impurity_solver], and [control].
The example below shows a parameter set for a square lattice:[system]
T = 0.1
n_iw = 1000
fix_mu = True
mu = 2.0[impurity_solver]
name = pomerol
exec_path{str} = pomerol2dcore
n_bath{int} = 3
fit_gtol{float} = 1e-6
[control]
max_step = 100
sigma_mix = 0.5
converge_tol = 1e-5
The [system] block includes T for temperature T and n_iw for the number of positive
Matsubara frequencies. One can use either the parameter β (inverse temperature) or T (temperature) to set the simulation temperature. We fix µ at 2 by specifying fix_mu = True andmu = 2.0 because the condition for half-filling is known. If fix_mu is not activated, µ is
adjusted every time the impurity problem is solved to make the occupation number equal tonelec.
The [impurity_solver] block specifies which impurity solver is used and its configuration. In this example, we use the ED solver pomerol [22], which is specified by name =pomerol. Other parameters in this block are solver-dependent and require the type specification, such as {str} and {int}. exec_path{str} is a full/absolute path to the executable(including its name). Here, we assume that the executable pomerol2dcore is in the PATH.
αβIn ED-based solvers, the hybridization function ∆σ (iωn ) in Eq. (13) is approximated by thefunction composed of Nbath poles:˜ αβ (iωn ) =
∆
σ

βb

N
bath
X

Vσαb (Vσ )∗b=1iωn − Eσb.

(19)

Hereafter, we drop the shell index s for simplicity. The energy level Eσb and the hybridizationαβ˜ αβparameter V αb are determined by fitting ∆σ (iωn ) to ∆σ (iωn ) [1]. The fitting tolerance isσ13shell=0, spin=up, 0
shell=0, spin=down, 0

0.66
0.64
0.62

shell=0, spin=up, 0
shell=0, spin=down, 0

10 2

diff

Renormalization factor

SciPost Phys. 10, 117 (2021)

10 3
10 4

0.60
0.58

10 5
2

4

6 8 10
iterations

12

2

4

6 8 10
iterations

12

Figure 6: One of the figures generated by dcore_check for a convergence check.
The renormalization factor z is plotted as a function of the iteration number (Leftpanel) and its absolute difference between iterations i and i − 1 is plotted in a logarithmic scale (Right panel).
specified by fit_gtol. The resultant finite-size system consisting of the correlated site plus
Nbath bath sites is solved using numerical diagonalization.
The [control] block includes parameters that control the self-consistent calculations.
The impurity problem is solved max_step = 100 times at maximum, but the loop is terminated if a convergence criterion is satisfied. The criterion we use is|O[i] − O[i − 1]| < converge_tol,(20)

where O[i] denotes a quantity at the i-th iteration. The quantities O include the chemical
potential µ and the renormalization factor z [see Eq. (21)]. sigma_mix specifies the mixing
parameter σmix in Eq. (14). The optimal value of sigma_mix depends on models and parameters. Generally, a larger value of sigma_mix (1 at the largest) can lead to a faster convergence,but may result in a bad convergence (e.g., oscillation and discontinuous change). A recommended strategy is that one starts from a large value of sigma_mix such as 1 and 0.5 anddecreases it if the convergence graph (see below) does not show a tendency of convergence.
The main program dcore is executed as$ dcore --np 4 dmft_square_pomerol.ini
Here, the number of MPI processes that are internally invoked is given after the –np option.
Note that dcore itself is launched as a single process. The results are stored in an HDF5 file
named square.out.h5.
After dcore is executed, one must check the convergence of the self-consistent calculations. For this purpose, the auxiliary program dcore_check is provided. dcore_check iscalled without the –np option as$ dcore_check dmft_square_pomerol.ini
Figure 6 shows one of the figures generated by dcore_check. The renormalization factor z
defined byzσ = [1 − ImΣσ (iω0 )/ω0 ]−1 ,
14

(21)

SciPost Phys. 10, 117 (2021)

1

CT-HYB

0.8
Z

0.6
0.4
0.2
0

0

1

2

3
Nbath

4

5

6

Figure 7: Renormalization factor z as a function of the number of bath sites Nbath
computed by the pomerol solver. The blue line shows the CT-QMC result z = 0.58.
is plotted as a function of the iteration number (left panel). z is bounded in the range 0 < z ≤ 1.
The right panel shows the difference of the renormalization factor between the i-th and (i −1)th iterations. The graph shows an exponential decay. The loop is terminated at the 13th
iteration, where the difference falls below converge_tol = 1e-5.
If the loop does not converge within max_step, one can continue the self-consistent calculations. To this end, one activates the restart option by adding a line in the input file:[control]restart = True
Then, dcore uses the final result for Σ(iωn ) in the previous run as the initial state and executes
additional max_step iterations. We note that the default is restart=False, which means
the self-consistent loop is started over even though the previous result exists in the current
directory. Other parameters in [control] block can be changed in the second run. For example,
sigma_mix may be reduced if the convergence graph does not show a converging behavior.
This procedure is repeated until convergence is reached.
At this point, the DMFT solutions have been obtained for static quantities and G(k, iωn )in the Matsubara domain. For example, the magnetization and the renormalization factor z
can be discussed. Figure 7 shows the convergence of z versus the number of bath sites Nbath .
Convergence to the CT-QMC result (explained in Sec. 6.1.4) is observed for Nbath ≥ 3.
6.1.3

Dynamical quantities

After the self-consistent calculations are finished, the post calculation dcore_post is conducted to compute dynamical quantities on the real frequency axis such as the density ofstates A(ω) and the single-particle excitation spectrum A(k, ω). The parameters for this step
are provided in the [tool] block as follows:[tool]
knode = [(G,0,0,0),(X,0.5,0,0),(M,0.5,0.5,0),(G,0,0,0)]nk_line = 100
omega_max = 6.0
omega_min = -6.0
Nomega = 401
broadening = 0.4

15

SciPost Phys. 10, 117 (2021)

The k-path is specified by the parameter knode. In this example, k starts from Γ and comesback to Γ through X and M. For each interval, the path is divided into nk_line = 100 points,
on which A(k, ω) is computed. The ω mesh is generated using the parameters omega_min,omega_max, and Nomega. The parameter broadening specifies the extent of an artificial
broadening δ in the spectrum. This is done by replacing ω by ω + iδ in the analytical continuation. broadening is set to a value of the order of T to obtain a smooth spectrum with the
ED solver. For other solvers that treat the thermodynamic limit, one should set broadening= 0 to avoid artificial broadening.
dcore_post is executed with the –np option because MPI parallel computation is used:$ dcore_post --np 4 dmft_square_pomerol.ini
After the computation is finished, numerical results are stored in post directory. The spectrum
A(k, ω) on a given k-path, for example, is saved in a text file square_akw.dat together witha plotting script square_akw.gp in gnuplot format. One can plot the data by$ gnuplot square_akw.gp
Figure 8(a) shows A(k, ω) plotted with this command. We note that the range of the colorbar
needs to be changed manually to obtain a distinct spectrum. Although the ED solver yielded
reasonable estimations for z, as shown in Fig. 7, the dynamical quantities clearly show artificial
features. Because ∆(ω) is composed only of n_bath = 3 poles located at ω = 0 and ±1.73,
A(k, ω) exhibits discrete features. However, the hybridization should take place equally in thewhole energy region. A larger n_bath is necessary to discuss dynamical properties.
We note that the Padé approximation used for the analytical continuation does not ensure
the causality, namely, A(k, ω) may becomes negative at some (k, ω). One can partly avoid
this unphysical behavior by increasing the value of broadening. Another risk is involved in
the analytical continuation from noisy Σ(iωn ), which will be described in 6.1.4.
6.1.4

CT-QMC solver

The CT-QMC method offers an unbiased simulation of general interacting models [36, 37]. In
particular, efficient calculations for the effective impurity problem are achieved by its hybridization-expansion formulation (CT-HYB), which performs Monte Carlo sampling in the expansion
with respect to hybridization between the correlated atom and the bath as perturbation [38].
In special cases without exchange interactions, an even more efficient algorithm based on

Figure 8: A(k, ω) for each spin component computed by (a) the ED solver with threebath sites and (b) the CT-QMC solver.

16

SciPost Phys. 10, 117 (2021)

a segment picture can be used [39]. Here, we use ALPS/CT-HYB-SEGMENT, which is the
ALPSCore [40] version of the implementation developed by Hafermann et al. [41]. A detailed
documentation with overview of parameters and sample scripts is provided in the subdirectory
Documentation in the repository [16].
The segment algorithm can handle only the density-density type interactions. This is thecase with the single-orbital Hubbard model and the multi-orbital Hubbard model with J = 0.
When the segment solver is attempted to use, DCore check if the model includes non-densitydensity interactions and if yes, the program stops before the impurity solver is invoked. Even insuch cases, one can still apply the segment solver by activating the density_density option
as[model]density_density = True
Then, the interactions generated by dcore_pre are restricted to density-density type. Notethat the density_density option changes the model (e.g. when J is finite) and hence theresults may differ from those obtained by other solvers such as ALPS/CT-HYB.
To invoke the CT-HYB solver, the parameters in the [impurity_solver] block of the
input file are changed as follows:[impurity_solver]
name = ALPS/cthyb-seg
exec_path{str} = /path/to/alps_cthyb
cthyb.TEXT_OUTPUT{int} = 1
MAX_TIME{int} = 60
cthyb.N_MEAS{int} = 50
cthyb.THERMALIZATION{int} = 100000
cthyb.SWEEPS{int} = 100000000
The CT-QMC solver is invoked by name = ALPS/cthyb-seg. The second parameter exec_
path{str} should be changed to point to the path to the alps_cthyb executable in theenvironment. Important parameters for the QMC simulation are• MAX_TIME{int},• cthyb.N_MEAS{int},• cthyb.THERMALIZATION{int},• cthyb.SWEEPS{int}.
The simulation finishes either when cthyb.SWEEPS{int} Monte Carlo sweeps are finished orwhen the elapsed time reaches MAX_TIME seconds. In this example, cthyb.SWEEPS{int} isso large that the computation time is controlled by MAX_TIME{int}. cthyb.N_MEAS{int}is the number of updates per measurement, and here we use 50 following the example in theofficial documentation [16]. cthyb.THERMALIZATION{int} specifies the number of updates
before the measurement starts. cthyb.THERMALIZATION{int} = 100000 should be large
enough to properly discard irrelevant samples. However, this value may not be large enoughat lower temperatures. This should be checked by, for example, plotting the histogram of the
expansion order [41]. The whole list of optional parameters can be found by running thecommand$ alps_cthyb --help
Note that solver-dependent parameters should be input with the type specification, e.g., {int}and {str}, as in the above example.
17

0.66

shell=0, spin=up, 0
shell=0, spin=down, 0

10 2

0.64

shell=0, spin=up, 0
shell=0, spin=down, 0

10 3
diff

Renormalization factor

SciPost Phys. 10, 117 (2021)

0.62

10 4

0.60
0.58

10 5
2 4 6 8 10 12 14 16 18 20
iterations

2 4 6 8 10 12 14 16 18 20
iterations

Figure 9: Case of CT-HYB solver for a figure generated by dcore_check for a convergence check. See the caption in Fig. 6.
Because the QMC solver is time-consuming, it is important to use a proper initial guessfor Σ(iωn ) to reduce the number of iterations. For example, we can use the converged result
obtained for similar parameters or the result obtained for the same parameters with a different
solver. This can be done by assigning the path of sigma.dat generated by dcore_check to
the initial_self_energy parameter as[control]initial_self_energy = /path/to/sigma.dat
max_step = 20
sigma_mix = 0.5
time_reversal = True
Here, /path/to/sigma.dat should be changed according to the environment. The parameter time_reversal = True activates the spin average of Σσ (iωn ). This is particularly useful for the QMC solver to improve statistics. Unlike the previous example using the ED solver,
converge_tol is not specified here, because statistical errors in the QMC method make theautomatic convergence check less reliable. The convergence is checked visually using thegraphs generated by dcore_check.
With this input file, all programs, namely dcore_pre, dcore, dcore_check, and dcore_post, are executed as in Sec. 6.1.1. Figure 9 shows graphs generated by dcore_check. The
right panel shows that the statistical error is on the order of 10−3 and that convergence isreached at around the 10th iteration. The result for A(k, ω) is shown in Fig. 8(b). Note that
artificial broadening is turned off by broadening = 0.0. The spectrum exhibits low-energy
quasiparticle excitations and high-energy broadening due to correlations.
We note that spectra computed using the Padé approximation is extremely sensitive tostatistical errors. For this reason, Fig. 8(b) might not be reproduced even with the same input
as above. In such cases, one can try improving the QMC accuracy by increasing the number of
MPI processes or by increasing MAX_TIME{int} parameter. An alternative way for analyticalcontinuation is discussed in Sec. 7 as a future development.
Here, a comment on the convergence check is in order. Although the automatic convergence check has been turned off in the above example, there is a way to employ the convergence check even for QMC solvers. In the case with Fig. 9, the following configuration works:

18

SciPost Phys. 10, 117 (2021)[control]converge_tol = 0.002
n_converge = 5
The DMFT loop is terminated when the convergence criterion (20) is satisfied n_converge
times consecutively. converge_tol is set larger than the magnitude of statistical errors, andthe other parameter n_converge=5 prevents the loop from being terminated prematurely.
This way of automatic convergence check can be attempted to use if the expected magnitudeof statistical errors is known in advance.
As discussed above, one could alternatively use a CT-HYB solver based on the matrix algorithm, such as ALPS/CT-HYB [see Sec. 6.2]. To obtain the same results with comparable
accuracy, one may have to run a CT-HYB solver based on the matrix algorithm longer in runtime, typically by one order of magnitude.
6.1.5

Wannier90 interface

DCore provides simple lattice models such as a one-dimensional chain, a two-dimensionalsquare lattice, and a three-dimensional cubic lattice. Arbitrary lattices that are not predefinedin DCore can be implemented using the Wannier90 format [42]. In the following, we show
how to construct the square lattice using Wannier90 input.
The only difference in dmft_square_pomerol.ini is in the [model] block, which includes[model]seedname = squarelattice = wannier90
norb = 1
nelec = 1.0
kanamori = [(4.0, 0.0, 0.0)]
nk = 8The difference between this block and that in Sec. 6.1.1 is lattice = wannier90. When
this is given, DCore reads a wannier90 file named seedname_hr.dat (square_hr.dat in
the present case) to set up the one-body part of the Hamiltonian, H(k). square_hr.dat
describing the square lattice is given bysquare lattice
1
9
1 1 1 1 1 1 1 1 1
0 0 0 1 1
0.0 0.00 1 0 1 1 -1.0 0.0
1 0 0 1 1 -1.0 0.0
0 -1 0 1 1 -1.0 0.0-1 0 0 1 1 -1.0 0.0
1 1 0 1 1
0.0 0.0
1 -1 0 1 1
0.0 0.0-1 1 0 1 1
0.0 0.0-1 -1 0 1 1
0.0 0.0
The first line can be any text. The second line specifies the number of orbitals in a unit cell. The
next line specifies the number of cells, Ncell , for which the transfer integrals are provided. Here,
we have Ncell = 9 cells, as shown in Fig. 10(a). The value of transfer integrals, H iW90
(R), in the
j
P
†
W90
0
one-body Hamiltonian H = RR0 i j H i j (R − R)cRi cR0 j is listed for all possible combinations
19

SciPost Phys. 10, 117 (2021)

(a)

(b)
1

1

1

1

2

1

2

1

2

1

1

1

3

4

3

4

3

4

1

1

1

1

2

1

2

1

2

3

4

3

4

3

4

1

2

1

2

1

2

3

4

3

4

3

4

Figure 10: Lattice configuration for the Wannier90 input for (a) the paramagnetic
state and (b) the AFM state. The circles show sites (atoms) and the number is thelabel of the site in the unit cell indicated by the square. The arrows indicate thehopping included in the input file.
of R, i, and j for fixed R0 = 0. Here, i and j denote the combined orbital and site indices in
the unit cell. The format is “l1 l2 l3 i j Re H iW90
(R) Im H iW90
(R)”, where R ≡ l1 a1 + l2 a2 + l3 a3jjand (a1 , a2 , a3 ) are the unit lattice vectors. The convention used by Wannier90 is related tothe one in Eq. (4) as H iW90(−R) = H i j (R).
j
6.1.6

Antiferromagnetic stateWith the Wannier90 format, one can construct arbitrary lattices. For example, we can enlarge
the unit cell to allow an antiferromagnetic (AFM) solution. Let us see how this is done for thesquare lattice example.
Figure 10(b) shows the lattice configuration for AFM calculations. Four “orbitals” are contained in a unit cell. The Wannier90 format file now includes the hopping inside the unit cell
(8 double-pointed arrows) and the hopping from the central cell to neighboring unit cells (8single-pointed arrows). The full input files as well as a script for generating the Wannier90
format file are available in the online tutorial [43].
The parameters in the [model] block of the input file are given as follow:[model]seedname = afm_dim2
lattice = wannier90
ncor = 4
norb = 1, 1
nelec = 4.0
kanamori = [(4.0, 0.0, 0.0), (4.0, 0.0, 0.0)]
corr_to_inequiv = 0, 1, 1, 0
nk = 8The total number of orbitals contained in the unit cell is specified by ncor = 4. nelec isthe total number of electrons in the unit cell, and hence needs to be multiplied by 4. In the
AFM state, we expect a staggered ordered state, where site 1 and site 4 are equivalent (e.g.,
spin-up state), and site 2 and site 3 are equivalent (spin-down state). This constraint can be
imposed by setting corr_to_inequiv = 0, 1, 1, 0, which indicates that there are only
20

SciPost Phys. 10, 117 (2021)

10 2

0.20
shell=0, Sx
shell=0, Sy
shell=0, Sz

0.15
0.10

diff

Spin moment

0.25

shell=0, Sx
shell=0, Sy
shell=0, Sz

10 4

0.05
0.00

10 3

4

8 12 16 20 24 28
iterations4

8 12 16 20 24 28
iterations

Figure 11: Spin moment mξ (s) (ξ = x, y, z) for inequivalent site s = 0 as a functionof the iteration number. The left panel shows mξ (s) itself and the right panel shows
the difference between iterations i and i − 1.
two inequivalent sites: site 1 and site 4 are labeled with 0, and site 2 and site 3 are labeledwith 1. Two independent impurity problems are solved in each step of the self-consistentcalculations. For each impurity model, the number of orbitals and the strength of interactions
are specified by the norb and kanamori parameters, respectively.
To obtain an AFM solution, we need to start the self-consistent calculation with a brokensymmetry initial guess. To this end, we set the initial self-energy toσΣσ
αβ (iωn , s) = vαβ (s),(22)σ
using local potential vαβ(s), which depends on the inequivalent site s. The [control] blockcontains an additional parameter that specifies the initial self-energy as[control]
max_step = 30
sigma_mix = 0.5
initial_static_self_energy = {0: ’init_se_up.txt’, 1: ’init_se_down.txt’}
A path to a text file is given in initial_static_self_energy for each inequivalent siteσusing the dictionary literal in Python. The text files describe non-zero components in vαβinσσ
the form “σ α β Re vαβ Im vαβ ”. To obtain a magnetic state (spin-up state), we can input themagnetic field in init_se_up.txt as

0 0 0 -1.0 0.0
1 0 0 1.0 0.0
and the magnetic field with the opposite direction is given in init_se_down.txt as

0 0 0 1.0 0.0
1 0 0 -1.0 0.0
The parameters in [system] and [impurity_solver] do not change from the paramagnetic
calculations.
Figure 11 shows the convergence of the spin moment mξ (s) (ξ = x, y, z) for inequivalentsite s = 0. A convergence of mz to a non-zero value was obtained around the 30th iteration.
The moment mz = 0.283 is about half the full moment (1/2). We confirmed that mz (s) for
21

SciPost Phys. 10, 117 (2021)

s = 1 converges to the negative side, namely mz (0) = −mz (1) (no figure). The parameters
U = 4, n = 1, and T = 0.1 are the same as those in the previous calculations in Fig. 8,
meaning that the paramagnetic solution was not a true solution of the DMFT equation. The
phase diagram is shown, for example, in Ref. [44].

6.2

Solving multi-orbital model using QMC solver: t 2g model on a Bethe lattice

In this subsection, we show how to solve a three-orbital model on a Bethe lattice by using
DCore. In this model, an interesting phenomenon called spin-freezing transition occurs [45].
Spin-freezing is signaled by a peculiar frequency dependence of the self-energy: ImΣ(iωn ) ∝ω0.5
n . We will reproduce this result using DCore.
6.2.1

Model definition

We first make the input file of dcore_pre for generating an HDF5 file that is necessary for
DMFT calculations. The Hamiltonian in Ref. [45] is defined asXX†
H=
εαk ckασ
ckασ +
Hloc (R),
(23)
R

k,σ

Hloc (R) = −−XµnRασ +

X

α,σ

X
α6=β

X

U nRα↑ nRα↓ +αU 0 nRασ nRβ−σ + (U 0 − J)nRασ nRβσ

α>β,σ

†
†
†
†
J(cRα↓cRβ↑
cRβ↓ cRα↑ + cRβ↑cRβ↓cRα↑ cRα↓ + H.c.),(24)

where α = 1, 2, 3 is the orbital index. We use a Kanamori interaction with U = 8, U 0 = U−2J =5.3333333, and J = 1.33333. We use the same model parameters as those used in the paper [45]: n = 1.6, t = 1.0, U/t = 8.0, J/U = 1/6. The input file for dcore_pre is as follows:[model]lattice = bethe
seedname = bethe
nelec = 1.6
t=1.0
norb = 3
kanamori = [(8.0, 5.3333333, 1.33333)]
nk = 1000
After running dcore_pre, an HDF5 file named bethe.h5 is generated. DCore discretizes
the density of states of the Bethe lattice, a semicircular density of states on [−2t, 2t], along a
virtual one-dimensional k axis with nk = 1000 points.
6.2.2

DMFT calculation

Next, we make the input file of dcore. The calculation in Ref. [45] was done using the matrix
formalism of the CT-HYB method. To make a direct comparison, we use an implementation of
the same algorithm, ALPS/CT-HYB, which was developed by one of the authors. An example
input file of dcore is shown below.

[mpi]
command = ’$MPIRUN -np #’[system]beta = 50.0
[impurity_solver]
22

SciPost Phys. 10, 117 (2021)

name = ALPS/cthyb
timelimit{int} = 300
exec_path{str} = hybmat[control]
max_step = 40
sigma_mix = 1.0
restart = False
In this example file, we omit the model block, which is the same as that in the input file ofdcore_pre. The inverse temperature β is set to 50t. We also define a command for MPI
parallelization using the mpi block. Here, # in command of the mpi block is replaced by the
number of processes specified at runtime. In the above example, we define the environment
variable MPIRUN and run the program with 24 MPI processes as follows.

$ export MPIRUN="mpirun"$ dcore dmft_bethe.ini --np 24
After running dcore, the results of the self-energy and Green’s functions in each iteration are
accumulated in an HDF5 file named *seedname*.out.h5 (bethe.out.h5 in the presentcase) and a text file that contains the self-energies (sigma.dat) is output in the check directory.
Before analyzing the computed results such as the self-energy, the user may check theconvergence of CT-QMC sampling. The impurity solver runs for 300 seconds at each iteration
step, and the standard output and error of the solver are redirected to the standard output
of dcore. The last part of the output of the solver at the last iteration looks as follows. The
perturbation orders measured just before and after the measurement steps are close to eachother (' 58), indicating that the thermalization time was long enough: ALPS/CT-HYB used 10% of the total simulation time (timelimit) for thermalization. One can see that the averagesign is close to 1, indicating that there is no severe sign problem.

==== Thermalization analysis ====
Perturbation orders just before and after measurement steps are 57.6833
and 57.9479.
==== Number of Monte Carlo steps spent in configuration spaces ====
Z function: 5.03099
G1: 4.97089==== Acceptance updates of operators hybridized with bath ====
1-pair_insertion_remover: 0.177554
2-pair_insertion_remover: 0.0142077
Single_operator_shift_updater: 0.196652
Operator_pair_flavor_updater: 0.166643==== Acceptance rates of worm updates ====
G1_ins_rem: 0.0339203
G1_mover: 0.214957
G1_flavor_changer: 0.0104972
G1_ins_rem_hyb: 0.139718
Average sign is 0.999800707488.
Total charge of impurity problem: 1.611984
23

SciPost Phys. 10, 117 (2021)

After checking the convergence of the self-consistent iterations using dcore_check, one
can plot the self-energy stored in sigma.dat as shown in Fig. 12. The solid line shows theresults taken from [45]. In the low-frequency region, ImΣ(iωn ) is proportional to ω1/2
n . Thethree symbols represent the orbital diagonal components of ImΣ(iωn ) computed by DCore.
The results obtained by dcore match those of [45].

3
dxy
dyz
dzx
PRL 101, 166405 (2008)−ImΣ(iωn )

2.5
2

1.5
1

0.5
0
0

0.5

1

1.5

2

ωn1/2
Figure 12: Matsubara frequency dependence of ImΣ(iωn ). The solid line indicatesthe numerical results reported in Ref. [45]. The three symbols indicate ImΣ(iωn ) ofd x y , d yz , and dz x orbitals, respectively, obtained by DCore.

6.3

DFT+DMFT example: SrVO3In this section, we demonstrate how to perform a multi-orbital DFT+DMFT calculation for a
realistic band structure. In particular, we use SrVO3 which has been studied by DMFT as atestbed in many previous studies.
6.3.1

Construction of Wannier functions

We fist construct maximally localized Wannier functions for the t 2g manifold using Quantum
ESPRESSO [46] and Wannier90 [42]. Any DFT programs that support Wannier90 may be usedalternatively. For our DFT calculation, we take a lattice constant of 7.29738 a.u and constructa three-orbital tight-binding model following a common procedure using the local density
approximation (LDA). We provide a more detailed description on how to construct Wannierfunctions using Quantum ESPRESSO or another DFT program, OpenMX [47] online [48].
There, we can download input files as well as the data file of the resultant tight-binding model.

6.3.2

Model definition

We show the model block of our input file for performing a DMFT calculation using theWannier functions below. We adopted the Kanamori interaction parameters (U = 3.44 eV,
U 0 = 2.49 eV, J = 0.46 eV) estimated in Ref. [49]. The number of grid points in the first
Brillouin zone along each reciprocal vector is set to 10. This number does not have to matchthe one used in the DFT calculation. The total number of k points is 103 .
24

Energy from EF [eV]SciPost Phys. 10, 117 (2021)

6
4
2
0
-2
-4
-6
-8

¡

X

M

¡

R

Figure 13: Comparison of LDA band structure and fit by the Wannier functions. The
circles denote the energies computed by LDA and the solid lines denote the Wannier
fit.

[model]lattice = wannier90
seedname = srvo3
nelec = 1.0
ncor = 1
norb = 3
kanamori = [(3.44, 2.49, 0.46)]
bvec=[(1.627091,0.0,0.0),(0.0,1.627091,0.0),(0.0,0.0,1.627091)]
nk = 10

6.3.3

Parameters for self-consistent calculations and post processing

We set input parameters for ALPS/CT-HYB in the impurity_solver block. The parameters
timelimit and exec_path are passed to ALPS/CT-HYB. The impurity solver runs for 100
seconds with 96 MPI processes at each self-consistent iteration. We set an initial guess for thechemical potential to an appropriate value and inverse temperature to β = 20 eV−1 ' 580 K.
In the input file shown below, we specify the temperature by using the parameter beta insteadof the parameter T for the sake of brevity in appearance.

[mpi]
command = ’$MPIRUN -np #’[system]beta = 20.0
mu = 12.290722[impurity_solver]
name = ALPS/cthyb
timelimit{int} = 400
exec_path{str} = hybmat[control]
max_step = 20
time_reversal = True
25

SciPost Phys. 10, 117 (2021)

sigma_mix = 0.8
The tool block includes input parameters for post processing. In this example, we compute A(k, ω) using dcore_post along the k-paths specified by symmetry points given in
knode. As set by the parameter omega_pade, the data of the self-energy in the energy windowof 0 < ωn ≤ 2 is used for analytic continuation to the real-frequency axis. Here, we introduce
a small broadening factor for computing spectral functions.

[tool]
broadening = 0.01
nk_line = 50
knode=[(G,0,0,0),(X,0.5,0,0),(M,0.5,0.5,0),(G,0,0,0),(R,0.5,0.5,0.5)]
omega_max = 2.0
omega_min = -2.0
Nomega = 400
omega_check = 30.0
omega_pade = 2.0

6.3.4

Results

In Fig. 13, one can see that the t 2g bands are fitted very well by the Wannier functions. Thet 2g manifold is separated in energy space from the eg bands above 1.5 eV and oxygen bands
below −2 eV. Since we constructed the t 2g model, the eg bands and oxygen bands are not taken
into account in the present DMFT calculation. Figure 14 shows A(k, ω) computed by DMFT
Compared to the DFT band structure in Fig. 13, one can see that the t 2g band is substantiallyrenormalized by strong correlation effects as expected.
There are several previous DMFT studies for this compounds [50–52]. These calculations
aimed at a quantitative description of the spectrum of the real material and therefore took into
account other bands such as the eg bands and oxygen bands. Thus, a direct comparison with
our result is not possible.
For DCore, one could take into account those uncorrelated bands by constructing Wannierfunctions for the full 3d orbitals and oxygen p orbitals in a wider energy window. In this case,
we may have to carefully adjust (manually) the level splitting between t 2g and the uncorrelated
orbitals because the current version of DCore does not support intra-shell (site) interactions
between the d and p shells.
We need extra care for quantum Monte Carlo simulations since they often suffer froma sign problem and many other technical issues. If the reader faces such a problem, it isrecommended to create an issue on GitHub at the issues page [53] or directly contact thedevelopers of the solver.

7

Conclusion

We introduced the open-source software package DCore, which is integrated DMFT software
for correlated electrons. We described the algorithm, the structure of the package, installation,and file format, and demonstrated usage for the single-/multi-orbital Hubbard models and theWannier90 model for SrVO3 .
We are planning to implement more functionality in a future version. The current version
supports only the Padé approximation for analytic continuation. However, the Padé approximation is sensitive to statistical noise. A future version will support external analytic continuation solvers based on more stable methods such as SpM [54] and Maxent [55]. Furthermore,
we will implement recently proposed efficient methods for computing static/dynamical lattice26

SciPost Phys. 10, 117 (2021)

2

10
8

1
0.5

6

0
4

-0.5
-1

A(k, !) (eV-1)

Energy (eV)

1.5

2

-1.5
-2

¡

X

M

¡

R¡

X

M

¡

R

0

Figure 14: A(k, ω) for each spin component computed by DMFT for the Wanniermodel of SrVO3 .
susceptibilities based on the Bethe-Salpeter equation [44, 56]. There are also plans to support
more external impurity solvers such as iQIST [57] and w2dynamics [18].

Acknowledgements
We acknowledge Takeo Kato and Yuichi Motoyama for supporting the DCore project. We
acknowledge the TRIQS community and the ALPS community for developing useful opensource libraries.
Author contributions For an updated record of individual contributions, consult the repository at https://issp-center-dev.github.io/DCore/master/index.html.
Funding information H.S., J.O., and K.Y. were supported by JSPS KAKENHI Grant No.
18H01158, 21H01003, 21H01041. H.S. was supported by JSPS KAKENHI Grant No. 16K17735.
J.O. was supported by JSPS KAKENHI Grant No. 18H04301 (J-Physics). N.T. was supported by
JSPS KAKENHI Grant No. 19H05817, and 19H05820. K.Y. was supported by JSPS KAKENHI
Grant No. 19K03649, and Building of Consortia for the Development of Human Resources inScience and Technology, MEXT, Japan.

References[1] A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg, Dynamical mean-field theory ofstrongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys.
68, 13 (1996), doi:10.1103/revmodphys.68.13.
[2] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet and C. A. Marianetti,
Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78,
865 (2006), doi:10.1103/revmodphys.78.865.

27

SciPost Phys. 10, 117 (2021)[3] D. Vollhardt, K. Held, G. Keller, R. Bulla, Th. Pruschke, I. A. Nekrasov and V. I. Anisimov,
Dynamical mean-field theory and its applications to real materials, J. Phys. Soc. Jpn. 74,
136 (2005), doi:10.1143/jpsj.74.136.
[4] C. Weber, K. Haule and G. Kotliar, Strength of correlations in electron- and hole-doped
cuprates, Nat. Phys. 6, 574 (2010), doi:10.1038/nphys1706.
[5] M. Aichhorn, L. Pourovskii and A. Georges, Importance of electronic correlations for structural and magnetic properties of the iron pnictide superconductor LaFeAsO, Phys. Rev. B 84,
054529 (2011), doi:10.1103/physrevb.84.054529.
[6] E. Razzoli et al., Tuning electronic correlations in transition metal pnictides: Chemistry beyond the valence count, Phys. Rev. B 91, 214502 (2015),
doi:10.1103/physrevb.91.214502.
[7] Z. P. Yin, K. Haule and G. Kotliar, Kinetic frustration and the nature of the magnetic andparamagnetic states in iron pnictides and iron chalcogenides, Nat. Mater. 10, 932 (2011),
doi:10.1038/nmat3120.
[8] J. H. Shim, K. Haule and G. Kotliar, Modeling the localized-to-itinerant electronic transition in the heavy fermion system CeIrIn5 , Science 318, 1615 (2007),
doi:10.1126/science.1149064.
[9] S. Y. Savrasov, G. Kotliar and E. Abrahams, Correlated electrons in δ-plutonium within adynamical mean-field picture, Nature 410, 793 (2001), doi:10.1038/35071035.
[10] K. Held, A. K. McMahan and R. T. Scalettar, Cerium volume collapse: Results from themerger of dynamical mean-field theory and local density approximation, Phys. Rev. Lett.
87, 276404 (2001), doi:10.1103/physrevlett.87.276404.
[11] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio and P. Seth, TRIQS:
A toolbox for research on interacting quantum systems, Comput. Phys. Commun. 196, 398(2015), doi:10.1016/j.cpc.2015.04.023.
[12] P. Seth, I. Krivenko, M. Ferrero and O. Parcollet, TRIQS/CTHYB: A continuous-time quantum Monte Carlo hybridisation expansion solver for quantum impurity problems, Comput.
Phys. Commun. 200, 274 (2016), doi:10.1016/j.cpc.2015.10.023.
[13] M. Aichhorn et al., TRIQS/DFTTools: A TRIQS application for ab initio calculations of correlated materials, Comput. Phys. Commun. 204, 200 (2016),
doi:10.1016/j.cpc.2016.03.014.
[14] B. Bauer et al., The ALPS project release 2.0: open source software for strongly correlated
systems, J. Stat. Mech. P05001 (2011), doi:10.1088/1742-5468/2011/05/p05001.
[15] H. Shinaoka, E. Gull and P. Werner, Continuous-time hybridization expansion quantum
impurity solver for multi-orbital systems with complex hybridizations, Comput. Phys. Commun. 215, 128 (2017), doi:10.1016/j.cpc.2017.01.003.
[16] ALPSCore/CT-HYB-SEGMENT, https://github.com/ALPSCore/CT-HYB-SEGMENT.
[17] H. Shinaoka, Y. Nomura and E. Gull, Efficient implementation of the continuous-time
interaction-expansion quantum Monte Carlo method, Comput. Phys. Commun. 252,
106826 (2020), doi:10.1016/j.cpc.2019.06.016.

28

SciPost Phys. 10, 117 (2021)[18] M. Wallerberger, A. Hausoel, P. Gunacker, A. Kowalski, N. Parragh, F. Goth, K. Held and
G. Sangiovanni, w2dynamics: Local one- and two-particle quantities from dynamical mean
field theory, Comput. Phys. Commun. 235, 388 (2019), doi:10.1016/j.cpc.2018.09.007.
[19] V. Singh, U. Herath, B. Wah, X. Liao, A. H. Romero and H. Park, DMFTwDFT:
An open-source code combining Dynamical Mean Field Theory with various density functional theory packages, Comput. Phys. Commun. 261, 107778 (2021),
doi:10.1016/j.cpc.2020.107778.
[20] K. Haule, DFT + embedded DMFT Functional, Rutgers (2007), http://hauleweb.rutgers.
edu/tutorials/index.html.
[21] hubbardI, https://triqs.github.io/hubbardI.
[22] A. E. Antipov, I. Krivenko and S. Iskakov, pomerol, Zenodo, v1.2 (2017),
doi:10.5281/zenodo.825870.
[23] M. Kawamura, K. Yoshimi, T. Misawa, Y. Yamaji, S. Todo and N. Kawashima,
Quantum lattice model solver HΦ, Comput. Phys. Commun. 217, 180 (2017),
doi:10.1016/j.cpc.2017.04.006.
[24] V. I. Anisimov, F. Aryasetiawan and A. I. Lichtenstein, First-principles calculations of theelectronic structure and spectra of strongly correlated systems: the LDA+U method, J. Phys.:
Condens. Matter 9, 767 (1997), doi:10.1088/0953-8984/9/4/002.
[25] V. I. Anisimov, F. Aryasetiawan and A. I. Lichtenstein, First-principles calculations of theelectronic structure and spectra of strongly correlated systems: the LDA+U method, J. Phys.:
Condens. Matter 9, 767 (1997), doi:10.1088/0953-8984/9/4/002.
[26] V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk and G. A. Sawatzky, Densityfunctional theory and NiO photoemission spectra, Phys. Rev. B 48, 16929 (1993),
doi:10.1103/physrevb.48.16929.
[27] https://issp-center-dev.github.io/DCore/v3.0.0/reference/input.html#system-block.
[28] https://issp-center-dev.github.io/DCore/v3.0.0/reference/input.html#
impurity-solver-block.
[29] https://issp-center-dev.github.io/DCore/v3.0.0/reference/input.html#control-block.
[30] https://issp-center-dev.github.io/DCore/v3.0.0/reference/input.html#tool-block.
[31] https://issp-center-dev.github.io/DCore/v3.0.0/reference/input.html#mpi-block.
[32] https://issp-center-dev.github.io/DCore/v3.0.0/reference/gf_format.html.
[33] https://issp-center-dev.github.io/DCore/v3.0.0/install.html.
[34] https://issp-center-dev.github.io/DCore/v3.0.0/impuritysolvers.html.
[35] https://issp-center-dev.github.io/DCore/v3.0.0/tutorial/square/square.html.
[36] A. N. Rubtsov, V. V. Savkin and A. I. Lichtenstein, Continuous-time quantum Monte Carlo
method for fermions, Phys. Rev. B 72, 035122 (2005), doi:10.1103/physrevb.72.035122.
[37] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer and P. Werner, Continuoustime Monte Carlo methods for quantum impurity models, Rev. Mod. Phys. 83, 349 (2011),
doi:10.1103/revmodphys.83.349.
29

SciPost Phys. 10, 117 (2021)[38] P. Werner and A. J. Millis, Hybridization expansion impurity solver: General formulation
and application to Kondo lattice and two-orbital models, Phys. Rev. B 74, 155107 (2006),
doi:10.1103/physrevb.74.155107.
[39] P. Werner, A. Comanac, L. de’ Medici, M. Troyer and A. J. Millis, Continuoustime solver for quantum impurity models, Phys. Rev. Lett. 97, 076405 (2006),
doi:10.1103/physrevlett.97.076405.
[40] A. Gaenko et al., Updated core libraries of the ALPS project, Comput. Phys. Commun. 213,
235 (2017), doi:10.1016/j.cpc.2016.12.009.
[41] H. Hafermann, P. Werner and E. Gull, Efficient implementation of the continuous-time
hybridization expansion quantum impurity solver, Comput. Phys. Commun. 184, 1280
(2013), doi:10.1016/j.cpc.2012.12.013.
[42] G. Pizzi et al., Wannier90 as a community code: new features and applications, J. Phys.:
Condens. Matter 32, 165902 (2020), doi:10.1088/1361-648x/ab51ff.
[43] https://issp-center-dev.github.io/DCore/v3.0.0/tutorial/afm/afm.html.
[44] J. Otsuki, K. Yoshimi, H. Shinaoka and Y. Nomura, Strong-coupling formula formomentum-dependent susceptibilities in dynamical mean-field theory, Phys. Rev. B 99,
165134 (2019), doi:10.1103/physrevb.99.165134.
[45] P. Werner, E. Gull, M. Troyer and A. J. Millis, Spin freezing transition and nonFermi-liquid self-energy in a three-orbital model, Phys. Rev. Lett. 101, 166405 (2008),
doi:10.1103/PhysRevLett.101.166405.
[46] P. Giannozzi et al., Advanced capabilities for materials modelling with Quantum ESPRESSO,
J. Phys.: Condens. Matter 29, 465901 (2017), doi:10.1088/1361-648x/aa8f79.
[47] T. Ozaki, Variationally optimized atomic orbitals for large-scale electronic structures, Phys.
Rev. B 67, 155108 (2003), doi:10.1103/physrevb.67.155108.
[48] https://issp-center-dev.github.io/DCore/v3.0.0/tutorial/srvo3/srvo3.html.
[49] Y. Nomura et al., Effective on-site interaction for dynamical mean-field theory, Phys. Rev.
B 86, 085117 (2012), doi:10.1103/physrevb.86.085117.
[50] I. A. Nekrasov, K. Held, G. Keller, D. E. Kondakov, Th. Pruschke, M. Kollar, O. K. Andersen,
V. I. Anisimov and D. Vollhardt, Momentum-resolved spectral functions of S r V O3 calculated
by LDA+DMFT, Phys. Rev. B 73, 155112 (2006), doi:10.1103/physrevb.73.155112.
[51] R. Sakuma, Ph. Werner and F. Aryasetiawan, Electronic structure of S r V O3 within
GW+DMFT, Phys. Rev. B 88, 235110 (2013), doi:10.1103/physrevb.88.235110.
[52] C. Taranto, M. Kaltak, N. Parragh, G. Sangiovanni, G. Kresse, A. Toschi and K. Held, Comparing quasiparticle GW+DMFT and LDA+DMFT for the test bed material S r V O3 , Phys.
Rev. B 88, 165119 (2013), doi:10.1103/physrevb.88.165119.
[53] https://github.com/issp-center-dev/DCore/issues.
[54] K. Yoshimi, J. Otsuki, Y. Motoyama, M. Ohzeki and H. Shinaoka, SpM: Sparse modeling
tool for analytic continuation of imaginary-time Green’s function, Comput. Phys. Commun.
244, 319 (2019), doi:10.1016/j.cpc.2019.07.001.

30

SciPost Phys. 10, 117 (2021)[55] R. Levy, J. P. F. LeBlanc and E. Gull, Implementation of the maximum entropy method for analytic continuation, Comput. Phys. Commun. 215, 149 (2017),
doi:10.1016/j.cpc.2017.01.018.
[56] H. Shinaoka, D. Geffroy, M. Wallerberger, J. Otsuki, K. Yoshimi, E. Gull and J. Kuneš,
Sparse sampling and tensor network representation of two-particle Green’s functions, SciPost Phys. 8, 012 (2020), doi:10.21468/scipostphys.8.1.012.
[57] L. Huang, iQIST v0.7:
An open source continuous-time quantum Monte
Carlo impurity solver toolkit, Comput. Phys. Commun. 221, 423 (2017),
doi:10.1016/j.cpc.2017.08.026.

31