Manual for the flexible DM-NRG code
Version 1.0.0

Team members:
Ö. Legeza
C. P. Moca
A. I. Tóth
I. Weymann
G. Zaránd1

2

Contents
1 Introduction
2

5

A short introduction to DM-NRG and to the use of symmetries
7
2.1 Wilson’s NRG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.1 The simplest example: The Kondo model . . . . . . . . . . . . . . . . . . 7
2.1.2 Extension to arbitrary symmetries . . . . . . . . . . . . . . . . . . . . . . 16
2.2 The DM-NRG method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 Installation and technical support
21
3.1 Memory and disk space requirements . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Compilation environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.3 Libraries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4 Automatic installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5 Manual installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.6 Generated binaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.7 BUG reports . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.8 Some hints for developers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.8.1 Organization of the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.8.2 On-line documentation of the code . . . . . . . . . . . . . . . . . . . . . . 24
3.8.3 Version management . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4 Using the DM-NRG code25
4.1 Main features of the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.2 Running the code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2.1 Flat density of states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2.2 Energy dependent density of states . . . . . . . . . . . . . . . . . . . . . 27
4.3 Initialization and the input file . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.3.1 <SECTION-PARAMETERS> . . . . . . . . . . . . . . . . . . . . . . . . 29
4.3.2 <SECTION-FLAGS> . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.3.3 <SECTION-SYMMETRIES> . . . . . . . . . . . . . . . . . . . . . . . . 31
4.3.4 <SECTION-BLOCK STATES> . . . . . . . . . . . . . . . . . . . . . . . 31
4.3.5 <SECTION-LOCAL STATES> . . . . . . . . . . . . . . . . . . . . . . . 31
4.3.6 <SECTION-LOCAL STATES SIGNS> . . . . . . . . . . . . . . . . . . . 32
4.3.7 <SECTION-BLOCK HAMILTONIAN> . . . . . . . . . . . . . . . . . . . 32
4.3.8 <SECTION-HOPPING OPERATORS> . . . . . . . . . . . . . . . . . . 32
4.3.9 <SECTION-SPECTRAL OPERATORS> . . . . . . . . . . . . . . . . . . 33
4.3.10 <SECTION-STATIC OPERATORS> . . . . . . . . . . . . . . . . . . . . 33
4.3.11 <SECTION-LOCAL HAMILTONIAN> . . . . . . . . . . . . . . . . . . . 34
4.3.12 <SECTION-LOCAL HOPPING OPERATORS> . . . . . . . . . . . . . 34
4.3.13 <SECTION-SPECTRAL FUNCTION> . . . . . . . . . . . . . . . . . . . 35
3

4.4

4.3.14 <SECTION-SPECTRAL FUNCTION BROADENING> . . . . . . . . .
Outputs of the DM-NRG code . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Output directory structure . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.2 Detailed description of the output files in the folder./results/automatic directory name/data . . . . . . . . . . . . . . . .
4.4.3 Analyzing the data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35
38
38
39
42

Acknowledgments

43

Bibliography

43

A Input file for the Kondo model with Ucharge (1) × Uspin (1) symmetries

47

B Input file for the Kondo model with Ucharge (1) × SUspin (2) symmetries

63

C Input file for the Kondo model with SUcharge (2) × SUspin (2) symmetries

75

D License agreements

87

4

Chapter 1

Introduction
Quantum impurity models describe interactions between some local degrees of freedom (e.g. aspin) and a continuum of non-interacting fermionic or bosonic states. The investigation of quantum impurity models is a starting point towards the understanding of more complex strongly
correlated systems, but quantum impurity models also provide the description of various correlated mesoscopic structures, biological and chemical processes, atomic physics and describe
phenomena such as dissipation or dephasing. Prototypes of these models are the Anderson
impurity model, or the single- and multi-channel Kondo models. The first two models are classic examples of Fermi liquid models, while the multi-channel Kondo model is the most basic
example of a non-Fermi liquid system, and as such, it serves possibly as the simplest realizationof a quantum critical state.
The solution of these models for low energies was a major issue in theoretical condensed matter
research and led to the development of various non-perturbative techniques. [The interested
readers are referred to the seminal book of Hewson [1] and to the extensive review of Cox and
Zawadowski [2].] However, despite the extensive effort, many of the methods developed areuncontrolled, while others can be applied only to a subclass of models or to restricted regionsof the parameter space. Wilson’s numerical renormalization method, originally developed for
the Kondo model, remained possibly the most reliable method to study dynamical correlationsof generic quantum impurity models as well as their thermodynamic properties and finite size
spectra. It is still one of the most popular methods to study quantum impurity models.
Although Wilson’s NRG has been used in its original form for a longtime, a number of new developments took place recently: First, a spectral sum-conserving density matrix NRG approach(DM-NRG) has been developed [3, 4], which has very recently been generalized for non-Abelian
symmetries [5]. We remark that using symmetries as much as possible is necessary in manycases to perform accurate enough calculations. NRG has been also restructured as a matrix
product state approach [6]. These new developments not only made NRG much more reliable
than Wilson’s original method [5, 7], but they made it possible to extend NRG to study nonequilibrium phenomena [8, 9], and opened the way to use methods familiar from the density
matrix renormalization group (DMRG) approach [10].
In this manual we do not intend to give a complete description of the NRG machinery. Rather,
we introduce some of the basic concepts that are needed to use NRG and to use the code we
provide for quantum impurity problems of interest. If you want to understand in detail, how
NRG and DM-NRG work, we advise you to consult Wilson’s original work [11, 12], and thereferences listed above.
The code we describe in this manual is a free density matrix numerical renormalization group(DM-NRG) code, which can be downloaded from the site http://www.phy.bme.hu/∼dmnrg.
This code is a flexible NRG code, which uses user-defined non-Abelian symmetries dynamically,
computes spectral functions, expectation values of local operators for user-defined impurity
5models. The code can use a uniform density of states as well as a user-defined density of states.
The current version of the code assumes fermionic bath’s. It uses any number of U (1), SU (2)charge SU (2) or Z2 symmetries, but the interested user can teach the code other symmetries
too. The code runs using a simple input file. We provide several example input files with thecode as well as a few Mathematica files with which these input files can easily be constructed.
An energy spectrum analyzer utility is also provided with the code to study finite size spectratoo.

6Chapter 2

A short introduction to
DM-NRG and to the use ofsymmetries
In this chapter, we give a short overview of Wilson’s numerical renormalization group (NRG)and the density matrix NRG (DM-NRG). As already mentioned in the introduction, here we
introduce only the basic concepts that are needed to use NRG and to use the code we provide
for quantum impurity problems, but we do not attempt/intend to give a complete overview ofthe existing literature. To learn more details about NRG and DM-NRG, we recommend to read
Wilson’s original work [11, 12].

2.1

Wilson’s NRG

2.1.1

The simplest example: The Kondo model

Before introducing the general concepts, let us discuss the simplest possible example, the Kondo
model. The Kondo model consists of a spin S interacting locally with a non-interacting conduction electron sea,J~X †HKondo = S(2.1)
ψσ ~σσ,σ′ ψσ′ + Hcond .
2′σ,σHere J denotes the Kondo coupling, ψσ† creates a conduction electron at the impurity site, and~σσ,σ′ is the vector of Pauli spin matrices. The term, Hcond , describes the conduction electron
bath, and its specific form is not very important for us.
Wilson’s mapping and iterative diagonalization
From the point of view of local dynamics, Hcond contains a lot of redundant information: In fact,
since our fermions are non-interacting, the local density of states ̺(ω), i.e. the spectral functionof the unperturbed Green’s function, Gψσ ,ψσ† (ω), is generated by Hcond with J = 0. Note that
̺(ω) determines completely the correlation functions of ψσ† and the spin dynamics even for theinteracting system, J 6= 0. The ingenious idea of Wilson was to discretize logarithmically ̺(ω)
using a discretization parameter Λ > 1, and thus map approximately the Hamiltonian (2.1) toa semi-infinite chain (see Fig. 2.1 and Fig. 2.2):
7

ρ(ω)

−DΛ

−1/2

−D

0

....

−1/2

DΛ

....

t0

t1

D

ω

t2
...

ξ1
ξ2

ξ0

S

ξ3

...

J

Figure 2.1: Discretization of the conduction band density of states on a logarithmic mesh andmapping onto the Wilson chain (see also Fig. 2.2).

W ilsonHKondo
=

∞ X
∞ X
X
X
J ~X †††
S
f0,σ ~σσ,σ′ f0,σ′ +ξn fn,σfn,σ +tn (fn,σfn+1,σ + h.c.) . (2.2)
2
′
n=0 σ
n=0 σ
σ,σ

†
Here the spin interacts only with the fermion f0,σat the end of the Wilson chain. The operator†
f0,σ creates a conduction electron at the impurity site, and is essentially identical to the operator
ψσ† . The on-site energies ξn and the hoppings tn depend solely on ̺(ω) and can be determinedrecursively [11]. Our code takes care of this part of the work: it determines numerically these
constants if a density of states ̺(ω) is provided (see section 4.4.3 and the description of theutility he).
For a flat and symmetrical density of states, ̺(ω) = 1/(2D), one has ξn = 0 and the tn ’s canbe determined analytically ( tn ∼ Λ−n/2 ) [11]. Longer and longer chains give more and moreaccurate description of the infinite chain. This observation led Wilson and his co-workers tosolve the Hamiltonian (2.2) iteratively. Introducing the operator
n X
n X
X
X
J ~X †††′′
Hn ≡ Sξm fm,σ fm,σ +tm (fm,σfm+1,σ + h.c.) .
f0,σ ~σσ,σ f0,σ +
2
′
m=0 σ
m=0 σ

(2.3)

σ,σ

one has the obvious recursion relation,
Hn+1 = Hn + τn,n+1 + Hn+1 ,

(2.4)

P
P
†
†
fn+1,σ (seefn+1,σ + h.c.) and Hn+1 = σ ξn+1 fn+1,σwith the notation, τn,n+1 = tn σ (fn,σalso Fig. 2.2). Then Wilson’s procedure consists of constructing from the low-energy eigenstates,|uin , of the operator Hn approximate eigenstates, |ũin+1 , of the operator Hn+1 . To do this,
one takes a definite number of the states |uin with the lowest energies and generates new states†from them by first adding an empty site, |uin → |uin+1 and then using the operators fn+1,σto8

τ0,1 τ1,2
H0 H1

τ N−1,NH2

HN−1 HN

H1
H2
Figure 2.2: Wilson chain of length N. The impurity site is at the origin and is represented bya square, while the sites are represented as circles. For the impurity the Hamiltonian H0 is
identical with H0 . As we add sites the iterative construction of the Hamiltonian at a given site
Hn is pictured.
create new states as|uin+1
 f†
|uin+1|uin →  n+1,↑
 .
† fn+1,↓ |uin+1
††
fn+1,↑
fn+1,↓|uin+1(2.5)

Then one diagonalizes the Hamiltonian Hn+1 in this new basis set. To construct the matrix
elements of Hn+1 in this new basis, one needs the following information
1. The eigenvalues Eun of Hn ,2. The matrix elements of Hn+1 between the local states, |µi, constructed from the vacuum
state |0i as


|0i

 f†
|0i,


(2.6)
{|µi} ≡  n+1,↑
†
 fn+1,↓ |0i††
fn+1,↑
fn+1,↓|0i†
3. The matrix elements of fn+1,σbetween the local states, |µi, and†
4. The matrix elements of fn,σbetween the block states |uin .

Diagonalizing the Hamiltonian Hn+1 one then obtains the new eigenstates |ũin+1 , their eigen†values, Eũn+1 , and one can trivially compute the matrix elements of fn+1,σbetween these states
too, to proceed to the next iteration.
If one is to compute the spectral function of a local operator A acting at site n = 0, then one
has to keep track of the matrix elements n hu|A|vin of this operator, too. Already the finite size
spectrum, i.e. the spectrum of Hn contains a lot of precious information [11]. However, oncethe matrix elements and the approximate spectrum of the Hamiltonian is at hand, one can goahead and also compute thermal expectation values or spectral functions from it [13, 14], andthus determine completely the properties of a quantum impurity problem.
9Symmetries in the Kondo model
The total spin operator∞
X~T = S
~+1
S
2 n=0

X

†
fn,σ
~σσ,σ′ fn,σ′(2.7)σ,σ′ ∈{↑,↓}

as well as the charge operator,
Q=

∞
X

n=0




X

σ={↑,↓}

†fn,σfn,σ − 1(2.8)

commute with Eq. (2.2). This implies that the eigenstates of the Hamiltonian can be classified
as multiplets. Every multiplet u will be characterized by a charge quantum number, Q(u) anda spin quantum number, S(u). We shall refer to these as representation indices. Internal states
within a multiplet are labeled by the z-component of the spin, that we shall refer to as labels.
Similarly, operators can also be characterized by quantum numbers. To give an example, in this††simple case, the components of {fn,↑, fn,↓ } and {fn,↓, −fn,↑ } transform under spin rotations as
the components | ± 1/2i of a spin Sf = 1/2 state, and have charges Qf † = 1/2 and Qf = −1/2,√√
respectively. The three operators formed from the impurity spin, Sm ≡ {−S + / 2, S z , S − / 2}
transform, on the other hand, as the three components of a spin SS = 1 state, |m = 0, ±i, while
they have trivially charge QS = 0. The operators discussed before provide simple examplesof irreducible tensor operator multiplets {A}. Notice the non-trivial prefactors in the previous
examples. These prefactors must always be worked out carefully, since the smallest sign mistake
may alter the NRG results significantly.
In the presence of symmetries, a powerful theorem, the Wigner–Eckart theorem tells us that,in order to compute the matrix elements of a member of an operator multiplet {A} between
two states within multiplets u and v a single matrix element is needed, the so-called reduced
matrix element [5]:hu||A||vi .
(2.9)
Wilson’s procedure thus gets a bit modified. First of all, one needs to classify not only theblock states (multiplets) but also the local states µ added at site n + 1 by symmetries. Alongthe iteration, one constructs from these new states using group-theoretical methods (Clebsch~ 2 and S z . Apart from this twist, theGordan coefficients), which are already eigenstates of S
T
discussion of the previous subsection is still valid, and gets only slightly modified. We can thus
make the following statement: To construct the (reduced) matrix elements of hi||Hn+1 ||ji inthis basis, one needs the following information
1. The eigenvalues Eun of Hn ,2. The matrix elements hµ||Hn+1 ||νi between the local multiplets.
†
3. The reduced matrix elements, hµ||fn+1||νi,

4. The reduced matrix elements n hu||fn† ||vin between the block multiplets, |uin .
5. Finally, to compute a spectral function of an operator multiplet A, one also needs n hu||A||vin .
At this point, the reader does not need to know in detail, how the construction of Hn+1 takes
place, since the code does that automatically, but for more details, the interested reader isreferred to Ref. [5].
10BlockstatesQ

Sz|1i = | ⇓i

−1

− 21

|2i = | ⇑i

−1

1
2

†
|3i = f0,↓
| ⇓i

0


†
†
|4i = √12 f0,↓
| ⇑i − f0,↑
| ⇓i

Added/local
statesqsz−1|1i = |0i

−1

0

0

0

†
|2i = fn+1,↓
|0i

0

− 21


†
†
|5i = √12 f0,↓
| ⇑i + f0,↑
| ⇓i

0

0

†
|3i = fn+1,↑
|0i

0

1
2

†
|6i = f0,↑
| ⇑i

0

1

†
†
|4i = fn+1,↑
fn+1,↓|0i1

0†|7i = f0,↑
f↓† | ⇓i

1

− 21

†
|8i = f0,↑
f↓† | ⇑i1

1
2

Table 2.1: Block states (left) and local/added states (right) for the single channel Kondo model
when Ucharge (1) × Uspin (1) symmetry is used. The first block states are formed from theimpurity spin and the conduction electron at site n = 0 of the Wilson chain.

Kondo model with Ucharge (1) × Uspin (1) symmetry

Let us first illustrate the concepts introduced above through the simplest case, when every
symmetry used is Abelian. In this case every multiplet consists of a single state, which ischaracterized by some quantum numbers. Let us thus consider the Kondo model with two
U (1) symmetries, generated by the z-component of the total spin, Sz and the charge operators,
Q. Now Sz and Q, provide the quantum numbers according to which the multiplets of the
Hamiltonian are classified. The classification of the states in this case is rather trivial andthe computation of the reduced matrix elements is also simple, since all the Clebsch-Gordan†
coefficients are 1 or 0. However, without spin SUspin (2) symmetry, the hoppings of fn,↑
and
†
fn,↓
are not related by symmetry. There are therefore two hopping operators which we need to††keep track of, fn,↑
and fn,↓.

The initial block states and the matrix elements of the hopping operators can be constructed
immediately. A Mathematica file, kondo U1 charge U1 spin.nb, with all the numerical detailsof this calculation as well as an input file is provided by the package. The corresponding blockand the local states are listed in Table 2.1.
In the first iteration, the reduced matrix elements of the Hamiltonian matrix between the block11

states (Table 2.1) are given by


0
 0

 0

 0

hu||H
||vi
=
J
0
0
0
 0

 0

 0
0


0
0 

0 

0 
 .
0 

0 

0 
0

0 0
0
0
0 0
0 0
0
0
0 0
0 1/4
0
0
0 0
0 0 −3/4 0
0 0
0 0
0
1/4 0 00 0
00 1/4 0
0 0
0
0
0 0
0 0
0
0
0 0

(2.10)

In general there might be off-diagonal matrix elements in this matrix. The code takes care ofthis and diagonalizes 0 hu||H0 ||vi0 before proceeding to the next iteration.
†
The matrix elements of the f0,↑
are


0
0 0
0
0
0 0 0

0
0 0
0
0
0 0 0 



0
0
0
0
0
0 0 0 √
 −1/ 2 0 0

0
0
0
0
0
†

 ,
√
(2.11)
0 hu||f0,↑ ||vi0 = 00
0 0 0 
 1/ 2 0 0


0
1 0
0
0
0 0 0 



0
0 1
0√
0√ 0 0 0 
0
0 0 1/ 2 1/ 2 0 0 0†
while those of f0,↓
are given as


0
0
 0
0

 1
0√

 0 1/ 2†√
0 hu||f0,↓ ||vi0 =  0 1/ 2
 0
0

 0
0
0
0

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0√
0√
0
0 1/ 2 −1/ 2 0
0
0
0
−1

0
0
0
0
0
0
0
0


0
0 

0 

0 
 .
0 

0 

0 
0

(2.12)

These matrix elements are needed to construct the hopping terms in the first iteration. In††
addition, one also needs the matrix elements of the hopping operators fn+1,↑
and fn+1,↓
betweenthe local states, listed in Table 2.1. These are given by
0 0 0 00 0 0 0
 1 0 0 0 
 0 0 0 0 ††
(2.13)
hµ||fn+1,↓||νi = hµ||fn+1,↑||νi =  0 0 0 0  .
 1 0 0 0  ,0 0 −1 0
0 1 0 0

Finally, if we want to determine the auto correlation function of Sz , e.g., then we also need tocompute the matrix elements this operator,−1/2 0
0
0
0
0
0
0
 0
1/2
0
0
0
0
0
0 


 0
0 −1/2 0
0
0
0
0 


 0
0
0
0
1/2 0
0
0 

 .
hu||S
||vi
=
(2.14)
0
z
0
 0
0
0
1/2 0
0
0
0 


 0
0
0
0
0 1/2
0
0 


 0
0
0
0
0
0 −1/2 0 
0
0
0
0
0
0
0
1/2
12

BlockstatesQS|1i = | ⇑i

−1

1
2

†
†
| ⇑i − f0,↑| ⇓i)
|2i = √12 (f0,↓00†|3i = f0,↑
| ⇑i

0

1

†
†
|4i = f0,↑
f0,↓
| ⇑i

1

1
2

Local
statesqs|1i = |0i−10†|2i = fn+1,↑
|0i

0

1
2

†
†
|3i = fn+1,↑
fn+1↓
|0i

1

0

Table 2.2: Block states (left) and local/added states (right) for the single channel Kondo model
when Ucharge (1) × SUspin (2) symmetry is used. The first block states are formed from theimpurity spin and the conduction electron at site n = 0 of the Wilson chain.

We remark at this point that, while Sz is an irreducible tensor operator and it has a spinquantum number 0, the operators Sx and Sy are not irreducible tensor operators with respect
to the Uspin (1) spin rotations, since they transform among each other. They are, however,linear combinations of the irreducible tensor operators, S± , which have spin quantum numbers±1 under rotations around the z-axis.
Once all this information typed into the input file, input.dat, the flexible DM-NRG code is
ready to compute the finite size spectrum, the spin’s spectral function, and the real part of
the retarded spin-spin correlation function using both traditional NRG as well as DM-NRG
methods! An example input file for this case is provided in the directory input files.

Kondo model with Ucharge (1) × SUspin (2) symmetry
Let us now show on the example of the Kondo model, how non-Abelian symmetries can beused. Let us thus enumerate all necessary information that is needed to perform a calculation
for the spin spectral function of the single channel Kondo model with Ucharge (1) × SUspin (2)symmetry. A Mathematica file (kondo U1 charge SU2 spin.nb) where these inputs have been
computed is also provided by the package.
In the first iteration we have four different block multiplets, u = 1, . . . , 4 formed from theimpurity spin and the conduction electron at site n = 0, while there are three added multiplets
in each iteration, µ = 1, 2, 3. These states and their quantum numbers are listed in Table
2.2 (only highest weight states are given). The reduced matrix elements of the Hamiltonian
between the block states then read


0
 0

0 hu||H0 ||vi0 = J 
0
0


0
0 0
−3/4 0 0  .
0
1/4 0 0
0 0

(2.15)††In this case, the matrix elements of f0,↑
and f0,↓are related with each other by symmetry, andtherefore they form a single hopping operator multiplet f0† of spin 1/2. The matrix elements of13

BlockstatesQ

S††| ⇑i − f0,↑| ⇓i)
|1i = √12 (f0,↓0†|2i = f0,↑
| ⇑i

0

1

†
†
|3i = f0,↑
f0,↓
| ⇑i

1
2

1
2

Local
statesqs†|1i = fn+1,↑
|0i

0

1
2

†
†
|2i = fn+1,↑
fn+1,↓
|0i

1
2

0

0

Table 2.3: Block states (left) and local/added states (right) for the single channel Kondo model
when SUcharge (2) × SUspin (2) symmetry is used. The first block states are formed from theimpurity spin and the conduction electron at site n = 0 of the Wilson chain.

f0† need to be determined using the Wigner–Eckart theorem, and are given by


0
 1
†

0 hu||f0 ||vi0 = 
1
0

0
0
0√
1/ 2

0
0
p0
− 3/2


0
0 
 .
0 
0

(2.16)†
To generate the hopping, we also need to know the reduced matrix elements of fn+1
betweenthe added (local) states. These are computed as†hµ||fn+1||νi 

0
1
0

0
0
√
− 2


0
0  .
0

(2.17)

Finally, to compute the spin-spin correlation function, we need to compute the reduced matrix
elements of the impurity spin. As discussed before, the components of the spin operator aregrouped in this case into a single tensor operator, S. The reduced matrix elements of thisoperator are computed as
 p

3/4 0
0
p0
 0
3/4
0 
0 − √

 .
0 hu||S||vi0 = 
0
1/2 1/ 2 p0 
0
0
0
3/4

(2.18)

To run the code with Ucharge (1) × SUspin (2) symmetry, one only needs to type this information
into the input file, input.dat, and then run the code. You can find the corresponding input
file, kondo model U1 c SU2 s.dat in the directory input files.
14

Kondo model with SUcharge (2) × SUspin (2) symmetry
The Hamiltonian (2.2) is invariant under the action of SUcharge (2) rotations in charge space,~
Uc = ei~ωQ Q , generated by the operators Qx = (Q+ + Q− )/2, Qy = (Q+ − Q− )/2i and Qz with
Q+

=

∞
X

†
†
(−1)n fn,↑
fn,↓
,

n=0

Qz
Q−

=
=

∞
1X



Q+

,

2 n=0



†

X

µ={↑,↓}

†fn,µfn,µ − 1 ,(2.19)and parametrized by real, three-component vectors ω~ Q . Since the spin symmetry generators
commute with the charge symmetry generators, the Hamiltonian (2.2) possesses a symmetry,
SUcharge (2) × SUspin (2).
Charge SU(2) symmetries must always be used with care, and one must always carefully check if
a given set of operators forms an irreducible tensor operator. One can check, e.g., that the cor††
rectly defined hopping operators read, γn† = {(γn )στ } = {fn,↑, fn,↓, (−1)n fn,↓ , (−1)n+1 fn,↑ },and that they indeed transform as a set of spin S = 1/2 and charge spin Q = 1/2 irreducible
tensor operators. Notice the curious sign-dependence of the last two operators and their relative
signs. Similar factors of (−1)n appear in the states generated from the highest weight states inTable 2.3 through the action of the operator Q− .
In the first iteration we have three different block multiplets, u = 1, . . . , 3 formed from theimpurity spin and the conduction electron at site n = 0. There are two added (local) multiplets
in each iteration, µ = 1, 2. These states and their quantum numbers are listed in Table 2.3 (onlyhighest weight states are given). The reduced matrix elements of the Hamiltonian between theblock states then read
−3/4 0 0
 0
1/4 0  .
(2.20)
0 hu||H0 ||vi0 = J
0
0 0
The matrix elements of γ0† need to be determined using the Wiegner-Eckart theorem, and are
given by
√ 

0
0
−√2
 0
†
0
− 2 
(2.21)
 .
0 hu||γ ||vi0 = 
q
0

√1
2

−

3
2

0

To generate the hopping, we also need to know the reduced matrix elements of γ1† between theadded (local) states. These are computed as
√ 

2
0
√
.
(2.22)
hµ||γ1† ||νi =− 2 0Finally, to determine the spin spectral function, we need to compute the reduced matrix elements of the impurity spin operator:
p


3/4
0
0 − √
 1/2 1/ 2
 .
(2.23)
0 hu||S||vi0 =
p0
3/4
0
0

To use the code with SUcharge (2)×SUspin (2) symmetry, one only needs to type this information
into the input file, input.dat, and then run the code. You can find the corresponding input file,
kondo model SU2 c SU2 s.dat in the directory input files, and a Mathematica notebook tocompute these matrix elements is also provided.
15

2.1.2

Extension to arbitrary symmetries

The procedure outlined above can be generalized to essentially any number and type of nonAbelian symmetries and to any Hamiltonian of the form:
H = H0 +∞Xn=0(τn,n+1 + Hn+1 ) .

(2.24)

The first term in Eq. (2.24), H0 , describes the quantum impurity coupled to the ‘environment’.
This is thus the ’interaction part’ of the Hamiltonian. The dynamics of the environment is
described by hopping terms that connect nearest-neighbors only, τn,n+1 , and the on-site terms
Hn . As already stated before, these on-site terms can account for the breaking of electron-holesymmetry or for the presence of superconducting correlations [15].
The above Hamiltonian can also be diagonalized iteratively using the recursion relation, (2.4).
The procedure is essentially the same as the one outlined in the previous subsection, with themodification that we now may have Γ different symmetries, which commute with the Hamiltonian and each other, and which we shall label by an index γ = 1, . . . , Γ. Correspondingly,
multiplets will be labeled by Γ different quantum numbers, {Qγ } = {Q1 , . . . , QΓ }, with Qγ thelabel of the irreducible representation of symmetry γ. States within a degenerate multiplet canthen be labeled by internal quantum numbers, which, by analogy to the spin operators, we shalldenote by Qγz . Finally, the representation indices and the internal quantum numbers can bothbe grouped into vectors, Q = {Q1 , Q2 , ..., QΓ } and Qz = {Q1z , Q2z , ..., QΓz }.
Similar to multiplets, a set of irreducible tensor operators, {A} is also characterized by somequantum numbers a, and az . The hopping terms, τn,n+1 , can also be decomposed using a setof such irreducible operators, which we shall call ‘hopping operators’, f † → C, and label themby λ. Thus the hopping operator multiplet, Cλ has quantum numbers cλ and members of thismultiplet are labeled by czλ . Ordering the members of a hopping operator multiplet into vectors,C λ ≡ {Ccλ ,czλ }, the hopping can be simply written asX(2.25)
τn,n+1 =tn,λ C n,λ · C †n+1,λ + h.c. .
λNote that here different hopping operator multiplets may have different hopping amplitudes,
and also, remember that now the C n,λ ’s denote creation operators. An example for having two
operators has already been for the Kondo model with Ucharge (1) × Uspin (1) symmetries, but inmost cases several hopping operators must be used for multichannel models too.
In our general framework the iteration step n → n + 1 reads now as follows:• Previous block states u can be characterized by their symmetry indices and labels, Qu , Qzu ,Eu, Qu , Qzu . Similarly, local states µ are also characterized by their symmetry indicesnEγq µ = {qµγ }, and internal labels q zµ = {qµ,z}: µ, qµ , q zµ .

• The new basis states can be constructed using Clebsch-Gordan coefficients as
E
E
E
X D
=
q µ q zµ ; Qu Qzu Qi Qzi u, Qu , Qzu ; µ, q µ , q zµi, Qi , Qzi
,
n

n−1

Qzu ,qzµ(2.26)

where we have also introduced a general notation for the Clebsch–Gordan coefficients as
a product of coefficients, one for each symmetry in the problemDq µ q zµ ; Qu Qzu Qi Qzi

E

=

Γ
Y

γ=1

16

γ
qµγ qµ,z; Qγu Qγu,z Qγi Qγi,z.

(2.27)

• Knowing the reduced matrix elements, n hu k Cn,λ k vin and hµ k Cn+1,λ k νi, and theeigenvalues Eun from iteration n one can construct the Hamiltonian in the new basis.
E• Diagonalizing this Hamiltonian, one can determine the new eigenstates, ũ, Qũ , Qzũn+1

their energies, Eun+1 , and the irreducible matrix elements, n+1 hũ k Cn+1,λ k ṽin+1 . If aspectral function is computed then the corresponding operator’s reduced matrix element,
n+1 hũ k A k ṽin+1 must also be computed.

• After the iteration some of the high-energy states are also discarded before one proceeds
to the next iteration.
Our flexible DM-NRG code implements this rather general procedure dynamically, and in sucha way that various symmetries can be taught to the code. If you want to teach a new symmetry
to the code, feel free to contact us for help.

2.2

The DM-NRG method

Wilson’s procedure gives accurate results for the finite size spectrum, and it can also be used
to compute spectral functions [13]. However, there are several problems with the simple extension of Wilson’s method. On one hand, in the original approach, spectral functions havebeen computed using only the approximate ground states in a given iteration. This may giveincorrect results if the fine structure of the ground state shows up at subsequent iterations only.
Hofstetter proposed to use a density matrix NRG (DM-NRG) method to treat this problem[7]. Furthermore, in the original method one simply chops off parts of the Hilbert space, andtherefore spectral sum rules are not respected. This problem has been solved only recently,
through the introduction of a complete basis set [8]. In this approach, rather than constructing
approximate eigenstates for a chain of length n, one constructs approximate eigenstates of Hnthat live on a chain of length N > n, by adding environment statesEEu, Qu , Qzu → u, Qu , Qzu ; e(2.28)
nn

where e labels all possible states of the environment, i.e., chain sites from site n + 1 to N .
The reduced density matrix
Physical quantities are then computed from the density matrix, which is approximated as
̺

=

N
X

̺n ,

(2.29)

n=0

̺n

=

1 X
Z e

X

X

u∈Discarded Qzu

ne−βEu u, Qu , Qzu ; eEnnwhere the partition function is defined as
X
n
−n
Z=
dimN
dim(u) e−βEu ,loc

Du, Qu , Qzu ; e ,(2.30)(2.31)uP
where dim(u) is the dimension of the multiplet u, and dimloc = µ dim(µ) denotes the dimension of the local basis at a given site of the chain. The second summation in Eq. (2.31) runs
over states that have been discarded in iteration n → n + 1.
17

The calculation of a spectral function then consists of a forward sweep, where we determine all
eigenstates in a given iteration, and a backward sweep, where we determine ̺ and the truncatedreduced density matrices, defined as
X
R[n] ≡ TrN −n {̺m } .
(2.32)
m≥n

The details of the corresponding recursion relations were given in Ref. [5].
Spectral function calculations in the flexible-DM-NRG framework
The main purpose of the code is to compute the retarded Green’s function of two local operators,
A and B. This is defined for two irreducible tensor operators as
iE
Dh
†
(0)
Θ (t) .
(2.33)
(t)
=
−i
A
(t),
B
GR
†
a,a
b,b
z
,B
A
a,az

z

b,bz

where the quantum numbers a, az and b, bz refer to the quantum numbers of operators A† andB † . For fermionic operators one must replace the commutator above by an anticommutator. By
symmetry, this Green’s function is non-zero only if the spectral operators A and B transform
accordingly to the same representation,GRA†a,az ,Bb,bz(t) = GRA,B † (t) δa,b δaz ,bz .

(2.34)

Although it is a cumbersome expression, for completeness, let us give here the expression
obtained generalizing the procedure of Ref.[8].
GRA,B † (z) =
×

"

N
X
X

X

n=0 i∈D,K (j,k)∈(K,K)
/
n

∗

n

hi k R[n] k ji n∗

hj k B † k ki n n hi k A† k ki n dim(i)hk k A† k ji n n hk k B † k ii n dim(k)n−ξdim(a)
dim(a)
z + 21 (Ein + Ejn ) − Ekn
z − 21 (Ein + Ejn ) − Ekn#.

(2.35)

Remarkably, this formula contains exclusively the reduced matrix elements and the dimensionsof the various multiplets and operator multiplets. Here the second sum is over all the multiplets
i, j, k of the given iteration subject to the restriction that j, k do not belong to kept states at thesame time. No summation is needed for states within the multiplets. The sign factorQΓ ξ is ξ = 1
for bosonic operators, while it is ξ = −1 for fermionic operators, and dim(a) = γ=1 dim(aγ )is the dimension of the operator multiplet A†a,az . The code computes the spectral function
using this formula, and the full Green’s function is then computed by doing a Hilbert transformnumerically.
Smoothing procedure
In general, the spectral functions are given as weighted sums of δ-functions of the form:
AA,B † (ω) ==1+− ℑm GRA,B † (ω + i 0 )πX
wi δ(ω − ωi )(2.36)(2.37)iwhere the weighting coefficients wi can be expressed from Eq. (2.35) after analytical continuation to the real axis. To get a smooth spectral function, however, the delta functions need to18

be replaced by some broadening kernels K(ω, ωi ) that fulfill the conditionZdω K(ω, ωi ) = 1 .

(2.38)

This condition guarantees that spectral sum rules remain valid after the broadening procedure.
In our code we use two types of kernels to generate the real axis spectral function: (a) the
log-Gaussian kernel [16] and (b) the interpolative kernel that we constructed to avoid certain
difficulties with the interpolation around ω = 0.
For the log-Gaussian method the following kernel is used
2
2
2 1
1
.
Klog (ω, ωi ) ≡ √ e−b /4 e−(ln ω−ln ωi ) /bb πωi

(2.39)

This kernel works rather well at T = 0 calculations, however, it extrapolates to the ω → ±0
values in a singular way. Therefore, it has problems at finite temperatures and also in cases
where spectral losses lead to an artificial jump in the spectral function at ω = 0.
The interpolative scheme avoids this difficulty by using a kernel
22 dx
1
Kint (ω, ωi ) ≡ √ e−[x(ω)−x(ωi)] /b
,
b π
dωi

(2.40)

with the smoothening function defined as1x(ω) = tanh
2



ω
TQ



ln



ω
TQ

2

+ eγ

!

,

(2.41)by choosing γ ≈ TQ e.g. This formula interpolates smoothly along the real axis for all frequencies. It has the property that for frequencies smaller than the ’quantum temperature’,|ω| ≪ TQ , the broadening is Gaussian, while for |ω| ≫ TQ it becomes log-Gaussian. In thisscheme, for zero temperature calculations the quantum temperature TQ must be the smallest
energy scale in the problem, on the other hand for finite temperature calculations, the quantum
temperature must be set to be in the range of the temperature itself, TQ ∼ 0.5T ÷ T . In bothinterpolative methods there is a fitting parameter b, which controls
the widths of kernels and in√general is Λ dependent, which should be in the range of b ∼ Λ/2 (for Λ = 2, b = 0.5 ÷ 0.9).

19

20

Chapter 3

Installation and technical support
3.1

Memory and disk space requirements

The resources required by the DM-NRG code depend a lot on the type of the ongoing calculation.
Increasing the number of block states that are kept during a run, the CPU time needed for
performing the calculation, the necessary memory, as well as the used disk space may increase
substantially. As an example in Table 3.1 we present some estimates for a simple model. The
following information was extracted by running the code on an Intel(R) Core(TM)2 Duo CPU
T7250 @ 2.00GHz processor with 4MB CPU cache. We have used the CentOS 5.1 distribution
for the Linux operating system. The code was run by using an input file for the Andersonmodel with U(1) × SU(2) symmetries that comes with the distribution. We made a run for aflat band and performed 50 iterations. The compiler used was gcc version 4.1.2.
Model name
Anderson model
U(1) × SU(2)symmetries

number ofblock states
50
100
200
500

CPU-typeIntel Core(TM)2 Duo
CPU T7250 @ 2.00GHz
with 4MB cache

CPU time
(sec)
1.27
3.56
14.39
156.26

disk space
(MB)
12
25
69
341

memory
(MB)
4
7
12
18

Table 3.1: CPU time, memory and disk space requirements.
Important note: when temporary files are also generated in ASCII format by setting text swap files flag
= ON in <SECTION-FLAGS> of the input file input.dat, approximately four times larger
disk space is required than the ones shown in the table.

3.2

Compilation environment

The code can be compiled with either g++ or with icc (intel c++ compiler). The Intel compiler, icc is available for free for non-commercial and non-academic use after signing a license
agreement.

3.3

Libraries

The Flexible-DM-NRG code is linked with state of the art libraries. The following list of
libraries (packages) are needed to have the flexible-dmnrg package properly installed. The
easiest way to install these packages is by using the rpm - RPM package manager method.
21

For further information and links to the websites from where these libraries can be downloaded
please refer to the web page of the code http://www.phy.bme.hu/∼dmnrg/index.html.
List of required libraries:• libg2c – Fortran 77 Libraries• blas – Basic Linear Algebra Subprograms• lapack – The LAPACK libraries for numerical linear algebra.
• gmp – A GNU arbitrary precision library• gsl – GNU Scientific Library• gsl-devel – GNU Scientific Library - development files

3.4

Automatic installation

Once the libraries were successfully installed, run the following shell script in the main directoryof the code:./configure
This command must run without errors. It checks your system, searches for the necessary
libraries, and sets the library paths. An error message returned while running ./configuremeans that some packages might be missing or not properly installed. If the ./configure has
run without errors then the make.sys and make.rules and config.h files have been generated.
These files now contain the system variables needed for the installation 1 .
Next run the following command:
make allIf everything goes smooth and no error messages were generated, then the code is successfully
installed. If for some reasons the automatic configuration does not work, please, follow the
instructions in the next subsection and install the code manually. If you do not succeed, feel
free to contact us.

3.5

Manual installation

Sometimes, for some mysterious reasons, the autoconf does not work properly and some of thelibraries are not detected. In this case, the make.sys file, that can be found in the main folderof the code, needs to be edited manually. Altogether there are ten shell variables that need tobe edited.
• CXX - the name of the C++ compiler, is usually g++ or icc.
• C - the name of the C compiler, usually gcc.
• CC FLAGS - flags for the g++ or icc compiler (this variable can be left empty at this stage).
You can set here the optimization level to speed up the code.
• GSL FLAGS - path to the gsl include files, usually -I/usr/include.
• GMP FLAGS - path to the gmp include files, usually -I/usr/include.
• FLIBS - fortran libraries, usually g2c or gfortran.
• LAPACK LIBS - lapack library, usually -llapack. For static libraries the full path name needs
to be provided sometimes: /lib/liblapack.a.
• BLAS LIBS - blas library, usually -lblas. For static libraries the full path name needs to beprovided sometimes: /lib/libblas.a.
• GSL LIBS - gsl library, usually -lgsl -lgslcblas. For static libraries the full path nameneeds to be provided sometimes: /lib/libgsl.a.
• GMP LIBS - gmp library, usually -lgmp. For static libraries the full path name needs to beprovided sometimes: /lib/libgmp.a.
1 The automatic installation checks for the g++ and gcc compilers and does not check for the icc. To run thecode with icc, manual intervention is necessary (see section Manual Installation), and the CXX environmental
variable needs to be modified from g++ to icc.

22

Once the make.sys file is properly edited run
make allIf the compilation runs without errors then the code is successfully installed and all the binaries
have been generated.

3.6

Generated binaries

After a successful installation the following binaries are generated:• fnrg - the main program that does the NRG/DM-NRG calculations.
• sfb - utility that does the spectral function broadening.
• he - utility that computes the hopping for an energy-dependent density of states.
• es - utility to generate/analyze the energy spectrum.
• cgc - utility to generate/compute the Clebsch-Gordan coefficients
To run the code the input.dat file must be set properly. For further details, refer to thenext sections of the manual. Example input files are provided with the code for the single
band Anderson model with U (1) charge and spin symmetries and for the Kondo model withcombinations of U (1) and SU (2) symmetries. These ’ready to run’ input files are available in the./input files folder. Furthermore, we also provided the Mathematica files which were used
to compute the content of these input files. Some of the technical details of these calculations
were summarized in Chapter 2.

3.7

BUG reports

We have made this code publicly available, so if you think you have found a bug in the code,
please let us know so we can fix it for future releases.
Please send any bug report to: dmnrg@neumann.phy.bme.hu. Use BUG REPORT in thesubject area.
Please include the followings in any report:• The version number of the code.
• input.dat file so we can test the code with it.
• A description of what is wrong. If the results are incorrect, in what way. If you get a crash,etc.
• Please do not send core dumps, of executables.
• The name of the compiler that was used and its version.
• The output from running ./configure• If the bug is related to configuration, then the contents of config.log.
If the bug report is good than we will do our best to fix it. If the bug report is poor we will
not do much about it, probably just ask for a better report.
If you think something in this manual is unclear, or downright incorrect, or if the language
needs to be improved, please send us a note to the same address.

3.8

Some hints for developers

3.8.1

Organization of the code

The components of the program are located in three main directories. Directory ./src/ contains
the *.cc source files, ./include/ the corresponding .h files and ./input files/ the user
23

defined input files. The directory ./doc/ contains the documentation of the code in LATEX
and HTML forms. Table 3.2 gives a hint on the structure of the code.
././

Makefile
DOXYFILE
src/
include/input files/
results/
doc/

*.cc
*.h
*.dat
nrg user guide.tex
html/
tex/

configuration for on-line documention
source code (source files)
source code (header files)
directory for input files
output directory of the code
documentation of the codeon-line documention
on-line documention

Table 3.2: Directory structure of the code.
The other files in the main directory are configurations files needed for properly installing the
code.

3.8.2

On-line documentation of the code

The documentation of the source code is managed with the use of the doxygen a free software
under Linux. At present version 1.4.6 is used. The configuration of doxygene is set in file
Doxyfile and the program is called from shell as doxywizard. Description of the classes, variables
and functions is included directly in the source code (in the *.cc and *.h files), therefore, allcomments about new variables, functions, classes etc. are documented right away. The on-line
documentation file ./doc/html/....html can be opened with any WEB-browser. The DMNRG reference manual in pdf form (refman.pdf) can be generated by the command: latex./doc/latex/refman.tex.
All the information about the functions (classes) i.e., what it does, how to call, description of
the input and output variables must be located in the main body of the functions, otherwise,
this information will not be accessible in the output files generated by doxygene! The hierarchyof the classes, dependencies and calling functions are generated automatically by doxygen.

3.8.3

Version management

Version management of the code is carried out with the use of Subversion system(svn) afree software under Linux. This program is similar to Control Version System (cvs) but
has several additional features. The WEB page of the svn interface can be accessed by themembers of the flexible DM-NRG developers group only.
The management of this manual – in LATEX form–, however, has to be done manually. There areno scripts yet, that would automatically update the manual once the code has been changed.
However, since all directories and files – including this manual as well – are checked in to the
svn it is straightforward to update the manual and all changes of it are documented by the svn
itself.

24

Chapter 4

Using the DM-NRG code
4.1

Main features of the code

If the installation was successful then, to run the code for your problem of interest, you firstneed an input file, input.dat. The structure and content of this file shall be described indetail later in this chapter, and a few examples are also provided with the code in the library,./input files. Furthermore, we provided some Mathematica files in the mathematica files
folder, which were used to generate these input files. For the generation of the irreducible
matrix elements it is sometimes important to know, what is the definition of the ClebschGordan coefficients used by the code. For this purpose, a utility called cgc is provided in the
package that computes the Clebsch-Gordan coefficients.
When running the code, all the results are written into the folder ./results. The content ofthese files shall also be described later in this chapter. Finally, these results can further beanalyzed with the utilities sfb – real axis spectral function generator and es – energy spectrum
generator, also provided with the code, and shall be described later.
Main features included in version 1.0.0:• Zero and finite temperature calculations of the spectral function of not necessarily thesame operators.
• Static averages.
• NRG type calculations or DM-NRG computations (using the complete basis set).
• The code uses a single input file that can be easily modified. Examples are available forAnderson as well as Kondo Hamiltonians in the presence of U (1) and SU (2) symmetries
for spin as well as charge symmetries.
• Flat density of states as well as arbitrary energy-dependent density of states can be considered. In the latter case arbitrary precision routines are used to generate the parametersof the Wilson chain.
• Utility for broadening the spectral function is provided to compute the imaginary and thereal parts of the retarded Green’s function
D
 E
†
GR
(t)
=
−i
Θ(t) ,
(4.1)
A(t),
B
(0)
†
A,B
±
with the operators A and B not necessarily identical.
25• Utility for computing the on-site energies as well as the hoppings for a non-uniform densityof states is provided.
• Utility to generate/analyze the finite size spectrum is provided.
Utilities provided in the package:• he - generates the hopping amplitudes along the chain and the on-site energies when an
energy-dependent density of states is used. (This utility provides input files that are readdirectly in the fnrg binary code).
• sfb - does the broadening of the spectral function and computes the real part of theGreen’s function by a Hilbert transform. This utility uses the outputs of the fnrg to
generate these data.
• notebook examples - there are some extra Mathematica 6.0 notebook examples where
the matrices for Hamiltonians as well as for some hopping and spectral operators arecomputed. These notebooks can easily be extended to other impurity models.
• es - generates the energy spectrum. This utility uses the output of the fnrg to generate
the data.
• cgc - an interactive utility, which generates the Clebsch-Gordan coefficients used by the
code.

4.2

Running the code

In this section we shall give some information on how exactly to run the code. At this point wesuppose that the code was successfully compiled and linked properly with the required libraries.
If this is the case there should be five binary codes in the main folder, fnrg, sfb, he, cgc, andes. Depending on the type of calculation that is performed, the following steps need to befollowed:

4.2.1

Flat density of states

Most NRG calculations are done for a flat density of states, ̺(ω) = 21 , with ω ∈ [−1, 1]. In thiscase there is no need for the he routine, since the hoppings along the chain are computed on
the fly. In this case the on-site energies are all zero as a consequence of electron-hole symmetry.
For this type of calculation, first the hoppings on site energies flag in the <SECTION-FLAGS>
needs to be set to ’OFF’ in the file input.dat. Then the fnrg must be run without any
arguments in the command line. The name of the input file must always be input.dat, since
the code looks for this particular file to read in the data.
Outputs are saved in a folder in results/. Results for the spectral function can be generated
by running the sfb binary without any argument after exiting normally from the run. Notethat sfb reads the same input file, input.dat, so do not modify it between the two runs
except for the part <SECTION-SPECTRAL FUNCTION BROADENING>. Depending on
the flags set in <SECTION-SPECTRAL FUNCTION BROADENING>, spectral functions,Green’s functions and static averages are computed. Note that the real part of a Green’s
function can only be computed if the corresponding spectral function has previously been
computed.
26

4.2.2

Energy dependent density of states

The user can also provide an arbitrary density of states (DOS) for the code. In the present
version of the code, the DOS is read from a file, dos.dat. This file must be in the folder,./dos mapping, and it must be provided by the user. This file must contain at least twocolumns, one of them being the energy, and the others being the corresponding density ofstates. There should be as many density of states columns as hopping operators are defined in
the input file.
There are some constraints that the DOS must satisfy• The DOS needs to be rescaled such that all energies are in the interval, E ∈ [−1, 1]. In
other words, the energy unit of the NRG calculation (=1) must be larger than the rangeof the density of states.
• The DOS needs to be normalized to 1.
• The energies must be ordered.
• Each density of states column corresponds to a given hopping operator in the input.dat
file. The columns must be ordered accordingly to the order of the hopping operators in
the input.dat file.
There is no further requirement for the mesh, both uniform and non-uniform meshes can beprovided.
When the NRG calculation is performed with an energy-dependent DOS, the hoppings andthe on-site energies must be computed first. This must be done interactively, using the binary
he. When running he, one needs to specify an iteration number. This can be larger than
the number of iterations specified in the input.dat file. After a successful run of he, thehopping amplitudes and the on-site energies are saved in files having descriptive names such as
hopping couplings *.dat and on site energies *.dat in the folder ./dos mapping/.
Next, the input.dat file needs to be set-up: There, make sure that the hoppings on site energies flagin the <SECTION-FLAGS> is set to ’ON’ for this type of run.
Once the hopping amplitudes and the on-site energies were computed, the same steps must be
followed as in the case of a flat density of states. The code fnrg must be run without arguments.
Note: The routine he is using arbitrary precision libraries for computing the hoppings andthe on-site energies, however due to the accumulation of the errors during the calculations
sometimes the he utility may give unexpected results for the hoppings. Therefore, we advise
the user to check the results for the hoppings before pluging them into fnrg for any kind of nrg
calculations.

4.3

Initialization and the input fileIn the initialization step – the zero’th NRG iteration step – one has to define the Hamiltonian.
Therefore, for the iterative process, among others, the user has to provide the following input
parameters:
1. One needs to specify the number (symmetry no) and type of symmetries.
2. One needs to specify the number (block state number) of initial block states, i.e., multipletsthat span the Hilbert space for n = 0 (we denote this number in this manual E
by K). One
also needs to define the quantum numbers Qu of the block states, u, Qu , Qzu

0.

3. Block Hamiltonian matrix - The matrix elements of the K × K Hamiltonian matrix H0 ,defined on the K block states.
27

4. block hopping operators - The reduced matrix elements of all hopping operators between
the block states. These are K × K matrices.
5. spectral operators - The reduced matrix elements of the operators, whose spectral functions
we compute. These are also K × K matrices.
6. local state number - The number of multiplets that span the local Hilbert space (we label
this number L).
E
7. local states - Quantum numbers q µ of the local states, µ, q µ , q zµ .

8. local hopping operators - The reduced matrix elements of the hopping operators. Theseare L × L matrices.
This list is rather incomplete, and explains only briefly the meaning of the most basic variables. In addition, many other, less relevant parameters must be specified for the NRG process(broadening, maximum iteration number, etc.). All these input parameters are read from file,
input.dat, which must be constructed for every model separately. Some examples of such
input files can be found in the folder, ./input files, and also in the appendix of this manual.
A detailed description of these parameters shall be provided in forthcoming subsections.

Structure of the input file
The input file is organized along sections, and contains a lot of comments. As a general rule,
any line that starts with the symbol ’#’ is treated as a comment line in these input files andis ignored by the code. You can thus add additional comment lines for yourself using thisconvention.
The input file has the following overall structure (that we shall describe in detail in subsequent
subsections):
1. <SECTION-PARAMETERS>
2. <SECTION-FLAGS>
3. <SECTION-SYMMETRIES>
4. <SECTION-BLOCK STATES>
5. <SECTION-LOCAL STATES>
6. <SECTION-LOCAL STATES SIGNS>
7. <SECTION-BLOCK HAMILTONIAN>(a) <BLOCK HAMILTONIAN TERM>
i. <BLOCK HAMILTONIAN NAME>
ii. <BLOCK HAMILTONIAN COUPLING>
iii. <BLOCK HAMILTONIAN REPRESENTATION INDEX>
iv. <BLOCK HAMILTONIAN MATRIX>
8. <SECTION-HOPPING OPERATORS>(a) <HOPPING OPERATOR>
i. <HOPPING OPERATOR NAME>
ii. <HOPPING OPERATOR REPRESENTATION INDEX>
iii. <HOPPING OPERATOR SIGN>
28

iv. <HOPPING OPERATOR MATRIX>
9. <SECTION-SPECTRAL OPERATORS>(a) <SPECTRAL OPERATOR>
i. <SPECTRAL OPERATOR NAME>
ii. <SPECTRAL OPERATOR REPRESENTATION INDEX>
iii. <SPECTRAL OPERATOR SIGN>
iv. <SPECTRAL OPERATOR MATRIX>
10. <SECTION-STATIC OPERATORS>(a) <STATIC OPERATOR>
i. <STATIC OPERATOR NAME>
ii. <STATIC OPERATOR REPRESENTATION INDEX>
iii. <STATIC OPERATOR SIGN>
iv. <STATIC OPERATOR MATRIX>
11. <SECTION-SPECTRAL FUNCTION>
12. <SECTION-SPECTRAL FUNCTION BROADENING>
13. <SECTION-LOCAL HOPPING OPERATORS>(a) <LOCAL HOPPING OPERATOR>
i. <LOCAL HOPPING OPERATOR NAME>
ii. <LOCAL HOPPING OPERATOR REPRESENTATION INDEX>
iii. <LOCAL HOPPING OPERATOR SIGN>
iv. <LOCAL HOPPING OPERATOR MATRIX>
v. <LOCAL ON SITE ENERGY MATRIX>

4.3.1<SECTION-PARAMETERS>In this section, the main parameters of the given model and also those of the DM-NRG method
are set. The order of the parameters does not matter as long as none of them is missing.
Model parameters• model: The name of the model under investigation. It can be any string. This name will be
automatically included in the name of the folder, where output data are saved.
• lambda: The value for the logarithmic discretization parameter, Λ. It is usually fixed between
2 and 3.
• symmetry no: The number of the symmetries used. The present version of the code handles
U (1), SU (2), Z(2), charge SU (2) symmetries. The properties of the symmetries are set in<SECTION-SYMMETRIES> by defining their names and setting upper and lower bounds ontheir representation indices. Later versions of the code will include SU (3) symmetry and crystal
field symmetries too.
• coupling no: The number of independent coupling constants. In general, the local Hamiltonian H0 can contain several couplings. The on-site part of the Anderson model, e.g., can be
written as
X
H0 = U (n − 1)2 + (ǫd + U/2) n + V(d†σ f0,σ + h.c.) ,σ29

i.e. it contains three different couplings, U , ǫd and V . These couple to the operators,(n − 1)2, n,and (d†σ f0,σ +h.c.). Thus the local Hamiltonian, H0 , can be written as a sum of coupling no = 3
independent terms,
coupling
X noGδ Hδ ,
(4.2)
H0 =
δ=1

where Hδ ’s denote the independent local terms and Gδ stand for the corresponding couplings.
The matrix elements of Hδ are specified in <SECTION-BLOCK HAMILTONIAN>.
Important note: The number of <BLOCK HAMILTONIAN TERM> items must agree withthe ’coupling no’ set in this section.
• spectral operator no: The number of the operators for which one wants to computethe spectral functions. The detailed parametrization of these operators is carried out in<SECTION-SPECTRAL OPERATORS>. The number of <SPECTRAL OPERATOR> items
must agree with the ’spectral operator no’ set in this section.
• static operator no: The number of operators for which static averages are computed. The
detailed parametrization of these operators is carried out in <SECTION-STATIC OPERATORS>.
The number of <STATIC OPERATOR> items there must agree with the ’static operator no’
set in this section.
• hopping operator no: The number of independent hopping operators. The detailed parametrization of these operators is carried out in <SECTION-HOPPING OPERATORS> and <SECTIONLOCAL HOPPING OPERATORS>.
• block state no: The number of multiplets that span the Hilbert space for iteration n = 0.
The detailed parametrization of the states is carried out in <SECTION-BLOCK STATES>
by setting the representation indices of all the states corresponding to all the symmetries usedduring the calculation.
• local state no: The number of local multiplets, added in each iteration. In the present
version of the code this is the same for all sites, 0 < i ≤ N . The detailed parametrization of
local states is carried out in <SECTION-LOCAL STATES>.
• local coupling no: The number of independent terms in the on-site Hamiltonian operators,
Hn . In the presence of electron-hole symmetry of the conduction band and without any othercorrelations (superconducting like) the on-site terms are all zero. If we have only on-site energies,ξn , then local coupling no= 0.
• spectral function no: The number of the spectral functions that should be calculated.
One can compute several spectral functions within the same run, at the expense of using
more memory. The detailed parametrization of the operators is carried out in <SECTIONSPECTRAL OPERATORS>.
DM-NRG parameters• max state no: The upper limit for the number of kept block multiplets.
• iteration no: The maximum number of the NRG iteration steps, N . The program terminates
once this value is reached.
• interval no: The energy grid used in each NRG iteration step when spectral functions arecalculated. Note that in DM-NRG a very large number of delta peaks is generated in every
iteration. Therefore, these delta peaks are stored on a finite mesh. In the present version of
the code, the difference between the maximum and minimum of the energy spectrum at each
iteration is divided by interval no to obtain the size of the mesh, dω. Typically interval no ∼1000.
• degeneracy threshhold: If two states are such that the difference between their energies
is less than degeneracy threshhold these states are considered to be degenerate.
When discarding states, those lowest in energy are kept while the higher energy states are
discarded, such that at the next iteration no more than max state no are kept. Degenerate
30

states are block discarded.
• temperature: The temperature used for the calculation. When the temperature is fixed
to zero, temperature = 0, the code runs by using the number of iterations provided by theuser. In case of a finite temperature, temperature > 0, the variable iteration no is determined
and modified during the run to correspond to the temperature specified, using the well-known
relationship, Tn ≃ (1 + Λ1 ) Λ−n/2 .

4.3.2<SECTION-FLAGS>In this section various flags controlling the DM-NRG method are set.
• dmnrg flag: This flag controls whether to perform the backward sweep and use the DM-NRG
method or not. Allowed tags are ’ON’ or ’OFF’. In case of ’ON’, the DM-NRG backward sweep
is carried out and the complete set of eigenstates is used to calculate the spectral functions. In
case of ’OFF’ the usual NRG procedure is performed only.
• text swap files flag: This flag controls whether to generate various temporary scratch files
as readable text files in addition to the binary files. Allowed tags are ’ON’ or ’OFF’.
• hoppings on site energy flag: This flag is important when energy-dependent density ofstates are considered. This flag controls whether the hoppings and the on-site energies are readfrom a file or not. If set to ’ON’, the hoppings as well as the on-site energies must be computed
using the ’he’ utility that comes with the code. This utility generates some data files in the./dos mapping folder, which are next read by the fnrg binary. When a flat band is considered,this flag can be set to ’OFF’ and then the hoppings are generated on the fly.

4.3.3

<SECTION-SYMMETRIES>In this section the types of symmetries used can be set. As many symmetries must be specified
as were set by the ’symmetry no’ variable. Each line must contain the name of a symmetry,
and a minimum and maximum value for the corresponding representation index. Those states
for which these bounds are exceeded are simply discarded, and to avoid such uncontrolled
truncation, a large value (≃ 20) should be used for most symmetries. The possible symmetry
types and the corresponding tags that the present version of the code can handle are as follows:’U (1)’, ’SU (2)’, ’Z(2)’, ’charge SU (2)’. Other symmetries will also be included in later versionsof the code.

4.3.4<SECTION-BLOCK STATES>In this section the quantum numbers of the block states are set. A matrix of dimensions
block state no×symmetry no with integer entries enumerating the representation indices of theblock states must be provided. The number of rows has to be equal to the value ’block state no’
set in <SECTION-PARAMETERS>. Integers in a row u (u = 1, . . . , block state no) correspond to the representation indices (quantum numbers) of the state u. Thus each row contains
symmetry no integer entries. Note that it is necessary to provide the double of the usual quantum numbers in case of SU (2) symmetries in order to be able to work with integer numbers.
It is important that the states must be ordered by symmetries: The ordering is done withrespect to the first symmetry, then with respect to the second symmetry used and so on. While
reading the input file, the code checks if the quantum numbers are in ascending order, and if
not, then an error message is generated and the run stops.

4.3.5<SECTION-LOCAL STATES>In this section the quantum numbers are set. An integer type matrix of dimensions local state no×symmetry no must be provided. The number of rows must be equal to the value
31

’local state no’ set in <SECTION-PARAMETERS>. Otherwise, the rules are the same as for<SECTION-BLOCK STATES>.

4.3.6<SECTION-LOCAL STATES SIGNS>

In this section the signs of the local states must be set, and a vector of length ’local state no’,containing ±1-s must be provided. The sign of the local state is -1 / +1 if it contains an odd/ even number of fermions.

4.3.7<SECTION-BLOCK HAMILTONIAN>

This section contains the definition of the interaction Hamiltonians, Hδ .
• <BLOCK HAMILTONIAN TERM> tag denotes the beginning of the definition ofa new interaction Hamiltonian term.
• <BLOCK HAMILTONIAN NAME> is a character string of maximum length char[256]to define the name of the Hδ term. It could be any string.
• <BLOCK HAMILTONIAN COUPLING> a real*8 number, which sets the coupling strength, (Gδ ), corresponding to Hδ . In this line you also specify the name (notation) of the coupling. By writing a line ”J = 0.736”, e.g., you teach the code that thiscoupling is called J, and its value is 0.736.
• <BLOCK HAMILTONIAN REPRESENTATION INDEX>
the quantum numbers of the Hamiltonian term, an integer vector of length ’symmetry no’.
Under any circumstances the quantum numbers must be 0 for all the symmetries, otherwise the problem is badly posed. The representation index of the Hamiltonian may looktherefore redundant at first, but to have a uniform treatment of all operators in the input
file, we decided to keep this line.
• <BLOCK HAMILTONIAN MATRIX> This is a matrix of dimension block state no× block state no. This matrix contains the matrix elements of the Hamiltonian Hδ .
between the states (multiplets) enumerated in <SECTION-BLOCK STATES>. Some
Mathematica 6.0 notebooks are provided in directory mathematica files which can beused to generate these matrices in a routine way.
The above series of entries must be repeated for every local Hamiltonian term, Hδ . The numberof these structures must be equal to coupling no defined in <SECTION-PARAMETERS>.

4.3.8<SECTION-HOPPING OPERATORS>This section provides information on the hopping operators. As a rule, we provide information
only on creation operators. No information on annihilation operators is needed.
• <HOPPING OPERATOR> tag denotes the beginning of the definition of a new
hopping operator.
• <HOPPING OPERATOR NAME> is a character string with maximum length char[256]to define the name of the corresponding hopping operator. (Could be something like“spin up hopping”, e.g.)• <HOPPING OPERATOR REPRESENTATION INDEX> the quantum numbersof the hopping operator, an integer vector of dimension 1 × ’symmetry no’ must be given.
32• <HOPPING OPERATOR SIGN> The sign is -1 for a fermionic operator and 1 fora bosonic operator.
• <HOPPING OPERATOR MATRIX> Contains the reduced matrix elements of the
hopping operator (∼ f0† ) between the block states specified in <SECTION-BLOCK STATES>.
It is a block state no × block state no matrix. We provided a few M athematica notebooks which you can use to generate these matrix elements in a simple way.
The above series of entries must be repeated as many times as many hopping operators one
has. The number of these structures must thus agree with the ’hopping operator no’ set in<SECTION-PARAMETERS>.

4.3.9<SECTION-SPECTRAL OPERATORS>This section provides information on the spectral operators, i.e., operators, whose spectral
functions are computed. As a rule, we provide information only on the ’creation operators’,
and no information on annihilation operators is given (see Section 2.2 for a precise definition).
The following structure must be repeated for every spectral operator. Clearly, the number ofstructures must be equal to the parameter spectral operator no set in section <SECTIONPARAMETERS>.
• <SPECTRAL OPERATOR> tag denotes the beginning of the definition of a newspectral operator.
• <SPECTRAL OPERATOR NAME> is a character string with maximum lengthchar[256] to define the name of the corresponding spectral operator. (E.g., “spin operator”)• <SPECTRAL OPERATOR REPRESENTATION INDEX> the quantum numbers of the spectral operator. An integer like vector of dimension 1 × symmetry no is
expected.
• <SPECTRAL OPERATOR SIGN> The sign is -1 for a fermionic-like operator and
1 for a bosonic operator.
• <SPECTRAL OPERATOR MATRIX> A matrix of dimensions block state no ×block state no is expected, containing the reduced matrix elements of the spectral operator
between the block states specified in <SECTION-BLOCK STATES>. We included in thispackage a few examples of M athematica notebooks, which can be used to generate thesematrices relatively easily.

4.3.10<SECTION-STATIC OPERATORS>This section provides information on the operators for which static averages are evaluated. The
following structure must be repeated for every static operator.
• <STATIC OPERATOR> tag denotes the beginning of the definition of a static operator.
• <STATIC OPERATOR NAME> is a character string with maximum length char[256],
defines the name of the static operator being defined.
• <STATIC OPERATOR REPRESENTATION INDEX> quantum numbers of thestatic operator, an integer type vector of dimension 1 × symmetry no is expected.
33

• <STATIC OPERATOR SIGN> This sign is -1 for a fermionic and 1 for a bosonicoperator.
• <STATIC OPERATOR MATRIX> Contains the reduced matrix elements of thestatic operator between the block states specified in <SECTION-BLOCK STATES>. It
is a block state no × block state no matrix. For the calculation of these matrix elements,
one can use the mathematica codes provided with the package.
The number of structures shown above must agree with the variable ’static operator no’ set in<SECTION-PARAMETERS>.

4.3.11<SECTION-LOCAL HAMILTONIAN>• <LOCAL HAMILTONIAN TERM> tag denotes the beginning of the definition ofa local Hamiltonian, i.e. a term that occurs in Hi ,Hi =local coupling
no
X
α=1

gα,i Hi,α ,with Hi,α = Hα ({C i,λ }). In other words, the local Hamiltonian Hi,α depends on the site
index i only through the creation operators acting at site i, otherwise its structure is fixed.
A trivial example of such operators are the energies, ξλ ↔ g, C i,λ C †i,λ ↔ Hi,α . In thepresent version of the code, we only handle uniform local couplings, gα,i = gα , however,later versions shall also be able to treat a user-specified coupling set, gα,i . Note that the onsite energy terms are specified under the section <SECTION-LOCAL HOPPING OPERATORS>,
not here. Here you can define additional local terms in the Hamiltonian, e.g., superconducting pairing terms.
• <LOCAL HAMILTONIAN NAME> is a character string with maximum lengthchar[256], defines the name of the Hi term.
• <LOCAL HAMILTONIAN COUPLING> is a real*8 number to set the couplingstrength, gα , corresponding to the local term, Hi,α . In the present version of the code thecoupling strength does not change along the Wilson chain.
• <LOCAL HAMILTONIAN REPRESENTATION INDEX> the quantum numbers of the Hamiltonian term, an integer type vector with dimension 1 × ’symmetry no’.
Similar to the Hamiltonian, these quantum numbers must be all 0, otherwise the Hamiltonian is not invariant under the symmetry assumed.
• <LOCAL HAMILTONIAN MATRIX> It is a local state no × local state no matrix, defined on the basis given in <SECTION-LOCAL STATES>. This contains the(reduced) matrix elements of the coupling Hamiltonians, Hα .

4.3.12<SECTION-LOCAL HOPPING OPERATORS>

Here we specify matrix elements of the local hopping operators between the local states enumerated in <SECTION-LOCAL STATES>. As a rule, we provide the information only for thecreation operators.
• <LOCAL HOPPING OPERATOR> denotes the beginning of the definition of a
hopping operator.
• <LOCAL HOPPING OPERATOR NAME> is a character string with maximum
length char[256]. Defines the name of the local hopping operator.
34

• <LOCAL HOPPING OPERATOR REPRESENTATION INDEX> the quantum numbers of the local hopping operator, an integer type vector with dimension 1 ×’symmetry no’ must be given.
• <LOCAL HOPPING OPERATOR SIGN> This sign is -1 for a fermionic and 1 fora bosonic operator.
• <LOCAL HOPPING OPERATOR MATRIX> Contains the reduced matrix elements of the static operator between the local states specified in <SECTION-LOCAL STATES>.
It is a local state no × local state no matrix. For the calculation of these matrix elements,
one can use the mathematica codes provided with the package.
• <LOCAL ON SITE ENERGY MATRIX> Contains the reduced matrix elementsof the on-site energy term between the local states specified in <SECTION-LOCAL STATES>.
It is a local state no × local state no matrix. For the calculation of these matrix elements,
one can use the mathematica codes provided with the package.
The number of structures shown above must agree with the variable ’hopping operator no’ set
in the <SECTION-PARAMETERS>.

4.3.13<SECTION-SPECTRAL FUNCTION>

In this section one specifies the pairs of spectral operators, whose spectral function is computed.
This section must contain ’spectral function no’ rows, each of them corresponding to a given
pair of operators.
To compute the spectral function of the Green’s function of the form
D
 E
A(t), B † (0) ξ Θ(t),GR(4.3)
A,B † (t) = −ione must add a row of the form

{name of A† ; name of B † } ,where ’name of A† ’ and ’name of B † ’ are the names of the ’creation operators’ set up under<SECTION-SPECTRAL OPERATORS>.

4.3.14<SECTION-SPECTRAL FUNCTION BROADENING>

Here you specify the parameters of the broadening, used to evaluate the spectral functions.
Note that the code fnrg does not generate the spectral functions, you need to run the utility
sfb for that immediately after the run. This part of the input file is thus used by sfb. Notethat one can improve/change the broadening parameters after running the code too, and thenrun sfb without running fnrg again. For details on the broadening, see Section 2.2).
The following parameters must be set in this section:• broadening method: the method that is used for the broadening of the spectral function along the real axis. Two methods are available: LOG GAUSS and INTERPOLATIVE LOG GAUSS.
• broadening parameter: The value of the broadening parameter (half-width of the log
Gauss function in both the log-Gauss and interpolative log-Gauss type of broadening).
This parameter is used no matter what method is used.
• broadening energy minim: The minimum energy for which the broadening is performed (broadening energy minim < |ω| < broadening energy maxim). This parameter is used for both methods. This quantity must be positive.
35

• broadening energy maxim: The maximum energy for which the broadening is performed (broadening energy minim < |ω| < broadening energy maxim). This quantity
must be positive.
• broadening grid mesh: The number of points along the real axis, for which the broadened spectral function is computed. The binary sfb generates a logarithmic mesh between[-broadening energy maxim, -broadening energy minim] and [broadening energy minim,
broadening energy maxim] for which the spectral function is computed.
• quantum temperature: This parameter is used only with the INTERPOLATIVE LOG GAUSS
method. It sets the scale below which a linear interpolation is used rather than a logarithmic interpolation.
• spectral function nrg even: This is a flag that is relevant for NRG runs only. Allowed
tags are ’YES’ or ’NO’: In case of ’YES’ the spectral function is computed from the even
NRG iterations, and a corresponding data file is generated. If set to ’NO’, the code skips
this calculation.
• spectral function nrg odd: This is a flag that is relevant for NRG runs only. Allowed
tags are ’YES’ or ’NO’: In case of ’YES’ the spectral function is also computed from odd
NRG iterations. If set to ’NO’, the code skips this calculation.
• spectral function dmnrg: This is a flag parameter that controls whether the spectral
function from the full NRG spectrum using the density matrix formalism is computed.
Allowed tags are ’YES’ or ’NO’. In case of ’YES’ the DM-NRG calculation is also performed.
• static average dmnrg: Flag parameter that controls whether the static average of theform hA† Bi is also computed by sfb for the two spectral operators, in addition to thecalculation of the spectral functions (see the notations in the previous section). Allowed
tags are ’YES’ or ’NO’.
• green function nrg even: This flag controls whether the real part of the Green’s function is computed from the even part of the spectral function generated from the nrg data,by performing a Hilbert transform. Allowed tags are ’YES’ or ’NO’. The real part of theGreen’s function can be computed only after the calculation of the spectral function, i.e.
if the spectral function nrg even flag is set to YES.
• green function nrg odd: This flag controls whether the real part of the Green’s function is computed from the odd part of the spectral function generated from the nrg data,by performing a Hilbert transform. Allowed tags are ’YES’ or ’NO’. The real part of theGreen’s function can be computed only after the calculation of the spectral function, i.e.
if the spectral function nrg odd flag is set to YES.
• green function dmnrg: This flag controls whether the real part of the Green’s function
is computed from the dmnrg spectral function by performing a Hilbert transform. Allowed
tags are ’YES’ or ’NO’. The real part of the Green’s function can be computed only after
the calculation of the spectral function, i.e. if the spectral function dmnrg flag is set
to YES.
• green function grid mesh: The number of points on a log-mesh used for doing theHilbert transform. Once the spectral function for a combination of spectral operators is computed, the real part of corresponding Green’s function is generated by performing a Hilbert transformation. To compute the Hilbert transform, a logarithmic
mesh is generated. The energy limits for the real/imaginary part of the Green’s functions are the same as those for the spectral function , i.e. [-broadening energy maxim,
36

-broadening energy minim] and [broadening energy minim, broadening energy maxim].
The energy values for which the Green’s function is computed correspond to those forthe spectral function. Cubic splines are used for the interpolation of the spectral function. Therefore, for too small green function grid mesh values sometimes oscillations areobserved in the energy dependence of the real part of the Green’s function. In this case,green f unction grid mesh must be increased until the oscillations vanish. Increasing thisparameter has only a small impact on computing time.

37

4.4

Outputs of the DM-NRG code

4.4.1

Output directory structure

The output of the code has the following directory structure:./results/./results/

all results of the code
all data for a given run

automatic directory name/(code generated)
log file.dat
input.dat
spectral function dmnrg *spectral function even *spectral function odd *green function dmnrg *green function even *green function odd *static average dmnrg*
energy levels even*energy levels odd*
data/

spectral operators/spectral operator *.dat
static operators/static operator *.dat
sectors/sectors *.dat
transformations/transformation *.dat
states/states *.dat
delta peaks/delta peaks *.dat
spectral functions/spectral function *.dat
mapping/dos.dat
mapping/hopping couplings*.dat
mapping/on site energies*.dat
mapping/density of states*.datm star.dat

Table 4.1: Output directory structure of the code. Directory automatic directory name/ isgenerated automatically from the parameters in the input file. The “*” in the file names also
stands for a list of the parameters (names of the couplings etc.) generated automatically.
All outputs of the code, which are saved on disk, are located in the directory results. In
this directory the program generates automatically a directory corresponding to the parameters of the given run. The name of this directory, automatic directory name/, contains thevalues of the following input variables: model, lambda, max state no, iteration no, temperature,symmetry no, symmetry types, coupling no, and a file log file.dat containing the log/error
messages during the run.
The automatic directory name/, folder has the following structure:• input.dat is a copy of the original input file.
• Various *.dat files contain the final results for the spectral functions, static averages,
green’s functions and the energy levels. The names of these files are descriptive, andcontain some additional tags such as, e.g., the names of quantum operators for which thespectral functions have been computed.
After running fnrg, only a directory data containing the raw data is generated. The spectral
functions, Green’s functions, static averages files are generated only later by running the sfb
utility. The energy level files are also generated separately, by running the es utility.
• Directory data contains the data files generated by fnrg. Since there are hundreds of files,
this directory is divided into subdirectories with descriptive names and containing the following
38

files:mapping/
sectors/sectors *.dat
states/states *.dat
delta peaks/delta peaks dmnrg *delta peaks/delta peaks nrg *spectral operators/spectral operator *static operators/static average dmnrg *static operators/static average nrg even *static operators/static average nrg odd *transformations/transformation *.dat
spectral functions/spectral function *.datm star.dat Here again, *’s stand for operator names and code-generated tags such as
the iteration number etc.

4.4.2

Detailed description of the output files in the folder./results/automatic directory name/data

This folder contains the original data generated/used by the code fnrg and the utilities he andsfb. It contains the following files/folders:• mapping/: This folder contains information on the hopping parameters and on site energies,
in case a user-defined density of states has been used. This directory is a copy of the directory
dos mapping, located in the main directory. It contains the following files:• mapping/hopping couplings*.dat: Contains the data for the hopping couplings along
the chain as computed by the he (hopping energies) utility. Each line contains information for
one iteration, the first column being the iteration number and subsequent columns representing
the values of the hopping amplitudes along the Wilson chain. Each column corresponds to a
hopping operator, so the number of columns has to agree with the ’hopping operator no’. The
tag “self-consistency” in the file name has no relevance here, it was included for later purposes.
• mapping/on site energies*.dat contains the data for the on-site energies along the chain
as computed by the he (hopping energies) utility. Each line contains information for each
iteration, the first column being the iteration number and the subsequent columns being thevalues of the on-site energies. Each column corresponds to a hopping operator, so the numberof columns has to agree with the ’hopping operator no’. The tag “self-consistency”in the filename has no relevance here, it was included for later purposes.
• sectors/sectors (iteration no).dat The code is organized using so-called sectors. A given
sector is characterized by a list of representation labels, specifying the symmetries of the statesof that sector. Similarly, operators are also organized using the notion of sectors. The files
sectors (iteration no).dat contain the information on the size and symmetry of sectors in
iteration iteration no.
This file has two separate parts. The first part describes the sectors in the new basis at iteration
iteration no. The first row of the file contains the total number of the sectors in the new basis(new block sector no). The subsequent new block sector no lines contain the information foreach sector (from 0 to new block sector no-1): the representation indices of the sectors, thedimensions of the sectors, the first states of the sectors among the list of all the basis states(multiplets), and the degeneracies of the sectors. Note that the sum of the dimensions multipliedwith the corresponding degeneracy values must be equal to the total number of basis states of
the new block state basis.
39

The second part of the files sectors (iteration no).dat contains the same information on
the states kept from the previous iteration (iteration iteration no-1).
• states/states (iteration no).dat: These files contain the information for the new basis
states at iteration iteration no, and the states kept from the previous iteration (iteration iteration no-1). States know about their parent states, and this information is also displayed inthese files.
The first row of the file contains the total number of the new basis states (new block state no)
The next new block state no lines contain the information on the multiplets after the diagonalization has been done. Each line contains the number (label) of the sector the state belongs
to, the index of the state within this sector, the list of representation indices, the labels of
the parent and local states the state was constructed from, the sign of the state (this is 0 for
block states), whether the state is kept (K) or eliminated (E) after truncation, and finally the
eigen-energy of the state. Note that the sector number cannot be less than zero and it must be
smaller than new block sector no, specified in file sectors (iteration no).dat. Also the variablesector index must be non-negative and less than the dimension of the sector given also in file
sectors (iteration no).dat.
The second part of the file states/states (iteration no).dat gives the same informationon the states kept from the previous iteration.
• transformations/transformation (iteration no).dat: These files contain the information on the matrix On that is used to diagonalize the Hamiltonian, and to transform thehopping operators at site iteration no, and the spectral operators (acting at site 0) to the newbasis asAn → On An On†(4.4)
These files are only generated if the flag text swap files flag in the file input.dat is set ON.
The first row of the file contains the total number of the sectors in the new basis in iteration
iteration no, and this is followed by three terms for each sector: the sector number, the size ofthe sector, and finally the real*8 matrix of size sector size × sector size, transforming the given
sector.
• transformations/transformation binary (iteration no).dat contains the information
for the O matrix performing the diagonalization at iteration iteration no in a binary form.
These files are only generated if the flag binary swap files flag in the file input.dat is set toON.
• spectral operators/spectral operator (operator name) (iteration no).dat contains
the matrix elements of a spectral operator acting at site 0 of the chain (for example d†0 or
Sz ) in a sector-decomposed form, after the diagonalization in iteration step iteration no has
been done. These operators (which are ’creation operators’) are defined and initialized in the
input.dat file, <SECTION-SPECTRAL OPERATORS>. These files are only generated if
the flag text swap files flag in the file input.dat is set ON.
The matrix of a spectral operator can be decomposed to sector no by sector no small blocks.
Depending on the quantum numbers of the operator, several sectors (actually most of them)
contain only zero elements. These are considered as non-active sectors. The first line of the filecontains the total number of the active sectors (those with proper change of quantum numbers).
Then for each active sector the row and column index of the sector is given and is followed by
the row and column dimension of the sector, and the matrix elements of the operator in that
sector.
• spectral operators/spectral operator binary (operator name) (iteration no).dat:
These files contain the same information as the files spectral operator (operator name) *.dat,40

excepting that the data are in a binary form. These files are only generated if the flag binary swap files flag in the file input.dat is set to ON.
• static operators/static operator (operator name) (iteration no).dat: They containthe matrix elements of a static operator acting at site 0 of the chain (for example n0 ) in asector-decomposed form, after the diagonalization in iteration step iteration no has been done.
These operators (which are ’creation operators’) are read from the input.dat file, <SECTIONSTATIC OPERATORS>. These files are only generated if the flag text swap files flag in thefile input.dat is set ON.
Static operators must have a block-diagonal structure by symmetry, therefore only these sectors
are active. These files are organized otherwise in the same way as the files spectral operator *.dat.
These files are only generated if the flag text swap files flag in the file input.dat is set to ON.
Furthermore corresponding binary files are also generated if the flag binary swap files flag inthe file input.dat is set to ON.
• delta peaks/delta peaks dmnrg (operator 1) (operator 2) (iteration no ).dat contain the raw spectral function peaks for the spectral function of operator 1 and operator 2,
generated from iteration iteration no. These files contain two columns: the first one is theenergy and the second one the contribution to the amplitude of the spectral function at thatenergy in iteration iteration no. The calculation is performed at a temperature T specified
in the input file, input.dat. These output files are used by the utility sfb to compute the
DM-NRG spectral functions, combining the data from all iterations.
• delta peaks/delta peaks nrg (operator name 1) (operator name 2) (iteration no).dat contains peaks for the spectral function of operator 1 and operator 2, generated fromiteration iteration no. These files contain two columns, the first one being the renormalizationenergy and the second one the amplitude of the spectral function at the corresponding energy
and iteration iteration no. The calculation is performed at a temperature T specified in the
input file, input.dat. These files are used to compute the NRG spectral function combining
data from either even or odd iterations.
• spectral functions/spectral function dmnrg (operator 1) (operator 2)(broadening par) (iteration) .dat contains two columns, the first one being the energy andthe second one the real axis value of the spectral function corresponding to the two operatorsfrom the filename at the iteration iteration no. Data for negative as well as positive energies are
saved into this file. These particular files are relatively important for analyzing the contributionof each iteration to the total spectral function.
• static operators/static average nrg even (operator name).dat andstatic average nrg odd (operator name).dat: These files contain average of a static operator as a function of temperature, computed by the standard NRG prescription, from theeven/odd iterations, respectively.
• static operators/static average dmnrg (operator name).dat andstatic average dmnrg (operator1) (operator2).dat: These files contain the expectation
values of the static operators and of the pairs of spectral operators specified in the file input.dat,
computed using the DM-NRG procedure at a temperature T , also specified in the file input.dat.
• m star.dat: m star is the value of the iteration index at which truncation of the states begins(new state no > max state no), i.e., for iterations larger then m star the states, transformation
matrices, spectral operators have to be saved. m star remains fixed for the rest of the NRG
steps. This file is set on input and by function get m star() in cut off.cc and read by various
functions when spectral functions are generated later.
41

4.4.3

Analyzing the data

There are several utilities that we provide to analyze the data you generated.
• First of all, spectral functions and their Hilbert transforms can be generated/analyzed by
the utility sfb. The broadening parameters must be set in the file input.dat sitting in
the main folder. You can set there the broadening method, as described in Chapter 4.3.
(See also the last section of Chapter 2.) Depending on the parameters set this utility
generates the following files.
spectral function dmnrg (operator1) (operator2) (broadening).dat
These files contain the broadened spectral function of operators operator1 and operator2. Correspondingly, in case the Hilbert transform flag is ON in the file input.dat,
the real part of the corresponding Green’s function is also generated in file green function dmnrg *.dat.
spectral function even (operator1) (operator2) (broadening).dat andspectral function odd (operator1) (operator2) (broadening).dat:
These files contain the broadened spectral functions of operator1 and operator2, ascomputed by the traditional NRG method. Data are collected from even and odd
iterations respectively. If the Hilbert transform flag is ON in the file input.dat,
then the real parts of the corresponding Green’s functions are also generated andstored in the files green function even *.dat and green function odd *.dat.
• The energy spectrum (finite size spectrum) can be analyzed by running the utility es.
While running this utility you can specify quantum number filtering as well as the number of states that you want to keep track of. This is a useful utility to obtain the fixed
point spectrum of a model, identify the corresponding conformal field theory, the dimensions of the irrelevant operators, or to extract Fermi liquid parameters. The results
are stored separately for even and odd iterations in the files energy levels even* andenergy levels odd*.

42

Acknowledgments
We are deeply indebted to a number of people for the useful discussions and comments that
helped us to create this flexible code. First of all, we would like to thank Laszlo Borda, withwhom we created the first NRG code in Hungary many years ago, and whose experience, advices
and constant support were crucial to write this flexible code. We benefited a lot from discussions
with other NRG experts too, including Mikito Koga, Ralf Bulla, Frithjof Anders, and Walter
Hofstetter. We are also indebted to Jan von Delft and Andreas Weichselbaum for many longand constructive discussions. Finally, we would like to thank Laszlo Udvardi for taking care of
our computer cluster at the TU Budapest (BUTE), where most of the work has been carriedout.
This work was supported by the Hungarian Research Fund (OTKA) under Grant Nos. NF61726,
T046303, K68340, and K73361. We are also grateful to the Institute of Mathematics of the
BUTE for providing access to their computer cluster (supported by OTKA under grant No.
63066). C.P.M. was partially supported by the Romanian Grant No. CNCSIS 780/2007. I.W.
acknowledges support from the Foundation for Polish Science.

43

44

Bibliography[1] A. C. Hewson, The Kondo Problem to Heavy Fermions (Cambridge University Press, Cambridge, 1993).
[2] D. L. Cox and A. Zawadowski, Adv. Phys. 47, 599 (1998).
[3] Robert Peters, Thomas Pruschke and Frithjof B. Anders, Phys. Rev. B 74, 245114 (2006).
[4] A. Weichselbaum and J. von Delft, Phys. Rev. Lett. 99, 076402 (2007).
[5] A. I. Toth, C. P. Moca, O. Legeza, G. Zarand, arXiv:0802.4332.
[6] H. Saberi, A. Weichselbaum and J. von Delft, arXiv:0804.0193.
[7] S. R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev. B 48, 10345 (1993); W. Hofstetter, Phys. Rev. Lett. 85, 1508 (2000).
[8] Frithjof B. Anders and Avraham Schiller Phys. Rev. Lett. 95, 196801 (2005).
[9] F. B. Anders, arXiv:0802.0371.
[10] U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
[11] K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
[12] Krishna-murthy, H. R., J. W. Wilkins, and K. G. Wilson, Phys. Rev. B 21, 1003 (1980);
Krishna-murthy, H. R., J. W. Wilkins, and K. G. Wilson, Phys. Rev. B 21, 1044 (1980).
[13] T. A. Costi and A. C. Hewson, Phil. Mag. B 65, 1165 (1992); T. A. Costi, A. C. Hewson
and V. Zlatić, J. Phys.: Condens. Matter 6, 2519 (1994).
[14] R. Bulla and A. C. Hewson, Z. Phys. B 104, 333 (1997).
[15] T. Hecht, A. Weichselbaum, J. von Delft and R. Bulla J. Phys.: Condens. Matter 20,
275213 (2008).
[16] R. Bulla, T. A. Costi and D. Vollhardt, Phys. Rev. B 64, 045103 (2001).

45

46

Appendix A

Input file for the Kondo modelwith Ucharge(1) × Uspin(1)symmetries# The comments must start with the symbol # .
# Empty lines are ignored.

#########################################################<SECTION-PARAMETERS>########################################################## Explanation of variables# ’model’ represents the type of the model that is# used for the calculations. It can be any name.
# ’lambda’ represents the value for the logarithmic# discretization parameter. Usually it is fixed between 2 and 3.
# ’max state no’ is the number of kept states# after an iteration.
# ’iteration no’ is the number of iterations performed# along the NRG run. Typically it is set to 40 - 70.
# ’symmetry no’ specifies how many symmetries# are used. The symmetries are characterized# in <SECTION-SYMMETRIES># ’coupling no’ specifies how many couplings# are present in the interaction Hamiltonian.
# ’local coupling no’ is the number of# Hamiltonians at the sites we add to the Wilson chain# ’spectral operator no’ is the number of spectral operators# acting on site 0, we keep track of.
# These spectral operators are characterized47

# in section <SECTION-SPECTRAL-OPERATORS># ’static operator no’ is the number of static operators# acting on site 0, we keep track of.
# These static operators are characterized# in section <SECTION-STATIC-OPERATORS># ’hopping operator no’ is the number of hopping operators.
# These hopping operators are characterized# in section <SECTION-HOPPING-OPERATORS># ’block state no’ specifies the number# of initial states/multiplets. These states are characterized# by the representation indices and are listed# in section <SECTION-BLOCK STATES>##################################################################################################################
model = kondo model # the name of the present model
lambda = 2.0 # discretization parametermax state no = 100 # maximum number of kept multiplets
iteration no = 20 # allowed number of iterations
symmetry no = 2 # the number of symmetries
coupling no = 2 # no of Hamiltonian terms, i.e.

J and B in this case

spectral operator no = 3 # no of operators for which correlation functions are computed
static operator no = 1 # no of operators for which the static average is computed.
hopping operator no = 2 # the number of hopping operators
block state no = 8 # the number of initial block states
local state no = 4 # the number of local states added in each iteration
local coupling no = 1 # the number of Hamiltonian terms on the added site
spectral function no = 3 # the number of spectral functions that need to be generated
interval no = 1000 # mash grid between the maximum and minimum energy at each iteration
degeneracy threshhold = 1e-6 # energy thresh-hold for discarding states
temperature = 0.0 # temperature used for the run########################################################</SECTION-PARAMETERS>

48

################################################################################################################<SECTION-FLAGS>######################################################### dmnrg flag controls whether the backward procedure# has to be done and the full set of eigen-states be# used for the calculation of the spectral function.
# ON - backwards iteration is done.
# OFF - only the nrg calculation is performed.
# text swap files flag = flag that fixes if the text mode is# used to save the files. The files will be saved in binary mode anyway.
# ON - the files are saved in the text mode and kept on the disk.
# OFF= the files are not saved in the text mode.
# binary swap files flag = flag that controls whether the# unnecessary binary files will be removed or not not after the calculation is done.
# ON - the files will be kept on the disk.
# OFF- the files will be removed from the disk.
# hoppings on site energies flag = flag that controls# the reading of the hoppings and the on-site energies.
# ON - the hoppings and the on site energies# are read from the files in the results/mapping/ folder# These files are generated by the he utility.
# OFF - the hoppings are computed on the fly# assuming a flat density of states (DOS = 0.5) on [-1, 1],# and the on-site energies are set to zero.
dmnrg flag = ON
text swap files flag = OFFbinary swap files flag = ON
hoppings on site energies flag = OFF########################################################</SECTION-FLAGS>#################################################################################################################<SECTION-SYMMETRIES>########################################################## As many symmetries must be defined below# as set by the variable ’symmetry no’# in the previous section.
# The possible symmetry types are the following:# U(1), SU(2), Z(2), charge SU(2).
## Next the limits for the representation index

49

# are set. All states with smaller/larger# representation indices will be cut off# For U(1) symmetries the representations are labeled# by integers (2 S z) and for SU(2) by non-negative# integers (2S).
# For the Z(2) symmetry the representations are labeled# by integers# 0 - even parity states# 1 - odd parity states########################################################## charge symmetry corresponding to Q
U(1) -10 10# spin symmetry corresponding to 2 S z
U(1) -20 20#########################################################</SECTION-SYMMETRIES>##################################################################################################################<SECTION-BLOCK STATES>########################################################## This section gives the representation indices of the# block states of the first iteration.
# The number of rows has to be equal# with the ’block state no’ in the first section,# and the number of columns to the ’symmetry no’.
# Each column contains indices corresponding to a symmetry.
#########################################################
−1 −1
−1 1
0
−2
0
0
0
0
0
2
1
−1
1
1
#########################################################</SECTION-BLOCK STATES>##################################################################################################################<SECTION-LOCAL STATES>########################################################## This section gives the representation indices# of the local states. The number of rows has to be equal# with the ’local state no’ in the first section50

# and the number of columns to the ’symmetry no’.
# Each column contains indices corresponding to a symmetry.
########################################################
−1 0
0
−1
0
1
1
0
########################################################</SECTION-LOCAL STATES>################################################################################################################<SECTION-LOCAL STATES SIGNS>######################################################### These are the signs of the added states.
# The number of signs has to be equal with# the number of ’local state no’.
# All the signs must be enumerated in one line.
# For even electron number, the sign is 1,# while for odd electron number the sign is -1.
#########################################################
1 -1 -1 1#########################################################</SECTION-LOCAL STATES SIGNS>##################################################################################################################<SECTION-BLOCK HAMILTONIAN># The total Hamiltonian will be written as a sum# of theP Hamiltonian terms.
H0 = δcoupling no Gδ Hδ
# There has to be as many <HAMILTONIAN TERM> items# as ’coupling no’ were set in the first section.
# The first tag is the name of the Hamiltonian term,# the second one gives the representation indices and# at last comes the matrix elements that will be# read into the sectors.

# first term in the Hamiltonian<BLOCK HAMILTONIAN TERM><BLOCK HAMILTONIAN NAME>
H Kondo</BLOCK HAMILTONIAN NAME><BLOCK HAMILTONIAN COUPLING># J - exchange couplingJ = 0.551

</BLOCK HAMILTONIAN COUPLING><BLOCK HAMILTONIAN REPRESENTATION INDEX>
0 0</BLOCK HAMILTONIAN REPRESENTATION INDEX><BLOCK HAMILTONIAN MATRIX>
0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0

0
0
0.25
0
0
0
0
0

0
0
0
−0.75
0
0
0
0

0
0
0
0
0.25
0
0
0

0
0
0
0
0
0.25
0
0

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0

</BLOCK HAMILTONIAN MATRIX></BLOCK HAMILTONIAN TERM># second term in the Hamiltonian<BLOCK HAMILTONIAN TERM><BLOCK HAMILTONIAN NAME>
H Zeeman</BLOCK HAMILTONIAN NAME><BLOCK HAMILTONIAN COUPLING># B - external magnetic field
B = 0.01</BLOCK HAMILTONIAN COUPLING><BLOCK HAMILTONIAN REPRESENTATION INDEX>
0 0</BLOCK HAMILTONIAN REPRESENTATION INDEX><BLOCK HAMILTONIAN MATRIX>
−0.5
0
0
0
0
0
0
0

0
0.5
0
0
0
0
0
0

0
0
−0.5
0
0
0
0
0

0
0
0
0
0.5
0
0
0

0
0
0
0.5
0
0
0
0

0
0
0
0
0
0.5
0
0

0
0
0
0
0
0
−0.5
0

0
0
0
0
0
0
0
0.5

</BLOCK HAMILTONIAN MATRIX>

52

</BLOCK HAMILTONIAN TERM>########################################################</SECTION-BLOCK HAMILTONIAN>################################################################################################################<SECTION-HOPPING OPERATORS># There must be as many <HOPPING OPERATOR> items# as ’hopping operator no’ were set in the first section# The matrices must contain the Reduced Matrix Elements.
######################################################### here is the first hopping operator<HOPPING OPERATOR><HOPPING OPERATOR NAME>f N dagger up</HOPPING OPERATOR NAME><HOPPING OPERATOR REPRESENTATION INDEX>
1 1</HOPPING OPERATOR REPRESENTATION INDEX><HOPPING OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</HOPPING OPERATOR SIGN><HOPPING OPERATOR MATRIX>
000−0.707106781186
0.707106781186
0
0
0

0
0
0
0
0
1
0
0

0
0
0
0
0
0
1
0

0
0
0
0
0
0
0
0.707106781186

0
0
0
0
0
0
0
0.707106781186</HOPPING OPERATOR MATRIX></HOPPING OPERATOR>

53

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0

# here is the second hopping operator<HOPPING OPERATOR><HOPPING OPERATOR NAME>f N dagger down</HOPPING OPERATOR NAME><HOPPING OPERATOR REPRESENTATION INDEX>
1 -1</HOPPING OPERATOR REPRESENTATION INDEX><HOPPING OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</HOPPING OPERATOR SIGN><HOPPING OPERATOR MATRIX>
0
0
1
0
0
0
0
0

0
0
0
0.707106781186
0.707106781186
0
0
0

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0.707106781186
0

0
0
0
0
0
0
−0.707106781186
0

0
0
0
0
0
0
0
−1

</HOPPING OPERATOR MATRIX></HOPPING OPERATOR>########################################################</SECTION-HOPPING OPERATORS>################################################################################################################<SECTION-SPECTRAL OPERATORS>######################################################### In this section we set up the matrices for the# operators for which the spectral function is computed# Matrices must contain the reduced matrix elements.
######################################################## here comes the first spectral operator.

54

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0

<SPECTRAL OPERATOR><SPECTRAL OPERATOR NAME>
Sz</SPECTRAL OPERATOR NAME><SPECTRAL OPERATOR REPRESENTATION INDEX>
0 0</SPECTRAL OPERATOR REPRESENTATION INDEX><SPECTRAL OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
1</SPECTRAL OPERATOR SIGN><SPECTRAL OPERATOR MATRIX>
−0.5
0
0
0
0
0
0
0

0
0.5
0
0
0
0
0
0

0
0
−0.5
0
0
0
0
0

0
0
0
0
0.5
0
0
0

0
0
0
0.5
0
0
0
0

0
0
0
0
0
0.5
0
0

0
0
0
0
0
0
−0.5
0

0
0
0
0
0
0
0
0.5

</SPECTRAL OPERATOR MATRIX></SPECTRAL OPERATOR># here comes the second spectral operator.
<SPECTRAL OPERATOR><SPECTRAL OPERATOR NAME>
f 0 dagger up</SPECTRAL OPERATOR NAME><SPECTRAL OPERATOR REPRESENTATION INDEX>
1 1</SPECTRAL OPERATOR REPRESENTATION INDEX><SPECTRAL OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</SPECTRAL OPERATOR SIGN><SPECTRAL OPERATOR MATRIX>

55

0
0
0
−0.707106781186
0.707106781186
0
0
0

0
0
0
0
0
1
0
0

0
0
0
0
0
0
1
0

0
0
0
0
0
0
0
0.707106781186

0
0
0
0
0
0
0
0.707106781186

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0

</SPECTRAL OPERATOR MATRIX></SPECTRAL OPERATOR># here comes the third spectral operator.
<SPECTRAL OPERATOR><SPECTRAL OPERATOR NAME>
f 0 dagger down</SPECTRAL OPERATOR NAME><SPECTRAL OPERATOR REPRESENTATION INDEX>
1 -1</SPECTRAL OPERATOR REPRESENTATION INDEX><SPECTRAL OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</SPECTRAL OPERATOR SIGN><SPECTRAL OPERATOR MATRIX>
0
0
1
0
0
0
0
0

0
0
0
0.707106781186
0.707106781186
0
0
0

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0.707106781186
0

0
0
0
0
0
0
−0.707106781186
0

0
0
0
0
0
0
0
−1

</SPECTRAL OPERATOR MATRIX></SPECTRAL OPERATOR>########################################################</SECTION-SPECTRAL OPERATORS>################################################################################################################<SECTION-STATIC OPERATORS>

56

0
0
0
0
0
0
0
0

0
0
0
0
0
0
0
0

######################################################### In this section we set up the matrices for the# operators for which the static values are computed# These operators need to have representation# indices always zero!
#######################################################<STATIC OPERATOR><STATIC OPERATOR NAME>Sz</STATIC OPERATOR NAME><STATIC OPERATOR REPRESENTATION INDEX>
0 0</STATIC OPERATOR REPRESENTATION INDEX><STATIC OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
1</STATIC OPERATOR SIGN><STATIC OPERATOR MATRIX>
−0.5
0
0
0
0
0
0
0

0
0.5
0
0
0
0
0
0

0
0
−0.5
0
0
0
0
0

0
0
0
0
0.5
0
0
0

0
0
0
0.5
0
0
0
0

0
0
0
0
0
0.5
0
0

0
0
0
0
0
0
−0.5
0

0
0
0
0
0
0
0
0.5

</STATIC OPERATOR MATRIX></STATIC OPERATOR>########################################################</SECTION-STATIC OPERATORS>################################################################################################################<SECTION-LOCAL HOPPING OPERATORS>########################################################<LOCAL HOPPING OPERATOR># here comes the first local hopping operator<LOCAL HOPPING OPERATOR NAME>
f N+1 dagger up

57

</LOCAL HOPPING OPERATOR NAME><LOCAL HOPPING OPERATOR REPRESENTATION INDEX>
1 1</LOCAL HOPPING OPERATOR REPRESENTATION INDEX><LOCAL HOPPING OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</LOCAL HOPPING OPERATOR SIGN><LOCAL HOPPING OPERATOR MATRIX>
0
0
1
0

0
0
0
1

0
0
0
0

0
0
0
0

</LOCAL HOPPING OPERATOR MATRIX><LOCAL ON SITE ENERGY MATRIX>
0
0
0
0

0
0
0
0

0
0
1
0

0
0
0
1

</LOCAL ON SITE ENERGY MATRIX></LOCAL HOPPING OPERATOR># here comes the second local hopping operator<LOCAL HOPPING OPERATOR><LOCAL HOPPING OPERATOR NAME>
f N+1 dagger down</LOCAL HOPPING OPERATOR NAME><LOCAL HOPPING OPERATOR REPRESENTATION INDEX>
1 -1</LOCAL HOPPING OPERATOR REPRESENTATION INDEX><LOCAL HOPPING OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</LOCAL HOPPING OPERATOR SIGN>

58

<LOCAL HOPPING OPERATOR MATRIX>
0
1
0
0

0
0
0
0

0
0
0
−1

0
0
0
0

</LOCAL HOPPING OPERATOR MATRIX><LOCAL ON SITE ENERGY MATRIX>
0
0
0
0

0
1
0
0

0
0
0
0

0
0
0
1

</LOCAL ON SITE ENERGY MATRIX></LOCAL HOPPING OPERATOR>########################################################</SECTION-LOCAL HOPPING OPERATORS>################################################################################################################<SECTION-LOCAL HAMILTONIAN># The total Local Hamiltonian of the added site
# will be written as a sum of Local Hamiltonian terms# acting
P on site n:
Hn = i gi,n Oi,n# The number of <LOCAL HAMILTONIAN TERM> items# must be equal with the variable ’local coupling no’# set in the first section# The first tag is the name of the local Hamiltonian term
# the second one gives the representation indices# last come the matrix elements.
######################################################### first term in the Local Hamiltonian#######################################################<LOCAL HAMILTONIAN TERM><LOCAL HAMILTONIAN NAME>
H1</LOCAL HAMILTONIAN NAME><LOCAL HAMILTONIAN COUPLING>#irrelevant
0.0</LOCAL HAMILTONIAN COUPLING><LOCAL HAMILTONIAN REPRESENTATION INDEX>

59

0 0
</LOCAL HAMILTONIAN REPRESENTATION INDEX><LOCAL HAMILTONIAN MATRIX>
1
0
0
0

0
0
0
0

0
0
1
0

0
0
0
0

</LOCAL HAMILTONIAN MATRIX></LOCAL HAMILTONIAN TERM></SECTION-LOCAL HAMILTONIAN>########################################################<SECTION-SPECTRAL FUNCTION>######################################################### In this section we set up the spectral function# calculation for different types of operators.
# All the time we provide in the input file the# matrices for the A† and B † operators# but the spectral function is computed for the# combination <A, B † ># spectral function for <f 0 up, f 0 dagger up>{f 0 dagger up; f 0 dagger up}# spectral function for <f 0 down, f 0 dagger down>{f 0 dagger down; f 0 dagger down}#spectral function for <S z, S z>
{S z; S z}########################################################</SECTION-SPECTRAL FUNCTION>######################################################### The section below is not needed by the fnrg binary.
# It is needed when doing the broadening of the spectral function# with the utility sfb.

########################################################<SECTION-SPECTRAL FUNCTION BROADENING>######################################################### In this section we include broadening parameters that# are used for computing the broaden spectral function# The methods that can be used are either LOG GAUSS or# INTERPOLATIVE LOG GAUSS.

60

# These procedures are described in the manual.
# The quantum temperature is relevant only
# for the INTERPOLATIVE LOG GAUSS method,# otherwise only the ’broadening parameter’ parameter# setting the width of the log-normal# distribution function is relevant.
# The quantum temperature must be always lower than T# for finite temperature calculations.
# A good choice is to set the quatum temperature in the range# 0.5T - T.
# The three parameters spectral function *# specify which spectral functions are being calculated.
# The accepted values are YES or NO only.
# The static average dmnrg flags specifies whether the expectation values# of the static operators are computed.
#The accepted values are YES or NO only.
# The green function dmnrg specifies whether the real part of the# Green’s function is computed from the spectral function data.
broadening method = INTERPOLATIVE LOG GAUSS # method used either LOG GAUSS or INTERPOLATIVE LOG GAUSS
broadening parameter = 0.70 # parameter used for the log gaussian distribution
broadening energy minim = 1e-6 # minimum energy for which the broadening is performed
broadening energy maxim = 4.0 # maximum energy for which the broadening is performed
broadening grid mesh = 100 # grid mesh on a log scale between the minimum and maximum energies
quantum temperature = 1e-7 # used for the INTERPOLATIVE LOG GAUSS method only.
spectral function nrg even = YES # spectral function for the even iterations is computed
spectral function nrg odd = YES # spectral function for the odd iterations part is computed
spectral function dmnrg = YES # spectral function for the dmnrg calculation is computed.
static average dmnrg = YES # average functions similar with the case of spectral function.
green function nrg even = YES # calculation of the real/imaginary part of the Green’s function
from the even part of the nrg spectral function
green function nrg odd = YES # calculation of the real/imaginary part of the Green’s function
from the odd part of the nrg spectral function
green function dmnrg = NO # calculation of the real/imaginary part of the Green’s function
from the spectral function
green function grid mesh = 400 # points for the mesh when doing Hilbert transform########################################################</SECTION-SPECTRAL FUNCTION BROADENING>########################################################

61

62

Appendix B

Input file for the Kondo modelwith Ucharge(1) × SUspin(2)symmetries# The comments must start with the symbol # .
# Empty lines are ignored.

#########################################################<SECTION-PARAMETERS>########################################################## Explanation of variables# ’model’ represents the type of the model that is# used for the calculations. It can be any name.
# ’lambda’ represents the value for the logarithmic# discretization parameter. Usually it is fixed between 2 and 3.
# ’max state no’ is the number of kept states# after an iteration.
# ’iteration no’ is the number of iterations performed# along the NRG run. Typically it is set to 40 - 70.
# ’symmetry no’ specifies how many symmetries# are used. The symmetries are characterized# in <SECTION-SYMMETRIES># ’coupling no’ specifies how many couplings# are present in the interaction Hamiltonian.
# ’local coupling no’ is the number of# Hamiltonians at the sites we add to the Wilson chain# ’spectral operator no’ is the number of spectral operators# acting on site 0, we keep track of.
# These spectral operators are characterized63

# in section <SECTION-SPECTRAL-OPERATORS># ’static operator no’ is the number of static operators# acting on site 0, we keep track of.
# These static operators are characterized# in section <SECTION-STATIC-OPERATORS># ’hopping operator no’ is the number of hopping operators.
# These hopping operators are characterized# in section <SECTION-HOPPING-OPERATORS># ’block state no’ specifies the number# of initial states/multiplets. These states are characterized# by the representation indices and are listed# in section <SECTION-BLOCK STATES>##################################################################################################################
model = kondo model # the name of the present model
lambda = 2.0 # discretization parametermax state no = 100 # maximum number of kept multiplets
iteration no = 20 # allowed number of iterations
symmetry no = 2 # the number of symmetries
coupling no = 1 # no of Hamiltonian terms, i.e.

J in this case

spectral operator no = 2 # no of operators for which correlation functions are computed
static operator no = 0 # no of operators for which the static average is computed.
hopping operator no = 1 # the number of hopping operators
block state no = 4 # the number of initial block states
local state no = 3 # the number of local states added in each iteration
local coupling no = 1 # the number of Hamiltonian terms on the added site
spectral function no = 2 # the number of spectral functions that need to be generated
interval no = 1000 # mash grid between the maximum and minimum energy at each iteration

degeneracy threshhold = 1e-6 # energy thresh-hold for discarding states
temperature = 0.0 # temperature used for the run########################################################

64

</SECTION-PARAMETERS>################################################################################################################<SECTION-FLAGS>######################################################### dmnrg flag controls whether the backward procedure# has to be done and the full set of eigen-states be# used for the calculation of the spectral function.
# ON - backwards iteration is done.
# OFF - only the nrg calculation is performed.
# text swap files flag = flag that fixes if the text mode is# used to save the files. The files will be saved in binary mode anyway.
# ON - the files are saved in the text mode and kept on the disk.
# OFF= the files are not saved in the text mode.
# binary swap files flag = flag that controls whether the# unnecessary binary files will be removed or not not after the calculation is done.
# ON - the files will be kept on the disk.
# OFF- the files will be removed from the disk.
# hoppings on site energies flag = flag that controls# the reading of the hoppings and the on-site energies.
# ON - the hoppings and the on site energies# are read from the files in the results/mapping/ folder# These files are generated by the he utility.
# OFF - the hoppings are computed on the fly# assuming a flat density of states (DOS = 0.5) on [-1, 1],# and the on-site energies are set to zero.
dmnrg flag = ON
text swap files flag = OFFbinary swap files flag = ON
hoppings on site energies flag =OFF########################################################</SECTION-FLAGS>#################################################################################################################<SECTION-SYMMETRIES>########################################################## As many symmetries must be defined below# as set by the variable ’symmetry no’# in the previous section.
# The possible symmetry types are the following:# U(1), SU(2), Z(2), charge SU(2).
#65# Next the limits for the representation index# are set. All states with smaller/larger# representation indices will be cut off# For U(1) symmetries the representations are labeled# by integers (2 S z) and for SU(2) by non-negative# integers (2S).
# For the Z(2) symmetry the representations are labeled# by integers# 0 - even parity states# 1 - odd parity states########################################################## charge symmetry corresponding to Q
U(1) -10 10# spin symmetry corresponding to S
SU(2) 0 20#########################################################</SECTION-SYMMETRIES>##################################################################################################################<SECTION-BLOCK STATES>########################################################## This section gives the representation indices of the# block states of the first iteration.
# The number of rows has to be equal# with the ’block state no’ in the first section,# and the number of columns to the ’symmetry no’.
# Each column contains indices corresponding to a symmetry.
#########################################################
−1 1
0
0
0
2
1
1
#########################################################</SECTION-BLOCK STATES>##################################################################################################################<SECTION-LOCAL STATES>########################################################## This section gives the representation indices# of the local states. The number of rows has to be equal# with the ’local state no’ in the first section# and the number of columns to the ’symmetry no’.
# Each column contains indices corresponding to a symmetry.
########################################################

66

−1 0
0
1
1
0
########################################################</SECTION-LOCAL STATES>################################################################################################################<SECTION-LOCAL STATES SIGNS>######################################################### These are the signs of the added states.
# The number of signs has to be equal with# the number of ’local state no’.
# All the signs must be enumerated in one line.
# For even electron number, the sign is 1,# while for odd electron number the sign is -1.
#########################################################
1 -1 1#########################################################</SECTION-LOCAL STATES SIGNS>##################################################################################################################<SECTION-BLOCK HAMILTONIAN># The total Hamiltonian will be written as a sum# of theP Hamiltonian terms.
H0 = δcoupling no Gδ Hδ
# There has to be as many <HAMILTONIAN TERM> items# as ’coupling no’ were set in the first section.
# The first tag is the name of the Hamiltonian term,# the second one gives the representation indices and# at last comes the matrix elements that will be# read into the sectors.

# first term in the Hamiltonian<BLOCK HAMILTONIAN TERM><BLOCK HAMILTONIAN NAME>
H Kondo</BLOCK HAMILTONIAN NAME><BLOCK HAMILTONIAN COUPLING># J - exchange couplingJ = 0.5</BLOCK HAMILTONIAN COUPLING>

67

<BLOCK HAMILTONIAN REPRESENTATION INDEX>
0 0</BLOCK HAMILTONIAN REPRESENTATION INDEX><BLOCK HAMILTONIAN MATRIX>
0
0
0
0

0
−0.75
0
0

0
0
0.25
0

0
0
0
0

</BLOCK HAMILTONIAN MATRIX></BLOCK HAMILTONIAN TERM>########################################################</SECTION-BLOCK HAMILTONIAN>################################################################################################################<SECTION-HOPPING OPERATORS># There must be as many <HOPPING OPERATOR> items# as ’hopping operator no’ were set in the first section# The matrices must contain the Reduced Matrix Elements.
######################################################### here is the first hopping operator<HOPPING OPERATOR><HOPPING OPERATOR NAME>f N dagger</HOPPING OPERATOR NAME><HOPPING OPERATOR REPRESENTATION INDEX>
1 1</HOPPING OPERATOR REPRESENTATION INDEX><HOPPING OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</HOPPING OPERATOR SIGN><HOPPING OPERATOR MATRIX>

68

0
1
1
0

0
0
0
0.707106781186

000−1.224744871392

0
0
0
0

</HOPPING OPERATOR MATRIX></HOPPING OPERATOR>########################################################</SECTION-HOPPING OPERATORS>################################################################################################################<SECTION-SPECTRAL OPERATORS>######################################################### In this section we set up the matrices for the# operators for which the spectral function is computed# Matrices must contain the reduced matrix elements.
######################################################## here comes the first spectral operator.
<SPECTRAL OPERATOR><SPECTRAL OPERATOR NAME>
S</SPECTRAL OPERATOR NAME><SPECTRAL OPERATOR REPRESENTATION INDEX>
0 2</SPECTRAL OPERATOR REPRESENTATION INDEX><SPECTRAL OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
1</SPECTRAL OPERATOR SIGN><SPECTRAL OPERATOR MATRIX>
0.866025403784
0
0
0

0
0
0.5
0

0
−0.866025403784
0.707106781186
00000.866025403784</SPECTRAL OPERATOR MATRIX></SPECTRAL OPERATOR># here comes the second spectral operator.
<SPECTRAL OPERATOR>

69

<SPECTRAL OPERATOR NAME>
f 0 dagger</SPECTRAL OPERATOR NAME><SPECTRAL OPERATOR REPRESENTATION INDEX>
1 1</SPECTRAL OPERATOR REPRESENTATION INDEX><SPECTRAL OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</SPECTRAL OPERATOR SIGN><SPECTRAL OPERATOR MATRIX>
0
1
1
0

0
0
0
0.707106781186

000−1.224744871392

0
0
0
0

</SPECTRAL OPERATOR MATRIX></SPECTRAL OPERATOR>########################################################</SECTION-SPECTRAL OPERATORS>################################################################################################################<SECTION-LOCAL HOPPING OPERATORS>########################################################<LOCAL HOPPING OPERATOR># here comes the first local hopping operator<LOCAL HOPPING OPERATOR NAME>
f N+1 dagger</LOCAL HOPPING OPERATOR NAME><LOCAL HOPPING OPERATOR REPRESENTATION INDEX>
1 1</LOCAL HOPPING OPERATOR REPRESENTATION INDEX><LOCAL HOPPING OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.

70

-1</LOCAL HOPPING OPERATOR SIGN><LOCAL HOPPING OPERATOR MATRIX>
0
1
0

0
0
−1.4142135623731

000</LOCAL HOPPING OPERATOR MATRIX><LOCAL ON SITE ENERGY MATRIX>
0
0
0

0
0
0

0
0
0

</LOCAL ON SITE ENERGY MATRIX></LOCAL HOPPING OPERATOR>########################################################</SECTION-LOCAL HOPPING OPERATORS>################################################################################################################<SECTION-LOCAL HAMILTONIAN># The total Local Hamiltonian of the added site
# will be written as a sum of Local Hamiltonian terms# acting
P on site n:
Hn = i gi,n Oi,n# The number of <LOCAL HAMILTONIAN TERM> items# must be equal with the variable ’local coupling no’# set in the first section# The first tag is the name of the local Hamiltonian term
# the second one gives the representation indices# last come the matrix elements.
######################################################### first term in the Local Hamiltonian#######################################################<LOCAL HAMILTONIAN TERM><LOCAL HAMILTONIAN NAME>
H1</LOCAL HAMILTONIAN NAME><LOCAL HAMILTONIAN COUPLING>#irrelevant
0.071</LOCAL HAMILTONIAN COUPLING><LOCAL HAMILTONIAN REPRESENTATION INDEX>
0 0</LOCAL HAMILTONIAN REPRESENTATION INDEX><LOCAL HAMILTONIAN MATRIX>
1
0
0

0
0
0

0
0
1

</LOCAL HAMILTONIAN MATRIX></LOCAL HAMILTONIAN TERM></SECTION-LOCAL HAMILTONIAN>########################################################<SECTION-SPECTRAL FUNCTION>######################################################### In this section we set up the spectral function# calculation for different types of operators.
# All the time we provide in the input file the# matrices for the A† and B † operators# but the spectral function is computed for the# combination <A, B † ># spectral function for <f 0, f 0 dagger>{f 0 dagger; f 0 dagger }#spectral function for spin <S, S>{S; S}########################################################</SECTION-SPECTRAL FUNCTION>######################################################### The section below is not needed by the fnrg binary.
# It is needed when doing the broadening of the spectral function# with the utility sfb.

########################################################<SECTION-SPECTRAL FUNCTION BROADENING>######################################################### In this section we include broadening parameters that# are used for computing the broaden spectral function# The methods that can be used are either LOG GAUSS or# INTERPOLATIVE LOG GAUSS.
# These procedures are described in the manual.

72

# The quantum temperature is relevant only
# for the INTERPOLATIVE LOG GAUSS method,# otherwise only the ’broadening parameter’ parameter# setting the width of the log-normal# distribution function is relevant.
# The quantum temperature must be always lower than T# for finite temperature calculations.
# A good choice is to set the quatum temperature in the range# 0.5T - T.
# The three parameters spectral function *# specify which spectral functions are being calculated.
# The accepted values are YES or NO only.
# The static average dmnrg flags specifies whether the expectation values# of the static operators are computed.
# The accepted values are YES or NO only.
# The green function dmnrg specifies whether the real part of the# Green’s function is computed from the spectral function data.
broadening method = INTERPOLATIVE LOG GAUSS # method used either LOG GAUSS or INTERPOLATIVE LOG GAUSS
broadening parameter = 0.70 # parameter used for the log gaussian distribution
broadening energy minim = 1e-6 # minimum energy for which the broadening is performed
broadening energy maxim = 4.0 # maximum energy for which the broadening is performed
broadening grid mesh = 100 # grid mesh on a log scale between the minimum and maximum energies
quantum temperature = 1e-7 # used for the INTERPOLATIVE LOG GAUSS method only.
spectral function nrg even = YES # spectral function for the even iterations is computed
spectral function nrg odd = YES # spectral function for the odd iterations part is computed
spectral function dmnrg = YES # spectral function for the dmnrg calculation is computed.
static average dmnrg = YES # average functions similar with the case of spectral function.
green function nrg even = YES # calculation of the real/imaginary part of the Green’s function
from the even part of the nrg spectral function
green function nrg odd = YES # calculation of the real/imaginary part of the Green’s function
from the odd part of the nrg spectral function
green function dmnrg = NO # calculation of the real/imaginary part of the Green’s function
from the spectral function
green function grid mesh = 400 # points for the mesh when doing Hilbert transform########################################################</SECTION-SPECTRAL FUNCTION BROADENING>########################################################

73

74

Appendix C

Input file for the Kondo modelwith SUcharge(2) × SUspin(2)symmetries# The comments must start with the symbol # .
# Empty lines are ignored.

#########################################################<SECTION-PARAMETERS>########################################################## Explanation of variables# ’model’ represents the type of the model that is# used for the calculations. It can be any name.
# ’lambda’ represents the value for the logarithmic# discretization parameter. Usually it is fixed between 2 and 3.
# ’max state no’ is the number of kept states# after an iteration.
# ’iteration no’ is the number of iterations performed# along the NRG run. Typically it is set to 40 - 70.
# ’symmetry no’ specifies how many symmetries# are used. The symmetries are characterized# in <SECTION-SYMMETRIES># ’coupling no’ specifies how many couplings# are present in the interaction Hamiltonian.
# ’local coupling no’ is the number of# Hamiltonians at the sites we add to the Wilson chain# ’spectral operator no’ is the number of spectral operators# acting on site 0, we keep track of.
# These spectral operators are characterized75

# in section <SECTION-SPECTRAL-OPERATORS># ’static operator no’ is the number of static operators# acting on site 0, we keep track of.
# These static operators are characterized# in section <SECTION-STATIC-OPERATORS># ’hopping operator no’ is the number of hopping operators.
# These hopping operators are characterized# in section <SECTION-HOPPING-OPERATORS># ’block state no’ specifies the number# of initial states/multiplets. These states are characterized# by the representation indices and are listed# in section <SECTION-BLOCK STATES>##################################################################################################################
model = kondo model # the name of the present model
lambda = 2.0 # discretization parametermax state no = 100 # maximum number of kept multiplets
iteration no = 20 # allowed number of iterations
symmetry no = 2 # the number of symmetries
coupling no = 1 # no of Hamiltonian terms, i.e.

J in this case

spectral operator no = 2 # no of operators for which correlation functions are computed
static operator no = 0 # no of operators for which the static average is computed.
hopping operator no = 1 # the number of hopping operators
block state no = 3 # the number of initial block states
local state no = 2 # the number of local states added in each iteration
local coupling no = 1 # the number of Hamiltonian terms on the added site
spectral function no = 2 # the number of spectral functions that need to be generated
interval no = 1000 # mash grid between the maximum and minimum energy at each iteration
degeneracy threshhold = 1e-6 # energy thresh-hold for discarding states
temperature = 0.0 # temperature used for the run########################################################</SECTION-PARAMETERS>

76

################################################################################################################<SECTION-FLAGS>######################################################### dmnrg flag controls whether the backward procedure# has to be done and the full set of eigen-states be# used for the calculation of the spectral function.
# ON - backwards iteration is done.
# OFF - only the nrg calculation is performed.
# text swap files flag = flag that fixes if the text mode is# used to save the files. The files will be saved in binary mode anyway.
# ON - the files are saved in the text mode and kept on the disk.
# OFF= the files are not saved in the text mode.
# binary swap files flag = flag that controls whether the# unnecessary binary files will be removed or not not after the calculation is done.
# ON - the files will be kept on the disk.
# OFF- the files will be removed from the disk.
# hoppings on site energies flag = flag that controls# the reading of the hoppings and the on-site energies.
# ON - the hoppings and the on site energies# are read from the files in the results/mapping/ folder# These files are generated by the he utility.
# OFF - the hoppings are computed on the fly# assuming a flat density of states (DOS = 0.5) on [-1, 1],# and the on-site energies are set to zero.
dmnrg flag = ON
text swap files flag = OFFbinary swap files flag = ON
hoppings on site energies flag = OFF########################################################</SECTION-FLAGS>#################################################################################################################<SECTION-SYMMETRIES>########################################################## As many symmetries must be defined below# as set by the variable ’symmetry no’# in the previous section.
# The possible symmetry types are the following:# U(1), SU(2), Z(2), charge SU(2).
## Next the limits for the representation index

77

# are set. All states with smaller/larger# representation indices will be cut off# For U(1) symmetries the representations are labeled# by integers (2 S z) and for SU(2) by non-negative# integers (2S).

# For the Z(2) symmetry the representations are labeled# by integers# 0 - even parity states# 1 - odd parity states########################################################## charge symmetry corresponding to Q
charge SU(2) 0 20
# spin symmetry corresponding to S
SU(2) 0 20#########################################################</SECTION-SYMMETRIES>##################################################################################################################<SECTION-BLOCK STATES>########################################################## This section gives the representation indices of the# block states of the first iteration.
# The number of rows has to be equal# with the ’block state no’ in the first section,# and the number of columns to the ’symmetry no’.
# Each column contains indices corresponding to a symmetry.
#########################################################
0 0
0 2
1 1
#########################################################</SECTION-BLOCK STATES>##################################################################################################################<SECTION-LOCAL STATES>########################################################## This section gives the representation indices# of the local states. The number of rows has to be equal# with the ’local state no’ in the first section# and the number of columns to the ’symmetry no’.
# Each column contains indices corresponding to a symmetry.
########################################################
0 1
1 0

78

########################################################</SECTION-LOCAL STATES>################################################################################################################<SECTION-LOCAL STATES SIGNS>######################################################### These are the signs of the added states.
# The number of signs has to be equal with# the number of ’local state no’.
# All the signs must be enumerated in one line.
# For even electron number, the sign is 1,# while for odd electron number the sign is -1.
#########################################################-1 1#########################################################</SECTION-LOCAL STATES SIGNS>##################################################################################################################<SECTION-BLOCK HAMILTONIAN># The total Hamiltonian will be written as a sum# of theP Hamiltonian terms.
H0 = δcoupling no Gδ Hδ
# There has to be as many <HAMILTONIAN TERM> items# as ’coupling no’ were set in the first section.
# The first tag is the name of the Hamiltonian term,# the second one gives the representation indices and# at last comes the matrix elements that will be# read into the sectors.

# first term in the Hamiltonian<BLOCK HAMILTONIAN TERM><BLOCK HAMILTONIAN NAME>
H Kondo</BLOCK HAMILTONIAN NAME><BLOCK HAMILTONIAN COUPLING># J - exchange couplingJ = 0.5</BLOCK HAMILTONIAN COUPLING><BLOCK HAMILTONIAN REPRESENTATION INDEX>

79

0 0
</BLOCK HAMILTONIAN REPRESENTATION INDEX><BLOCK HAMILTONIAN MATRIX>
−0.75
0
0

0
0.25
0

0
0
0

</BLOCK HAMILTONIAN MATRIX></BLOCK HAMILTONIAN TERM>########################################################</SECTION-BLOCK HAMILTONIAN>################################################################################################################<SECTION-HOPPING OPERATORS># There must be as many <HOPPING OPERATOR> items# as ’hopping operator no’ were set in the first section# The matrices must contain the Reduced Matrix Elements.
######################################################### here is the first hopping operator<HOPPING OPERATOR><HOPPING OPERATOR NAME>f N dagger</HOPPING OPERATOR NAME><HOPPING OPERATOR REPRESENTATION INDEX>
1 1</HOPPING OPERATOR REPRESENTATION INDEX><HOPPING OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</HOPPING OPERATOR SIGN><HOPPING OPERATOR MATRIX>
000.707106781186

00−1.224744871391−1.4142135623731−1.4142135623731
080

</HOPPING OPERATOR MATRIX></HOPPING OPERATOR>########################################################</SECTION-HOPPING OPERATORS>################################################################################################################<SECTION-SPECTRAL OPERATORS>######################################################### In this section we set up the matrices for the# operators for which the spectral function is computed# Matrices must contain the reduced matrix elements.
######################################################## here comes the first spectral operator.
<SPECTRAL OPERATOR># here comes the spectral operator.
<SPECTRAL OPERATOR><SPECTRAL OPERATOR NAME>
f 0 dagger</SPECTRAL OPERATOR NAME><SPECTRAL OPERATOR REPRESENTATION INDEX>
1 1</SPECTRAL OPERATOR REPRESENTATION INDEX><SPECTRAL OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</SPECTRAL OPERATOR SIGN><SPECTRAL OPERATOR MATRIX>
00.707106781186

00−1.224744871391−1.4142135623731−1.4142135623731
0</SPECTRAL OPERATOR MATRIX></SPECTRAL OPERATOR># here comes the second spectral operator.
<SPECTRAL OPERATOR>

81

<SPECTRAL OPERATOR><SPECTRAL OPERATOR NAME>
S</SPECTRAL OPERATOR NAME><SPECTRAL OPERATOR REPRESENTATION INDEX>
0 2</SPECTRAL OPERATOR REPRESENTATION INDEX><SPECTRAL OPERATOR SIGN># The sign is 1 for bosonic operators and -1 for fermionic operators.
1</SPECTRAL OPERATOR SIGN><SPECTRAL OPERATOR MATRIX>
00.5
0−0.866025403784
0.707106781186
0000.866025403784</SPECTRAL OPERATOR MATRIX></SPECTRAL OPERATOR>########################################################</SECTION-SPECTRAL OPERATORS>################################################################################################################<SECTION-LOCAL HOPPING OPERATORS>########################################################<LOCAL HOPPING OPERATOR># here comes the first local hopping operator<LOCAL HOPPING OPERATOR NAME>
f N+1 dagger</LOCAL HOPPING OPERATOR NAME><LOCAL HOPPING OPERATOR REPRESENTATION INDEX>
1 1</LOCAL HOPPING OPERATOR REPRESENTATION INDEX><LOCAL HOPPING OPERATOR SIGN>

82

# The sign is 1 for bosonic operators and -1 for fermionic operators.
-1</LOCAL HOPPING OPERATOR SIGN><LOCAL HOPPING OPERATOR MATRIX>
0−1.4142135623731

1.4142135623731
0</LOCAL HOPPING OPERATOR MATRIX><LOCAL ON SITE ENERGY MATRIX>
0
0

0
0

</LOCAL ON SITE ENERGY MATRIX></LOCAL HOPPING OPERATOR>########################################################</SECTION-LOCAL HOPPING OPERATORS>################################################################################################################<SECTION-LOCAL HAMILTONIAN># The total Local Hamiltonian of the added site
# will be written as a sum of Local Hamiltonian terms# acting
P on site n:
Hn = i gi,n Oi,n# The number of <LOCAL HAMILTONIAN TERM> items# must be equal with the variable ’local coupling no’# set in the first section# The first tag is the name of the local Hamiltonian term
# the second one gives the representation indices# last come the matrix elements.
######################################################### first term in the Local Hamiltonian#######################################################<LOCAL HAMILTONIAN TERM><LOCAL HAMILTONIAN NAME>
H1</LOCAL HAMILTONIAN NAME><LOCAL HAMILTONIAN COUPLING>#irrelevant
0.0</LOCAL HAMILTONIAN COUPLING>

83

<LOCAL HAMILTONIAN REPRESENTATION INDEX>
0 0</LOCAL HAMILTONIAN REPRESENTATION INDEX><LOCAL HAMILTONIAN MATRIX>
1
0

0
0

</LOCAL HAMILTONIAN MATRIX></LOCAL HAMILTONIAN TERM></SECTION-LOCAL HAMILTONIAN>########################################################<SECTION-SPECTRAL FUNCTION>######################################################### In this section we set up the spectral function# calculation for different types of operators.
# All the time we provide in the input file the# matrices for the A† and B † operators# but the spectral function is computed for the# combination <A, B † ># spectral function for <f 0, f 0 dagger>{f 0 dagger; f 0 dagger }# spectral function for spin <S, S>{S; S }########################################################</SECTION-SPECTRAL FUNCTION>######################################################### The section below is not needed by the fnrg binary.
# It is needed when doing the broadening of the spectral function# with the utility sfb.

########################################################<SECTION-SPECTRAL FUNCTION BROADENING>######################################################### In this section we include broadening parameters that# are used for computing the broaden spectral function# The methods that can be used are either LOG GAUSS or# INTERPOLATIVE LOG GAUSS.
# These procedures are described in the manual.
# The quantum temperature is relevant only
# for the INTERPOLATIVE LOG GAUSS method,84

# otherwise only the ’broadening parameter’ parameter# setting the width of the log-normal# distribution function is relevant.
# The quantum temperature must be always lower than T# for finite temperature calculations.
# A good choice is to set the quatum temperature in the range# 0.5T - T.
# The three parameters spectral function *# specify which spectral functions are being calculated.
# The accepted values are YES or NO only.
# The static average dmnrg flags specifies whether the expectation values# of the static operators are computed.
# The accepted values are YES or NO only.
# The green function dmnrg specifies whether the real part of the# Green’s function is computed from the spectral function data.
broadening method = INTERPOLATIVE LOG GAUSS # method used either LOG GAUSS or INTERPOLATIVE LOG GAUSS
broadening parameter = 0.70 # parameter used for the log gaussian distribution
broadening energy minim = 1e-6 # minimum energy for which the broadening is performed
broadening energy maxim = 4.0 # maximum energy for which the broadening is performed
broadening grid mesh = 100 # grid mesh on a log scale between the minimum and maximum energies
quantum temperature = 1e-7 # used for the INTERPOLATIVE LOG GAUSS method only.
spectral function nrg even = YES # spectral function for the even iterations is computed
spectral function nrg odd = YES # spectral function for the odd iterations part is computed
spectral function dmnrg = YES # spectral function for the dmnrg calculation is computed.
static average dmnrg = YES # average functions similar with the case of spectral function.
green function nrg even = YES # calculation of the real/imaginary part of the Green’s function
from the even part of the nrg spectral function
green function nrg odd = YES # calculation of the real/imaginary part of the Green’s function
from the odd part of the nrg spectral function
green function dmnrg = NO # calculation of the real/imaginary part of the Green’s function
from the spectral function
green function grid mesh = 400 # points for the mesh when doing Hilbert transform########################################################</SECTION-SPECTRAL FUNCTION BROADENING>########################################################

85

86

Appendix D

License agreements
GNU LESSER GENERAL PUBLIC LICENSE Version 3, 29 June 2007
Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/> Everyone is permitted
to copy and distribute verbatim copies of this license document, but changing it is not allowed.
This version of the GNU Lesser General Public License incorporates the terms and conditionsof version 3 of the GNU General Public License, supplemented by the additional permissions
listed below.
0. Additional Definitions.
As used herein, ”this License” refers to version 3 of the GNU Lesser General Public License,
and the ”GNU GPL” refers to version 3 of the GNU General Public License.
”The Library” refers to a covered work governed by this License, other than an Application ora Combined Work as defined below.
An ”Application” is any work that makes use of an interface provided by the Library, but which
is not otherwise based on the Library. Defining a subclass of a class defined by the Library isdeemed a mode of using an interface provided by the Library.
A ”Combined Work” is a work produced by combining or linking an Application with the
Library. The particular version of the Library with which the Combined Work was made isalso called the ”Linked Version”.
The ”Minimal Corresponding Source” for a Combined Work means the Corresponding Source
for the Combined Work, excluding any source code for portions of the Combined Work that,considered in isolation, are based on the Application, and not on the Linked Version.
The ”Corresponding Application Code” for a Combined Work means the object code and/or
source code for the Application, including any data and utility programs needed for reproducing
the Combined Work from the Application, but excluding the System Libraries of the CombinedWork.
1. Exception to Section 3 of the GNU GPL.
You may convey a covered work under sections 3 and 4 of this License without being bound bysection 3 of the GNU GPL.
2. Conveying Modified Versions.
If you modify a copy of the Library, and, in your modifications, a facility refers to a function ordata to be supplied by an Application that uses the facility (other than as an argument passed
when the facility is invoked), then you may convey a copy of the modified version:a) under this License, provided that you make a good faith effort to ensure that, in the event
an Application does not supply the function or data, the facility still operates, and performs
whatever part of its purpose remains meaningful, orb) under the GNU GPL, with none of the additional permissions of this License applicable tothat copy.
3. Object Code Incorporating Material from Library Header Files.
87

The object code form of an Application may incorporate material from a header file that is partof the Library. You may convey such object code under terms of your choice, provided that,
if the incorporated material is not limited to numerical parameters, data structure layouts andaccessors, or small macros, inline functions and templates (ten or fewer lines in length), you doboth of the following:a) Give prominent notice with each copy of the object code that the Library is used in it andthat the Library and its use are covered by this License.
b) Accompany the object code with a copy of the GNU GPL and this license document.
4. Combined Works.
You may convey a Combined Work under terms of your choice that, taken together, effectively
do not restrict modification of the portions of the Library contained in the Combined Workand reverse engineering for debugging such modifications, if you also do each of the following:a) Give prominent notice with each copy of the Combined Work that the Library is used in it
and that the Library and its use are covered by this License.
b) Accompany the Combined Work with a copy of the GNU GPL and this license document.
c) For a Combined Work that displays copyright notices during execution, include the copyrightnotice for the Library among these notices, as well as a reference directing the user to the copiesof the GNU GPL and this license document.
d) Do one of the following:0) Convey the Minimal Corresponding Source under the terms of this License, and the Corresponding Application Code in a form suitable for, and under terms that permit, the user torecombine or relink the Application with a modified version of the Linked Version to produce a
modified Combined Work, in the manner specified by section 6 of the GNU GPL for conveying
Corresponding Source.
1) Use a suitable shared library mechanism for linking with the Library. A suitable mechanism
is one that (a) uses at run time a copy of the Library already present on the user’s computer
system, and (b) will operate properly with a modified version of the Library that is interfacecompatible with the Linked Version.
e) Provide Installation Information, but only if you would otherwise be required to provide suchinformation under section 6 of the GNU GPL, and only to the extent that such informationis necessary to install and execute a modified version of the Combined Work produced byrecombining or relinking the Application with a modified version of the Linked Version. (Ifyou use option 4d0, the Installation Information must accompany the Minimal Corresponding
Source and Corresponding Application Code. If you use option 4d1, you must provide theInstallation Information in the manner specified by section 6 of the GNU GPL for conveying
Corresponding Source.)
5. Combined Libraries.
You may place library facilities that are a work based on the Library side by side in a single
library together with other library facilities that are not Applications and are not covered bythis License, and convey such a combined library under terms of your choice, if you do both of
the following:a) Accompany the combined library with a copy of the same work based on the Library, uncombined with any other library facilities, conveyed under the terms of this License.
b) Give prominent notice with the combined library that part of it is a work based on the
Library, and explaining where to find the accompanying uncombined form of the same work.
6. Revised Versions of the GNU Lesser General Public License.
The Free Software Foundation may publish revised and/or new versions of the GNU Lesser
General Public License from time to time. Such new versions will be similar in spirit to thepresent version, but may differ in detail to address new problems or concerns.
Each version is given a distinguishing version number. If the Library as you received it specifies
that a certain numbered version of the GNU Lesser General Public License ”or any later version”
applies to it, you have the option of following the terms and conditions either of that published88version or of any later version published by the Free Software Foundation. If the Library as
you received it does not specify a version number of the GNU Lesser General Public License,
you may choose any version of the GNU Lesser General Public License ever published by theFree Software Foundation.
If the Library as you received it specifies that a proxy can decide whether future versions of the
GNU Lesser General Public License shall apply, that proxy’s public statement of acceptance ofany version is permanent authorization for you to choose that version for the Library.

89