SciPost Phys. 3, 013 (2017)

The ALF (Algorithms for Lattice Fermions) project release 1.0
Documentation for the auxiliary field quantum Monte Carlo code
Martin Bercx, Florian Goth, Johannes S. Hofmann and Fakher F. Assaad
Institut für Theoretische Physik und Astrophysik, Universität Würzburg,
97074 Würzburg, Germany
alf@physik.uni-wuerzburg.de

Abstract
The Algorithms for Lattice Fermions package provides a general code for the finite temperature auxiliary field quantum Monte Carlo algorithm. The code is engineered to beable to simulate any model that can be written in terms of sums of single-body operators, of squares of single-body operators and single-body operators coupled to an
Ising field with given dynamics. We provide predefined types that allow the user tospecify the model, the Bravais lattice as well as equal time and time displaced observables. The code supports an MPI implementation. Examples such as the Hubbard modelon the honeycomb lattice and the Hubbard model on the square lattice coupled to atransverse Ising field are provided and discussed in the documentation. We furthermore discuss how to use the package to implement the Kondo lattice model and theSU(N )-Hubbard-Heisenberg model. One can download the code from our Git instanceat https://alf.physik.uni-wuerzburg.de and sign in to file issues.
Copyright the ALF Project.
This work is licensed under the Creative Commons
Attribution-ShareAlike 4.0 International License.
Published by the SciPost Foundation.

Received 06-04-2017
Accepted 17-07-2017
Check forupdatesPublished 16-08-2017
doi:10.21468/SciPostPhys.3.2.013

You are free to share and benefit from this documentation as long as this license is preserved
and proper attribution to the authors is given. For details, see the ALF project homepage
alf.physik.uni-wuerzburg.de. Contact address: alf@physik.uni-wuerzburg.de.

Contents
1 Introduction
1.1 Motivation
1.2 Definition of the Hamiltonian
1.3 Outline

1
1
3
5

2 Auxiliary Field Quantum Monte Carlo
2.1 Formulation of the method
2.1.1 The partition function
2.1.2 Observables
2.1.3 Reweighting and the sign problem
2.2 Updating schemes
2.2.1 The default: sequential single spin flips

6
6
7
8
9
10
10

1

SciPost Phys. 3, 013 (2017)

2.2.2 Sampling of e−S0,I
Stabilization - a peculiarity of the BSS algorithm
Monte Carlo sampling
2.4.1 The Jackknife resampling method
2.4.2 An explicit example of error estimation
Pseudo code description

10
11
12
13
14
17

3 Data Structures and Input/Output
3.1 Implementation of the Hamiltonian and the lattice3.1.1 The Operator type
3.1.2 Specification of the model
3.1.3 The Lattice type
3.2 The observable types Obser_Vec and Obser_Latt
3.2.1 Scalar observables
3.2.2 Equal time and time displaced correlation functions
3.3 File structure
3.3.1 Input files
3.3.2 Output: Observables
3.3.3 Output: Precision
3.4 Scripts
3.5 Analysis programs
3.6 Running the code
3.6.1 Compilation
3.6.2 Starting a simulation
3.6.3 Error analysis

18
18
18
19
21
21
23
23
24
24
26
27
28
28
30
30
31
32

2.3
2.4

2.5

4 Examples
32
4.1 The SU(2)-Hubbard model on a square lattice32
4.1.1 Setting the Hamiltonian: Ham_set
32
4.1.2 Observables
34
4.1.3 Numerical precision
36
4.2 The Mz -Hubbard model on a square lattice37
4.2.1 The interaction term: Call Ham_V37
4.2.2 The measurements: Call Obser, Call ObserT
38
4.2.3 Numerical precision38
4.3 The SU(2)-Hubbard model on the honeycomb lattice38
4.3.1 Working with multi-orbital unit cells: Call Ham_Latt38
4.3.2 The hopping term: Call Ham_Hop
39
4.3.3 Observables: Call Obser, Call ObserT
39
4.4 The SU(2)-Hubbard model on a square lattice coupled to a transverse Ising field 40
4.4.1 The Ising term
40
4.4.2 The interaction term: Call Ham_V
40
4.4.3 The function Real (Kind=8) function S0(n,nt)
41
5 Miscellaneous
5.1 Other models
5.1.1 Kondo lattice model
5.1.2 SU(N )-Hubbard-Heisenberg models
5.2 Performance, memory requirements and parallelization

42
42
42
42
43

6 Conclusions and Future Releases

442

SciPost Phys. 3, 013 (2017)

A License

45

References

46

1

Introduction

1.1

Motivation

The auxiliary field quantum Monte Carlo (QMC) approach is the algorithm of choice to simulate thermodynamic properties of a variety of correlated electron systems in the solid stateand beyond [1–8]. Apart from the physics of the canonical Hubbard model [9, 10], the topics
one can investigate in detail include correlation effects in the bulk and on surfaces of topological insulators [11, 12], quantum phase transitions between Dirac fermions and insulators[13–17], deconfined quantum critical points [18,19], topologically ordered phases [19], heavy
fermion systems [20, 21], nematic [22] and magnetic [23] quantum phase transitions in metals, antiferromagnetism in metals [24], superconductivity in spin-orbit split bands [25], SU(N )symmetric models [26, 27], long-ranged Coulomb interactions in graphene systems [28, 29],
cold atomic gasses [30], low energy nuclear physics [31], entanglement entropies and spectra [32–36], etc. This ever growing list of topics is based on algorithmic progress and on recent
symmetry related insights [37–41] enabling one to find negative sign problem free formulations of a number of model systems with very rich phase diagrams.
Auxiliary field methods can be formulated in very different ways. The fields define the configuration space C. They can stem from the Hubbard-Stratonovich (HS) [42] transformation
required to decouple the many-body interacting term into a sum of non-interacting problems,or they can correspond to bosonic modes with predefined dynamics such as phonons or gauge
fields. In all cases, the result is that the grand-canonical partition function takes the form, Xˆ
Z = Tr e−β H =e−S(C) ,(1)
C

where S is the action of non-interacting fermions subject to a space-time fluctuating auxiliary
field. The high-dimensional integration over the fields is carried out stochastically. In thisformulation of many body quantum systems, there is no reason for the action to be a real
number. Thereby e−S(C) cannot be interpreted as a weight. To circumvent this problem one
can adopt re-weighting schemes and sample |e−S(C) |. This invariably leads to the so called
negative sign problem with associated exponential computational scaling in system size andinverse temperature [43, 44]. The sign problem is formulation dependent, and as mentioned
above there has been tremendous progress at identifying an ever growing class of negative
sign problem free models covering a rich domain of collective emergent phenomena. For
continuous fields, the stochastic integrations can be carried out with Langevin dynamics or
hybrid methods [45]. However, for many problems one can get away with discrete fields [46].
In this case, Monte Carlo importance sampling will often be put to use [47]. We note that
due to the non-locality of the fermion determinant, see below, cluster updates, such as in theloop or stochastic series expansion algorithms for quantum spin systems [48–50], are hard toformulate for this class of problems. The search for efficient updating schemes that enable tomove quickly within the configuration space defines ongoing challenges.
Formulations do not differ only by the choice of the fields, continuous or discrete, and thesampling strategy, but also by the formulation of the action itself. For a given field configura3

SciPost Phys. 3, 013 (2017)

tion, integrating out fermionic degrees of freedom generically leads to a fermionic determinantof dimension β N where β corresponds to the inverse temperature and N to the volume of thesystem. Working with this determinant leads to the Hirsch-Fye approach [51] and its timecomplexity which quantifies the computational effort is given by O (β N )3 . 1 The Hirsch-Fyealgorithm is the method of choice for impurity problems, but has generically been outperformed by a class of so-called continuous-time quantum Monte Carlo approaches [52–54].
One key point of continuous-time methods is that they are action based and thereby allow tohandle retarded interactions obtained when integrating out fermion or boson baths. In highdimensions and/or at low temperatures, the cubic scaling originating from the fermionic determinant is expensive. To circumvent this, the hybrid Monte-Carlo approach [5,55] expresses
the fermionic determinant in terms of a Gaussian integral thereby introducing a new variablein the Monte Carlo integration. The resulting algorithm is the method of choice for latticegauge theories in 3+1 dimensions and has been used to provide ab-inito estimates of light
hadron masses starting from quantum chromo dynamics [56].
The algorithm implemented in the ALF project lies between the above two extremes. We
will keep the fermionic determinant, but formulate the problem so as to work only with N × N
matrices. This Blankenbecler, Scalapino, Sugar (BSS) algorithm scales linearly in imaginary
time β, but remains cubic in the volume N . Furthermore, the algorithm can be formulated
either in a projective manner [3, 4], adequate to obtain zero temperature properties in the
canonical ensemble, or at finite temperatures in the grand-canonical ensemble [2].
The aim of the ALF project is to introduce a general formulation of the finite temperature
auxiliary field QMC method with discrete fields so as to quickly be able to play with different
model Hamiltonians at minimal programming cost. We have summarized the essential aspectsof the auxiliary field QMC approach in this documentation, and refer the reader to Refs. [7,57]
for complete reviews. We will show in all details how to implement a variety of models, run
the code, and produce results for equal time and time displaced correlation functions. The
program code is written in Fortran according to the 2003 standard and is able to natively
utilize MPI for massively parallel runs on todays supercomputing systems.
The ALF package is not the first open source project aimed at providing simulation tools
for correlated quantum matter. The most notable package is certainly the ALPS library [58]. It
is actively maintained and features a whole set of algorithms for strongly correlated quantum
lattice models including Monte Carlo, exact diagonalization, and density matrix renormalization group codes. It however does not include the auxiliary field QMC algorithm offered by
the ALF package. Other projects include QUEST [59], TRIQS [60], w2dynamics [61] andiQist [62]. IQist, TRIQS and w2dynamics focus on approximate solutions via the CT-HYB [52]algorithm within the dynamical mean field approximation. The QUEST project implements
the same algorithm as in the ALF project but is currently restricted to the Hubbard model andit does not allow to easily incorporate different Hamiltonians.
The ALF source code is placed under the GNU GPL license. The project is currently hosted
on servers of the university of Würzburg where we have set up a GitLab instance (https://alf.
physik.uni-wuerzburg.de) aimed at encouraging community outreach. Each potential user cansign in, receive space for his ALF related projects and share them with others. This site serves
the GitLab issue tracker as well as a wiki so that members can collect information they consider
useful for the project. We have set up an E-Mail address for reaching the core developers at
alf@physik.uni-wuerzburg.de.
1Here we implicitly assume the absence of a negative sign problem.

4

SciPost Phys. 3, 013 (2017)

1.2

Definition of the Hamiltonian

The first and most fundamental part of the project is to define a general Hamiltonian which
can accommodate a large class of models. Our approach is to express the model as a sum of
one-body terms, a sum of two-body terms each written as a perfect square of a one body term,
as well as a one-body term coupled to an Ising field with dynamics to be specified by the user.
The form of the interaction in terms of sums of perfect squares allows us to use generic formsof discrete approximations to the HS transformation [63, 64]. Symmetry considerations areimperative to enhance the speed of the code. We therefore include a color index reflecting an
underlying SU(N ) color symmetry as well as a flavor index reflecting the fact that after the HStransformation, the fermionic determinant is block diagonal in this index.
The class of solvable models includes Hamiltonians Hˆ that have the following generalform:
Hˆ = HˆT + HˆV + HˆI + Hˆ0,I ,HˆT =Ncol X
MT X
Nfl N
dim
X
X
k=1 σ=1 s=1 x, y

HˆV =
HˆI =

MV
X
k=1
MI
X
k=1

Uk

ĉ †xσs Tx(ks)y ĉ yσs ≡N N N
col X
dim
fl X
X
σ=1 s=1 x, y

(2)

MT
X

T̂ (k) ,(3)

k=1¨ N N N
col X
dim
fl
X
X
σ=1 s=1

Ẑk

where

ĉ †xσs Vx(ks)y ĉ yσs

«2


+ αks

x, y

≡

MV
X

Uk V̂ (k)

2

,

(4)

k=1

ĉ †xσs I x(ks)y ĉ yσs


≡

MI
X

Ẑk Î (k) .

(5)

k=1

The indices and symbols have the following meaning:• The number of fermion flavors is set by Nfl . After the HS transformation, the action willbe block diagonal in the flavor index.
• The number of fermion colors is set by Ncol . The Hamiltonian is invariant under SU(Ncol )rotations. 2• Both the color and the flavor index can describe the spin degree of freedom, the choicedepending on the spin symmetry of the simulated model and the HS transformation.
This point is illustrated in the examples, see Secs. 4.1 and 4.2.
• Ndim is the total number of spacial vertices: Ndim = Nunit cell Norbital , where Nunit cell isthe number of unit cells of the underlying Bravais lattice and Norbital is the number of(spacial) orbitals per unit cell.
• The indices x and y label lattice sites where x, y = 1, · · · , Ndim .
• Therefore, the matrices T (ks) , V (ks) and I (ks) are of dimension Ndim × Ndim .
• The number of interaction terms is labelled by MV and M I . M T > 1 would allow for acheckerboard decomposition.
• ĉ †yσs is a second quantized operator that creates an electron in a Wannier state centeredaround lattice site y, with color σ, and flavor index s. The operators satisfy the anticommutation relations:¦©
¦©
ĉ †yσs , ĉ y 0 σ0 s0 = δ y, y 0 δs,s0 δσ,σ0 , and ĉ yσs , ĉ y 0 σ0 s0 = 0.
(6)2

Note that in the code Ncol ≡ N_SUN.

5

SciPost Phys. 3, 013 (2017)

The Ising part of the general Hamiltonian (2) is Hˆ0,I + HˆI and has the following properties:• Ẑk is an Ising spin operator which corresponds to the Pauli matrix σ̂z . It couples to a
general one-body term.
• The dynamics of the Ising spins is given by Hˆ0,I . This term is not specified here; it has tobe specified by the user and becomes relevant when the Monte Carlo update probability
is computed in the code (see Sec. 4.4 for an example).
Note that the matrices T (ks) , V (ks) and I (ks) explicitly depend on the flavor index s but not on
the color index σ. The color index σ only appears in the second quantized operators such thatthe Hamiltonian is manifestly SU(Ncol ) symmetric. We also require the matrices T (ks) , V (ks)and I (ks) to be Hermitian.
As we will detail below, the definition of the above Hamiltonian allows to tackle several
non-trivial models and phenomena. There are however a number of model Hamiltonians thatcannot be simulated with ALF. Since we have opted for discrete fields, the electron-phonon
interaction is not included. Furthermore, continuous HS transformations, that turn out to be
extremely useful to include long-range Coulomb interactions [28, 65, 66], are not accessible
in the present form of the package.3 In many cases such as in 3 He, three - or more bodyinteractions should be included to capture relevant exchange mechanisms [67, 68]. Thesehigher order processes are not captured in the ALF since it is limited to two-body interactions.
The formulation of the Hamiltonian, has an explicit global U(1) symmetry corresponding toparticle number conservation. Hence using the ALF for a given model implies the existence ofa canonical transformation where a particle number is conserved. Imaginary time dependent
Hamiltonians, required to compute Renyi entropies and entanglement spectra [33, 35, 36]
are not yet in the scope of ALF. Finally, one should also mention that auxiliary field QMC
simulations are Hamiltonian based such that retarded interactions are not included in the ALF.
For this set of problems, CT-INT type approaches are the method of choice [53, 54, 69]. The
above short comings partially define a set of future directions that will be discussed in the
concluding part of this documentation.

1.3

OutlineTo use the code, a minimal understanding of the algorithm is necessary. In Sec. 2, we go very
briefly through the steps required to formulate the many-body imaginary-time propagation interms of a sum over HS and Ising fields of one-body imaginary-time propagators. The user has
to provide this one-body imaginary-time propagator for a given configuration of HS and Ising
fields. We equally discuss the Monte Carlo updates, the strategies for numerical stabilizationof the code, as well as the Monte Carlo sampling.
Section 3 is devoted to the data structures that are needed to implement the model, as
well as to the input and output file structure. The data structure includes an Operator typeto optimally work with sparse Hermitian matrices, a Lattice type to define one- and twodimensional Bravais lattices, and two Observable types to handle scalar observables (e.g.
total energy) and equal time or time displaced two- point correlation functions (e.g. spin-spincorrelations).
The Monte Carlo run and the data analysis are separated: the QMC run dumps the resultsof bins sequentially into files which are then analyzed by analysis programs. In Sec. 3.5, we
provide a brief description of the analysis programs for our observable types. The analysis
programs allow for omitting a given number of initial bins in order to account for warmup
times. Also, a rebinning analysis is included to a posteriori take account of long autocorrelation
times. Finally, Sec. 3.6 provides all the necessary details to compile and run the code.
3

Note however that one can readily add short ranged interactions by including terms such as (n̂i + n̂ j − 2)2 .

6

SciPost Phys. 3, 013 (2017)

In Sec. 4, we give explicit examples on how to use the code for the Hubbard model onsquare and honeycomb lattices, for different choices of the Hubbard-Stratonovich transformation (see Secs. 4.1, 4.2 and 4.3) as well as for the Hubbard model on a square lattice coupled
to a transverse Ising field (see Sec. 4.4 ). Our implementation is rather general such that avariety of other models can be simulated. In Sec. 5 we provide some information on how tosimulate the Kondo lattice model as well as the SU(N ) symmetric Hubbard-Heisenberg model.
Finally, in Sec. 6 we list a number of features that are considered for future releases of the
ALF program package.

2

Auxiliary Field Quantum Monte Carlo

2.1

Formulation of the method

Our aim is to compute observables for the general Hamiltonian (2) in thermodynamic equilibrium as described by the grand-canonical ensemble. We will show below how the grandcanonical partition function is rewritten as Xˆ
Z = Tr e−β H =e−S(C) + O (∆τ2 )(7)
Cand define the space of configurations C. Note that the chemical potential term is alreadyincluded in the definition of the one-body term HˆT , see eq. (3), of the general Hamiltonian.
The outline of this section is as follows. First, we derive the detailed form of the partition
function and outline the computation of observables (Sec. 2.1.1 - 2.1.3). Next, we present thepresent update strategy, namely local updates (Sec. 2.2). We equally discuss the measures wehave implemented to make the code numerically stable (Sec. 2.3). Finally, we discuss the autocorrelations and associated time scales during the Monte Carlo sampling process (Sec. 2.4).
The essential ingredients of the auxiliary field quantum Monte Carlo implementation in
the ALF package are the following:• We will discretize the imaginary time propagation: β = ∆τLTrotter . Generically thisintroduces a systematic Trotter error of O (∆τ)2 [70]. We note that there has been considerable effort at getting rid of the Trotter systematic error and to formulate a genuine
continuous-time BSS algorithm [71]. To date, efforts in this direction are based on a
CT-AUX type formulation [72, 73] and face two issues. The first issue is that they are
restricted to a class of models with Hubbard-type interactions(n̂i − 1)2 = (n̂i − 1)4 ,(8)such that the basic CT-AUX equation [74]1+U
1 X αs(n̂i −1)U(n̂i − 1)2 =e
with= cosh(α) − 1 and K ∈ R
K
2 s=±1K(9)holds. The second issue is that in the continuous-time approach it is hard to formulate acomputationally efficient algorithm. Given this situation it turns out that the multi-grid
method [75–77] is an efficient scheme to extrapolate to small imaginary-time steps soas to eliminate the Trotter systematic error if required.
• Having isolated the two-body term, we will use the discrete HS transformation [63, 64]:
X
p
2
(10)
e∆τλÂ =γ(l)e ∆τλη(l)Â + O (∆τ4 ) ,l=±1,±2

7

SciPost Phys. 3, 013 (2017)

where the fields η and γ take the values:r p γ(±1) = 1 + 6/3, η(±1) = ± 2 3 − 6 ,
r 
p
p 
γ(±2) = 1 − 6/3, η(±2) = ± 2 3 + 6 .
p(11)

Since the Trotter error is already of order (∆τ2 ) per time slice, this transformation isnext to exact.
• We will work in a basis for the Ising spins where Ẑk is diagonal: Ẑk |sk 〉 = sk |sk 〉, where
sk = ±1.
• From the above it follows that the Monte Carlo configuration space C is given by the
combined spaces of Ising spin configurations and of HS discrete field configurations:C = si,τ , l j,τ with i = 1 · · · M I , j = 1 · · · MV , τ = 1 · · · LTrotter .
(12)Here, the Ising spins take the values si,τ = ±1 and the HS fields take the values
l j,τ = ±2, ±1.
2.1.1

The partition function

With the above, the partition function of the model (2) can be written as follows.
ˆZ = Tr e−β H LTrotter

MV
MI
MT
Y
Y
Y
(k)
(k)
(k) 2
ˆ
+ O (∆τ2 )= Tr e−∆τH0,Ie−∆τUk (V̂ )e−∆τσ̂k Î
e−∆τ T̂
k=1

=

M L
Trotter
V Y
X Y
C

k=1

k=1



γk,τ e−S0,I ({si,τ }) ×

k=1 τ=1

¨L
M
«
MT
MI
Trotter Y
V
Y
Y
Y
p
(k)
(k)
(k)
TrF
e −∆τUk ηk,τ V̂
+ O (∆τ2 ) .
e−∆τ T̂e−∆τsk,τ Î
τ=1

k=1

(13)

k=1

k=1

In the above, the trace Tr runs over the Ising spins as well as over
 the fermionic degrees of
freedom, and TrF only over the fermionic Fock space. S0,I si,τ is the action corresponding
to the Ising Hamiltonian, and is only dependent on the Ising spins so that it can be pulled out of
the fermionic trace. We have adopted the short hand notation ηk,τ = η(l k,τ ) and γk,τ = γ(l k,τ ).
At this point, and since for a given configuration C we are dealing with a free propagation, we
can integrate out the fermions to obtain a determinant:
¨L
M
«
MI
MT
Trotter Y
V
Y
Y
Y
p
(k)
(k)
(k)
e−∆τsk,τ Îe−∆τ T̂=
TrF
e −∆τUk ηk,τ V̂
τ=1

k=1

k=1
Nfl
Y



k=1
M
PV LTrotterPe k=1 τ=1

Ncol

p

−∆τUk αk,s ηk,τ



×

s=1
Nfl
Y
s=1





det 1 +LY
MV
Trotter Y

ep−∆τUk ηk,τ V

τ=1 k=1

(ks)

MI
Y
k=1

8

e−∆τsk,τ I

(ks)

MT
Y
k=1

e−∆τT

(ks)

Ncol, (14)

SciPost Phys. 3, 013 (2017)

where the matrices T(ks) , V(ks) , and I(ks) define the Hamiltonian [Eq. (2) - (5)]. All in all, thepartition function is given by:
Z =

X

e

−S0,I ({si,τ })

≡

s=1
X


γk,τ eNcol

Nfl MV LTrotter pP
P P−∆τUk αk,s ηk,τs=1 k=1

τ=1×

k=1 τ=1

C
Nfl
Y

M L
Trotter
V Y
Y





det 1 +LY
MV
Trotter Y

ep−∆τUk ηk,τ V (ks)

τ=1 k=1

e

−S(C)

MI
Y

e

−∆τsk,τ I (ks)

k=1

MT
Y

e

−∆τT (ks)Ncol+ O (∆τ2 )k=1

2+ O (∆τ ) .

(15)

C

In the above, one notices that the weight factorizes in the flavor index. The color index raises
the determinant to the power Ncol . This corresponds to an explicit SU(Ncol ) symmetry for eachconfiguration. This symmetry is manifest in the fact that the single particle Green functions
are color independent, again for each given configuration C.
2.1.2

Observables

In the auxiliary field QMC approach, the single-particle Green function plays a crucial role. It
determines the Monte Carlo dynamics and is used to compute observables:X
Tr e−β Ĥ Ô
e−S(C)
 =〈Ô〉 =
P(C)〈〈Ô〉〉(C) , with P(C) = P
.
(16)
−S(C)
Tr e−β Ĥ
CeC〈〈Ô〉〉(C) corresponds to the expectation value of Ô for a given configuration C. For a given
configuration C one can use Wick’s theorem to compute 〈〈Ô〉〉(C) from the knowledge of thesingle-particle Green function:
G(x, σ, s, τ|x 0 , σ0 , s0 , τ0 ) = 〈〈T ĉ xσs (τ)ĉ †x 0 σ0 s0 (τ0 )〉〉C ,(17)

where T corresponds to the imaginary-time ordering operator. The corresponding equal timequantity reads:
G(x, σ, s, τ|x 0 , σ0 , s0 , τ) = 〈〈ĉ xσs (τ)ĉ †x 0 σ0 s0 (τ)〉〉C .
(18)
Since for a given HS field translation invariance in imaginary-time is broken, the Green function has an explicit τ and τ0 dependence. On the other hand it is diagonal in the flavor index,
and independent on the color index. The latter reflects the explicit SU(N ) color symmetry
present at the level of individual HS configurations. As an example, one can show that theequal time Green function at τ = 0 reads [7]:G(x, σ, s, 0|x 0 , σ, s, 0) = 1 +LY
Trotter
τ=1

withBτ(s) =

MT
Y
k=1

e

−∆τT (ks)

Bτ(s)

−1
(19)
x,x 0

MV p
MI
Y
Y
(ks)
−∆τUk ηk,τ V (ks)ee−∆τsk,τ I .

(20)

k=1

k=1

To compute equal time as well as time displaced observables, one can make use of Wick’s

9

SciPost Phys. 3, 013 (2017)

theorem. A convenient formulation of this theorem for QMC simulations reads:〈〈T c †x (τ1 )c x 0 (τ01 ) · · · c †x (τn )c x 0 (τ0n )〉〉C =
1
n
n
1

†
0
〈〈T c x (τ1 )c x 0 (τ1 )〉〉C 〈〈T c †x (τ1 )c x 0 (τ02 )〉〉C
1
1
1
2

0
†
0
†
〈〈T c x 2 (τ2 )c x 0 (τ1 )〉〉C 〈〈T c x 2 (τ2 )c x 0 (τ2 )〉〉C
1
2

det 
..
..

.
.
〈〈T c †x (τn )c x 0 (τ01 )〉〉C 〈〈T c †x (τn )c x 0 (τ02 )〉〉C
n

n

1

...

〈〈T c †x (τ1 )c x 0 (τ0n )〉〉C...

〈〈T c †x (τ2 )c x 0 (τ0n )〉〉C 2n . (21)
..


.
†
0
〈〈T c x (τn )c x 0 (τn )〉〉C

..

.
...

2

1

n



n

n

Here, we have defined the super-index x = {x, σ, s}. In the subroutines Obser and ObserT of
the module Hamiltonian_Examples.f90 (see Sec. 3.2) the user is provided with the equal
time and time displaced correlation function. Using the above formulation of Wick’s theorem,
arbitrary correlation functions can be computed. We note however, that the program is limited
to the calculation of observables that contain only two different imaginary times.
2.1.3

Reweighting and the sign problem

In general, the action S(C) will be complex, thereby inhibiting a direct Monte Carlo samplingof P(C). This leads to the infamous sign problem. The sign problem is formulation dependentand as noted above, much progress has been made at understanding the class of models thatcan be formulated without encountering this problem [37–40]. When the average sign is not
too small, we can nevertheless compute observables within a reweighting scheme. Here we
adopt the following scheme. First note that the partition function is real such that:
X
X
X 

Z=
e−S(C) =e−S(C) =ℜ e−S(C) .
(22)
C

C

C

Thereby4 , and with the definition
ℜ e−S(C) ,sign (C) =ℜ e−S(C)the computation of the observable [Eq. (16)] is re-expressed as follows:
P −S(C)〈〈Ô〉〉(C)CeP〈Ô〉 =−S(C)Ce −S(C)  e−S(C)P〈〈Ô〉〉(C)Cℜ eℜ[e−S(C) ]


=
P
−S(C)
Cℜ e
nP
o P
 −S(C)  −S(C) e−S(C)
ℜ
e
〈〈
Ô〉〉
sign
(C)
(C) /
C
C ℜ e
ℜ[e−S(C) ]


P


=
P
−S(C)
sign (C) / C ℜ e−S(C)
C ℜ e
¬
¶
e−S
sign ℜ[e−S ] 〈〈Ô〉〉P=.
〈sign〉 P

(23)

(24)

The average sign isP〈sign〉 P = −S(C) ℜesign (C)
C


,
P
−S(C)
C ℜ e

(25)

4
The attentive reader will have noticed that for arbitrary Trotter decompositions, the imaginary time propagatoris not necessarily Hermitian. Thereby, the above equation is correct only up to corrections stemming from thecontrolled Trotter systematic error.

10

SciPost Phys. 3, 013 (2017)

and we have 〈sign〉 P ∈ R per definition. The Monte Carlo simulation samples the probability
distributionℜ e−S(C) .
P(C) = P
(26)
−S(C)
C ℜ e
such that the nominator and denominator of Eq. (24) can be computed.
The negative sign problem is an issue since the average sign is a ratio of two partition
functions such that one can argue that [43]〈sign〉 P ∝ e−∆N β .

(27)∆ is intensive positive quantity andp N β denotes the Euclidean volume. In a Monte Carlosimulation, the error scales as 1/ TCPU where TCPU corresponds to the computational time.
Since the error on the average sign has to be much smaller than the average sign itself, one
sees that:
TCPU  e2∆N β .
(28)
Two comments are in order. First, the presence of a sign problem invariably leads to an exponential increase of CPU time as a function of the Euclidean volume. And second, ∆ isformulation dependent. For instance, at finite doping, the SU(2) invariant formulation of the
Hubbard model presented in Sec. 4.1 has a much more severe sign problem than the formulation presented in Sec. 4.2 where the HS field couples to the z-component of the magnetization.
Typically one can work with average signs down to 〈sign〉 P ' 0.1.

2.2

Updating schemes

The program allows for different types of updating schemes. Given a configuration C wepropose a new one, C 0 , with probability T0 (C → C 0 ) and accept it according to the MetropolisHastings acceptance-rejection probability,
T0 (C 0 → C)W (C 0 )0P(C → C ) = min 1,,(29)
T0 (C → C 0 )W (C)
so as to guarantee the stationarity condition. Here, W (C) = ℜ e−S(C) .
VariableType

Description

Propose_S0

LogicalIf true, proposes local moves according to the probability e−S0,I .

Table 1: Variable required to control the updating scheme.

2.2.1

The default: sequential single spin flips

The default updating scheme is a sequential single spin flip algorithm. Consider the Ising spinsi,τ . We will flip it with probability one such that for this local move the proposal matrix issymmetric. If we are considering the Hubbard-Stratonovich field l i,τ we will propose withprobability 1/3 one of the other three possible fields. Again, for this local move, the proposal
matrix is symmetric. Hence in both cases we will accept or reject the move according to
W (C 0 )0.
(30)
P(C → C ) = min 1,
W (C)It is worth noting that this type of sequential spin flip updating does not satisfy detailed balance
but the more fundamental stationarity condition [47].
11

SciPost Phys. 3, 013 (2017)

2.2.2

Sampling of e−S0,I

Consider an Ising spin at space-time i, τ and the configuration C. Flipping this spin will generate the configuration C 0 and we will propose the move according to0

e−S0,I (C )0

T0 (C → C ) =e−S0,I (C 0 )+e−S0,I (C)

=1−

1
1+e

−S0,I (C 0 )/e−S0,I (C)

.

(31)

Note that the function S0 in the Hamitonian_example.f90 module computes precisely the0
ratio e−S0,I (C ) /e−S0,I (C) so that T0 (C → C 0 ) does not require any further programming. Thereby
one will accept the proposed move with the probability:
e−S0,I (C) W (C 0 )0P(C → C ) = min 1, −S (C 0 ).
(32)
e 0,I W (C)With Eq. 15 one sees that the bare action S0,I (C) determining the dynamics of the Ising spinin the absence of coupling to the fermions does not enter the Metropolis acceptance-rejection
step. This sampling scheme is used if the logical variable Propose_S0 is set to true.

2.3

Stabilization - a peculiarity of the BSS algorithm

From (15) it can be seen that for the calculation of the Monte Carlo weight and for the observables a long product of matrix exponentials has to be formed. On top of that we need to
be able to extract the single-particle Green function for a given flavor index at say time slice
τ = 0. As mentioned above in Eq. (19), this quantity is given by:−1
LY
Trotter
G = 1+
Bτ
.
(33)
τ=1

To boil this down to more familiar terms from linear algebra we remark that we can recast thisproblem as the task to find the solution of the linear systemY(1 +Bτ )x = b.
(34)
τ

The Bτ ∈ Cn×n depend on the lattice size as well as other physical parameters that can bechosen such that a matrix norm of Bτ can be unbound in size. From standard perturbation
theory for linear systems it is known that the computed solution x̃ would contain a relative
error ofY |x̃ − x|Bτ.
(35)
= O εκ p 1 +|x|τHere ε denotes the machine precision, which is 2−53 for IEEE double precision numbers, andκ p (M) is the condition number of the matrix M with respect to the matrix p-norm. The
importantproperty that makes straight-forward inversion so badly suited stems from the fact
Q
that τ Bτ contains exponentially large and small scales as can be seen in Eq. (15). Thereby, as
a function of increasing inverse temperature, the condition number will grow exponentially sothat the computed solution x̃ will often contain no correct digits at all. To circumvent this, moresophisticated methods have to be employed. We will first of all assume that the multiplicationof NWrap B matrices has an acceptable condition number. Assuming for simplicity that LTrotter
is an integer multiple of NWrap, we can write:−1LTrotterNWrap
Y−1 NWrapY
G = 1 +Bi·NWrap+τ 
i=0|τ=1 {z
}

≡Bi

12

.

(36)

SciPost Phys. 3, 013 (2017)

Within the auxiliary field QMC implementation of the ALF project, we are by default employing the strategy of forming a product of QR-decompositions which was proven to be weakly
backwards stable in Ref. [78]. The key idea is to efficiently separate the scales of a matrix fromthe orthogonal part of a matrix. This can be achieved using a QR decomposition of a matrix
A in the form Ai = Q i Ri . The matrix Q i is unitary and hence in the usual 2-norm it holds that
κ2 (Q i ) = 1. To get a handle on the condition number of Ri we will form the diagonal matrix(Di )n,n = |(Ri )n,n | ,(37)

and set R̃i = Di−1 Ri This gives the decomposition
Ai = Q i Di R̃i .

(38)

Di now contains the row norms of the original Ri matrix and hence attempts to separate off thetotal scales of the problem from Ri . This is similar in spirit to the so-called matrix equilibration
which tries to improve the condition number of a matrix by suitably chosen column and row
scalings. Due to a theorem by van der Sluis [79] we know that the choice in Eq. (37) is almost
optimal among all diagonal matrices D from the space of diagonal matrices D in the sense that
κ p ((Di )−1 Ri ) ≤ n1/p min κ p (D −1 Ri ).
D∈DNow, given an initial decomposition of A j−1 =formed in the following three steps:

Q

i Bi = Q j−1 D j−1 T j−1 an update B j A j−1 is

1. Form M j = (B j Q j−1 )D j−1 . Note the parentheses.
2. Do a QR decomposition of M j = Q j D j R j . This gives the final Q j and D j .
3. Form the updated T matrices T j = R j T j−1 .
While this might seem like quite an effort that has to be performed for every multiplication
it has to be noted that even with this stabilization scheme the algorithm preserves the timecomplexity class of O (β N 3 ) expressed in the physical parameters inverse temperature β andlattice size N . While there is no analytical expression for the dependence of the stability on
the physical parameters our experience has been that for a given number of stabilization steps
along the imaginary time axis [in the notation of Eq. (36) this number is LTrotter /NWrap],
the precision will be largely invariant of the system size N , whereas with increasing inverse
temperature β the number of stabilization steps often has to be increased to maintain a given
precision. The effectiveness of the stabilization has to be judged for every simulation from theoutput file info (Sec. 3.3.3). For most simulations there are two values to look out for:• Precision Green• Precision Phase
The Green function as well as the average phase are usually numbers with a magnitude ofO (1). For that reason we recommend that NWrap is chosen such that the mean precision is of
the order of 10−8 or better. We have included typical values of Precision Phase and of themean and the maximal values of Precision Green in the discussion of example simulations,
see Sec. 4.1.3 and Sec. 4.2.3.

13

SciPost Phys. 3, 013 (2017)

2.4

Monte Carlo sampling

Error estimates in Monte Carlo simulations can be delicate and are based on the central limit
theorem [80]. This theorem requires independent measurements and a finite variance. In thissubsection we will give examples of the issues that a user will have to look out for while usinga Monte Carlo code. Those effects are part of the common lore of the field and we can only
touch on them briefly in this text. For a deeper understanding of the inherent issues of Markov
chain Monte Carlo methods we refer the reader to the pedagogical introduction in chapter
1.3.5 of Krauth [81], the overview article of Sokal [47], the more specialized literature byGeyer [82] and chapter 6.3 of Neal [83].
In general, one distinguishes local from global updates. As the name suggest, the local
update corresponds to a small change of the configuration, e.g. a single spin flip of one of the
LTrotter (M I + MV ) field entries (see Sec. 2.2), whereas a global update changes a significant
part of the configuration. The default update scheme of the implementation at hand are local
updates such that a minimum amount of moves is required to generate a independent configuration. The associated time scale is called the autocorrelation time, Tauto , and is genericallydependent upon the choice of the observables.
Our unit of sweeps is defined such that each field is visited twice in a sequential propagationfrom τ = 0 to τ = L Trotter and back. A single sweep will generically not suffice to produce an
independent configuration. In fact, the autocorrelation time Tauto characterizes the required
time scale to generate an independent values of 〈〈Ô〉〉C for the observable O. This has several
consequences for the Monte Carlo simulation:• First of all, we start from a randomly chosen field configuration such that one has toinvest at least one, but generically much more, Tauto to generate relevant, equilibrated
configurations before reliable measurements are possible. This phase of the simulation
is known as the warm-up or burn-in phase. In order to keep the code as flexible as possible (different simulations might have different autocorrelation times), measurements
are taken from the very beginning. Instead, we provide the parameter n_skip for the
analysis to ignore the first n_skip bins.
• Secondly, our implementation averages over a given amount of measurements set by thevariable NSWEEPS before storing the results, known as one bin, on the disk. A bin corresponds to NSWEEPS sweeps. The error analysis requires statistically independent bins
to generate reliable confidence estimates. If bins are to small (averaged over a period
shorter then Tauto ), the error bars are then typically underestimated. Most of the time,
the autocorrelation time is unknown before the simulation is started. Sometimes theused compute cluster does not allow single runs long enough to generate appropriately
sized bins. Therefore, we provide the N_rebin parameter that specifies how many bins
are combined into a new bin during the error analysis. In general, one should checkthat a further increase of the bin size does not change the error estimate (For an explicit
example, the reader is referred to Sec. 2.4.2 and the appendix of Ref. [57]).
The N_rebin variable can be used to control a second issue. The distribution of theMonte Carlo estimates 〈〈Ô〉〉C is unknown. The result in the form (mean ± error) assumes a Gaussian distribution. Every original distribution with a finite variance turns
into a Gaussian one, once it is folded often enough (central limit theorem). Due to theinternal averaging (folding) within one bin, many observables are already quite Gaussian. Otherwise one can increase N_rebin further, even if the bins are already independent [84].
• The third issue concerns time displaced correlation functions. Even if the configurations
are independent, the fields within the configuration are still correlated. Hence, the data
14

SciPost Phys. 3, 013 (2017)

for Sα,β (~k, τ) (see Sec. 3.2; Eqn. 65) and Sα,β (~k, τ + ∆τ) are also correlated. Setting
the switch N_Cov=1 triggers the calculation of the covariance matrix in addition to theusual error analysis. The covariance is defined by
COVττ0 =

1
NBin



Sα,β (~k, τ) − 〈Sα,β (~k, τ)〉 Sα,β (~k, τ0 ) − 〈Sα,β (~k, τ0 )〉 .

(39)

An example where this information is necessary is the calculation of mass gaps extracted
by fitting the tail of the time displaced correlation function. Omitting the covariancematrix will underestimate the error.
2.4.1

The Jackknife resampling method

For each observable Â, B̂, Ĉ · · · the Monte Carlo program computes a data set of NBin (ideally)
independent values where for each observable the measurements belong to the same statistical distribution. In the general case, we would like to evaluate a function of expectation
values, f (〈Â〉, 〈B̂〉, 〈Ĉ〉 · · · ) – see for example the expression (24) for the observable including
reweighting – and are interested in the statistical estimates of its mean value and the standard
error of the mean. A numerical method for the statistical analysis of a given function f which
properly handles error propagation and correlations among the observables is the Jackknife
method, which is, like the related Bootstrap method, a resampling scheme [85]. Here webriefly review the delete-1 Jackknife scheme which is based on the idea to generate Nbin new
data sets of size Nbin −1 by consecutively removing one data value from the original set. By A(i)we denote the arithmetic mean for the observable Â, without the i-th data value Ai , namely
NBin
X
1
A(i) ≡Ak .
NBin − 1 k=1, k6=i(40)

As the corresponding quantity for the function f (〈Â〉, 〈B̂〉, 〈Ĉ〉 · · · ), we define
f(i) (〈Â〉, 〈B̂〉, 〈Ĉ〉 · · · ) ≡ f (A(i) , B(i) , C(i) · · · ) .

(41)

Following the convention in the literature, we will denote the final Jackknife estimate of themean by f(·) and its standard error by ∆ f . The Jackknife mean is given by
NBin
1 X
f(·) (〈Â〉, 〈B̂〉, 〈Ĉ〉 · · · ) =f(i) (〈Â〉, 〈B̂〉, 〈Ĉ〉 · · · ) ,NBin i=1(42)

and the standard error, including bias correction, is given by(∆ f )2 =

NBin
2

NBin − 1 X
f(i) (〈Â〉, 〈B̂〉, 〈Ĉ〉 · · · ) − f(·) (〈Â〉, 〈B̂〉, 〈Ĉ〉 · · · ) .
NBin i=1(43)

In case of f = 〈Â〉, the results (42) and (43) reduce to the plain sample average and thestandard, bias corrected, estimate of the error.
2.4.2

An explicit example of error estimation

In the following we use one of our examples, the Hubbard model on a square lattice in the Mz
Hubbard-Stratonovich decoupling (see Sec. 4.2), to show explicitly how to estimate errors. We
will equally show that the autocorrelation time is dependent upon the choice of the observable.
In fact, different observables within the same run can have different autocorrelation times and
15

SciPost Phys. 3, 013 (2017)

1

b)

SŜ z

0.025

0.8

(SŜ x + SŜ y )/2

0.02

0.6

(SŜ x + SŜ y + SŜ z )/3

0.015

Error

Auto orrelation

a)

0.4

0.01

SŜ z

0.2

0.005

(SŜ x + SŜ y )/2

0

0

0

50

100
tQM C

150

200

(SŜ x + SŜ y + SŜ z )/3
0

100

200
300
NRebin400

500

Figure 1: The autocorrelation function Aut oÔ (t QMC ) (a) and the scaling of the error witheffective bin size (b) of three equal time spin-spin correlation functions Ô of the Hubbard
model in the Mz decoupling (see Sec. 4.2). Simulations were done on a 6 × 6 square lattice,with U/t = 4 and β t = 6. The original bin contained only one sweep and we calculatedaround one million bins on a single core. The different autocorrelation times for the x y-planecompared to the z-direction can be detected from the decay rate of the autocorrelation function(a) and from the point where saturation of the error sets in (b), which defines the required
effective bin size for independent measurements. Apparently and as argued in the text, theimproved estimator (SŜ x + SŜ y + SŜ z )/3 has the smallest autocorrelation time.
of course, this time scale depends on the parameter choice. Hence, the user has to checkautocorrelations of individual observables for each simulation! Typical regions of the phase
diagram that require special attention are critical points where length scales diverge.
To determine the autocorrelation time, we calculate the correlation function
NBin −t QMC
XOi − ÔOi+t QMC − Ô ,Aut oÔ (t QMC ) =(44)
Oi − Ô Oi − Ô
i=0
where Oi refers to the Monte Carlo estimate of the observable Ô in the i th bin. This function typically shows an exponential decay and the decay rate defines the autocorrelation time.
Figure 1 (a) shows the autocorrelation functions Aut oÔ (t QMC ) for three spin-spin-correlation
functions [Eq. (65)] at momentum ~k = (π, π) and at τ = 0: Ô = SŜ z for the z spin direction,Ô = (SŜ x + SŜ y )/2 for the x y plane, and Ô = (SŜ x + SŜ y + SŜ z )/3 for the total spin. The Hubbard model has a SU(2) spin symmetry. However, we chose a HS field which couples to thez-component of the magnetization, Mz , such that each configuration breaks this symmetry. Of
course, after Monte Carlo averaging one expects restoration of the symmetry. The model, on
bipartite lattices, shows spontaneous spin-symmetry breaking at T = 0 and in the thermodynamic limit. At finite temperatures, and within the so-called renormalized classical regime,
quantum antiferromagnets have a length scale that diverges exponentially with decreasing
temperatures [86]. The parameter set chosen for Fig. 1 is non-trivial in the sense that it places
the Hubbard model in this renormalized classical regime where the correlation length is substantial. Figure 1 clearly shows a very short autocorrelation time for the x y-plane whereas
we detect a considerably longer autocorrelation time for the z-direction. This is a direct consequence of the long magnetic length scale and the chosen decoupling. The physical reason
for the long autocorrelation time corresponds to the restoration of the SU(2) spin symmetry.
This insight can be used to define an improved, SU(2) symmetric estimator for the spin-spin
16

SciPost Phys. 3, 013 (2017)

correlation function, namely (SŜ x + SŜ y + SŜ z )/3. Thereby, global spin rotations are no longer
an issue and this improved estimator shows the shortest autocorrelation time as seen clearlyin Fig. 1 (b). Other ways to tackle large autocorrelation can be global updates or paralleltempering.
Using the time series of Monte Carlo samples we would like to obtain estimates of themean and the standard error of the mean. A simple method which we will describe in thistutorial is the rebinning method, also known in the literature as rebatching, where a fixed
number (denoted by N_rebin) of adjacent original bins are aggregated to form a new effective bin. In addition to measuring the decay rate of the autocorrelation function (44), ameasure for the autocorrelation time can be also obtained by the rebinning method. For acomparison to other methods of estimating the autocorrelation time we refer the reader to theliterature [82, 83, 87]. A reliable error analysis requires independent bins, otherwise the error
is typically underestimated. This behavior is observed in Fig. 1 (b), where the effective bin
size has been systematically increased by rebinning. If the effective bin size is smaller than the
autocorrelation time the error will be underestimated. When the effective bin size becomeslarger than the autocorrelation time converging behavior sets in and in this region the error
estimate will be correct.
For the analysis of the Monte Carlo data (see Sec. 3.5), the user can provide a finite value
for N_auto to trigger the computation of autocorrelation functions Aut oÔ (t QMC ) in the ranget QMC = [0, N_auto]. Since these computations are quite time consuming and require many
Monte Carlo bins the default value is N_auto=0 if unspecified. To produce Fig. 1, we set
N_auto = 500 and used a total of approximately one million bins.

17

SciPost Phys. 3, 013 (2017)

2.5

Pseudo code description

Algorithm 1 Basic structure of the auxiliary field QMC implementation in Prog/main.f90
1: call ham_set. Set the Hamiltonian and the lattice
2: call confin. Read in an auxiliary-field configuration or generate it randomly
3: for n = L Trotter to 1 do
4:call wrapul. Fill the storage, needed for the first actual Monte Carlo sweep. Compute propagation matrices and store them at stabilization points

5: end for
6: for nbc = 1 to Nbin do. Loop over bins. The bin defines the unit of Monte Carlo time7:

for nsw = 1 to Nsweep do. Loop over sweeps. Each sweep updates twice. (upward and downward in imag. time) the space-time lattice of auxiliary fields

8:

for nτ = 1 to LTrotter do. Upward sweep
call wrapgrup. Propagate Green fct. from n t au − 1 to nτ , and compute new. estimate of Green fct. at nτ , using sequential updates

9:

10:
11:

12:
13:
14:
15:
16:
17:
18:
19:

20:

. Stabilization:
if nτ = stabilization point in imaginary time thencall wrapur . Compute propagation from previous stabilization point to nτ. Storage management:. Read from storage: propagation from LTrotter to nτ. Write to storage : the just computed propagation
call cgr. Recalculate the Green function at time nτ in a stable way
call control_precisionG. Compare propagated and recalculated Green fct.
end ifif nτ ∈ [LOBS_S T, LOBS_EN ] thencall obser
end ifend for. Measure the equal time observables

for nτ = LTrotter to 1 do
. Downward sweep. Repeat the above steps (update, propagation, stabilization, equal time
. measurements) for the downward direction in imaginary timeend forif nτ = 1 then. Measure the time displaced observables
call tau_m23:
end if24:
end for25:
call pr_obs. Calculate measurement averages for current bin and write them to disk
26:
call confout. Write auxiliary-field configuration to disk
27: end for
21:

22:

18

SciPost Phys. 3, 013 (2017)3

Data Structures and Input/Output

3.1

Implementation of the Hamiltonian and the lattice

The module Hamiltonian, contained in the file Hamiltonian.f90, defines the model
Hamiltonian, the lattice under consideration and the desired observables (Table 2). We
have collected a number of example Hamiltonians, lattices and observables in the file
Hamiltonian_Examples.f90. The examples are described in Sec. 4. To implement a userdefined model, only the module Hamiltonian has to be set up. Accordingly, this documentation focusses almost entirely on this module and the subprograms it includes. The remaining
parts of the code may hence be treated as a black box.
To specify the Hamiltonian, one needs an Operator and a Lattice type as well as a type
for the observables. These three data structures will be described in the following sections.
Subprogram

Description

Ham_SetReads in model and lattice parameters from the file parameters
and it sets the Hamiltonian by calling Ham_latt, Ham_hop,and Ham_V.
Sets the hopping term HˆT by calling Op_make and Op_set.
Sets the interaction terms HˆV and HˆI by calling Op_make
and Op_set.
Sets the lattice by calling Make_Lattice.
A function which returns an update ratio for the Ising term
HˆI,0 .
Assigns memory storage to the observables
Computes the scalar observables and equal-time correlation
functions.
Computes time-displaced correlation functions.
Initializes the observables to zero.
Writes the observables to the disk by calling Print_bin.

Ham_hop
Ham_V
Ham_Latt
S0
Alloc_obs
Obser
ObserT
Init_obs
Pr_obsSection

3.1.1, 3.1.2
3.1.1, 3.1.2
3.1.3
4.4.3

3.2
3.2

Table 2: Overview of the subprograms of the module Hamiltonian to define the Hamiltonian,
the lattice and the observables. The highlighted subroutines have to be modified by the user.

3.1.1

The Operator typeThe fundamental data structure in the code is the data structure Operator. It is implemented
as a Fortran derived data type. This type is used to define the Hamiltonian (2). In general,the matrices T(ks) , V(ks) and I(ks) are sparse Hermitian matrices. Consider the matrix X ofdimension Ndim × Ndim , as a representative for each of the above three matrices. Let us denotewith {z1 , · · · , zN } a subset of N indices, for which
6= 0 if x, y ∈ {z1 , · · · zN }X x, y(45)=0otherwise
Usually, we have N  Ndim . We define the N × Ndim matrices P as
Pi,x = δzi ,x ,(46)

where i ∈ [1, · · · , N ] and x ∈ [1, · · · , Ndim ]. The matrix P selects the non-vanishing entries of
X , which are contained in the rank-N matrix O:
X = P T OP ,
19

(47)

SciPost Phys. 3, 013 (2017)

and
X x, y =N
X

Pi,x Oi, j P j, y =

N
X

i, j

δzi ,x Oi j δz j , y .

(48)

i, j

Since the P matrices have only one non-vanishing entry per column, they can conveniently be~ , with entries
stored as a vector P
Pi = zi .
(49)
There are many useful identities which emerge from this structure. For example:
n∞XP T OPXP T OPe =e== 1 + P T eO − 1 P ,
n!
n=0

(50)

since
P P T = 1N ×N .

(51)

In the code, we define a structure called Operator to capture the above. This typeOperator bundles several components that are needed to define and use an operator matrix in the program.
3.1.2

Specification of the model

VariableType

Description

Op_X%N
Op_X%O
Op_X%P
Op_X%g
Op_X%alpha
Op_X%type

IntegerComplexIntegerComplex
ComplexInteger

Op_X%U
Op_X%E
Op_X%N_non_zero

ComplexReal
Integer

Effective dimension N
Matrix O of dimension N × N
Matrix P encoded as a vector of dimension N
Coupling strength gConstant αParameter to set the type of HS transformation(1 = Ising, 2 = discrete HS for perfect-square term)
Matrix containing the eigenvectors of O
Eigenvalues of O
Number of non-vanishing eigenvalues of O

Table 3: Member variables of the Operator type. In the left column, the letter X is a placeholder for the letters T and V, indicating hopping and interaction operators, respectively. The
highlighted variables have to be specified by the user.
In this section we show how to specify the Hamiltonian (2) in the code. More pre(ks)cisely,to set the matrix representation of the imaginary-time propagators – e−∆τT ,p we have(ks)(ks)
e −∆τUk ηkτ V , and e−∆τskτ I – that appear in the partition function (15). For each pair ofindices (k, s), these terms have the general form
Matrix Exponential = e g φ(type) X .

(52)

In case of the perfect-square term, we additionally have to set the constant α, see the definitionof the operators V̂ (k) in Eq. (4). The data structures which hold all the above information arevariables of the type Operator (see Table 3). For each pair of indices (k, s), we store thefollowing parameters in an Operator variable:~ and O defining the matrix X [see Eq. (47)]• P

20

SciPost Phys. 3, 013 (2017)• the constants g, α
• optionally: the type type of the discrete fields φ
In case of the Ising term, we store type=1 which sets φkτ = skτ . In case of the perfect-squareterm, the field results from the discrete HS transformation (10) and we store type=2 whichsets φkτ = ηkτ . Note that we have dropped the color index σ, since the implementation uses
the SU(Ncol ) invariance of the Hamiltonian.
Accordingly, the following data structures fully describe the Hamiltonian (2):• For the hopping Hamiltonian (3), we have to set the exponentiated hopping matrices(ks)e−∆τT :In this case X(ks) = T(ks) . Precisely, a single variable Op_T describes the operator matrix
N
dim
X

ĉ †x Tx(ks)y ĉ y


,

(53)

x, y

where k = [1, M T ] and s = [1, Nfl ]. To make contact with the general expression (52) we
set g = −∆τ (and α = 0). In case of the hopping matrix, the type variable Op_T%type
is neglected by the code. All in all, the corresponding array of structure variables is
Op_T(M T ,N f l ).
• For the interaction Hamiltonian(4), which is of perfect-square type, we have to set thep−∆τUk ηkτ V (ks)exponentiated matrices e:
In this case, X = V(ks) . A single variable Op_V describes the operator matrix:
N
dim
X

ĉ †x Vx,(ks)y ĉ y




+ αks

,

(54)

x, y

where k = [1, MV ] and s = [1, Nfl ]. To makep contact with the general expression (52)and to set the constant α, we choose g = −∆τUk and α = αks . The discrete HubbardStratonovich decomposition which is used for the perfect-square interaction, is selectedby setting the type variable to Op_V%type = 2. All in all, the required structure variables
Op_V are defined using the array Op_V(MV ,N f l ).
• For the Ising interaction Hamiltonian (5), we have to set the exponentiated matrices(ks)e−∆τskτ I :In this case, X = I (k,s) . A single variable Op_V then describes the operator matrix:
N
dim
X

ĉ †x I x(ks)y ĉ y


,

(55)

x, y

where k = [1, M I ] and s = [1, Nfl ]. To make contact with the general expression (52),
we set g = −∆τ (and α = 0). The Ising interaction is specified by setting the typevariable Op_V%type=1. All in all, the required structure variables are contained in thearray Op_V(M I ,N f l ).
• In case of a full interaction [perfect-square term (4) and Ising term (5)], we define thecorresponding doubled array Op_V(MV +M I ,N f l ) and set the variables separately forboth ranges of the array according to the above.

21

SciPost Phys. 3, 013 (2017)

3.1.3

The Lattice type

We have a lattice module which can generate one- and two-dimensional Bravais lattices. Notethat the orbital structure of each unit cell has to be specified by the user in the Hamiltonian
module. The user has to specify unit vectors a~1 and a~2 as well as the size of the lattice. The
size is characterized by two vectors ~L1 and ~L2 and the lattice is placed on a torus (periodic
boundary conditions):ĉ~i+~L1 = ĉ~i+~L2 = ĉ~i .
(56)
The function call

Call Make_Lattice( L1, L2, a1,

a2, Latt )will generate the lattice Latt of type Lattice. Note again that the orbital structure of theunit cell has to be provided by the user. The reciprocal lattice vectors are defined by:
a~i · g~i = 2πδi, j ,(57)

and the Brillouin zone corresponds to the Wigner-Seitz cell of the lattice. With ~k =the k-space quantization follows from:

 
 
~L1 · g~1 ~L1 · g~2 α1n~L2 · g~1 ~L2 · g~2 α2 = 2π mP~i ,i αi g(58)

such that~k = n~b1 + m~b2 ,
~b1 =

with
2π




( ~L2 · g~2 )~g1 − ( ~L2 · g~1 )~g2 ,( ~L1 · g~1 )( ~L2 · g~2 ) − ( ~L1 · g~2 )( ~L2 · g~1 )


2π
~b2 =
( ~L1 · g~1 )~g2 − ( ~L1 · g~2 )~g1 .
( ~L1 · g~1 )( ~L2 · g~2 ) − ( ~L1 · g~2 )( ~L2 · g~1 )

and
(59)

The Lattice module equally handles the Fourier transformation. For example the subroutine Fourier_R_to_K carries out the transformation:
X ~
1
~ ~
S(~k, :, :, :) =e−i k·(i− j ) S(~i − ~j, :, :, :) ,(60)
Nunit cell ~ ~i, jand Fourier_K_to_R the inverse Fourier transformX ~1
S(~r , :, :, :) =e i k·~r S(~k, :, :, :).
Nunit cell ~

(61)

k∈BZ

In the above, the unspecified dimensions of the structure factor can refer to imaginary-timeand orbital indices.

3.2

The observable types Obser_Vec and Obser_LattOur definition of the model includes observables [Eq. (24)]. We have defined two observabletypes: Obser_vec for an array of scalar observables such as the energy, and Obser_Lattfor correlation functions that have the lattice symmetry. In the latter case, translation symmetry can be used to provide improved estimators and to reduce the size of the output.
We also obtain improved estimators by taking measurements in the imaginary-time interval[LOBS_ST,LOBS_EN] (see the parameter file in Sec. 3.3.1) thereby exploiting the invariance22

SciPost Phys. 3, 013 (2017)

VariableType

Description

Latt%a1_p, Latt%a2_p
Latt%L1_p, Latt%L2_p

Real
Real

Latt%N
Latt%list

Integer
Integer

Latt%invlist
Latt%nnlist

Integer
Integer

Latt%imj
Latt%BZ1_p, Latt%BZ2_p
Latt%b1_p, Latt%b1_p
Latt%listk

Integer
Real
Real
Integer

Latt%invlistk
Latt%b1_perp_p,
Latt%b2_perp_p

Integer

Unit vectors a~1 , a~2
Vectors ~L1 , ~L2 that define the topology of the lattice.
Tilted lattices are thereby possible to implement.
Number of lattice points, Nunit cell
Maps each lattice point i = 1, · · · , Nunit cell to a real
space vector denoting the position of the unit cell:~ i = list(i,1) a~1 + list(i,2) a~2 ≡ i1 a~1 + i2 a~2R
Invlist(i1 , i2 ) = i
j = nnlist(i, n1 , n2 ), n1 , n2 ∈ [−1, 1]~j = R~ i + n1 a~1 + n2 a~2R~~i − R~ j . im j, i, j ∈ 1, · · · , Nunit cell
R im j(i, j) = R
Reciprocal space vectors g~i (See Eq. 57)
k-quantization (See Eq. 59)
Maps each reciprocal lattice point k = 1, · · · , Nunit cellto a reciprocal space vector~kk = listk(k,1)~b1 + listk(k,2)~b2 ≡ k1 ~b1 + k2 ~b2
Invlistk(k1 , k2 ) = k

RealOrthonormal vectors to ~bi . For internal use.

Table 4: Components of the Lattice type for two-dimensional lattices using as example thedefault lattice name Latt. The highlighted variables have to be specified by the user. Other
components of the Lattice are generated upon calling: Call Make_Lattice( L1, L2,
a1, a2, Latt ).
under translation in imaginary-time. Note that the translation symmetries in space and in time
are broken for a given configuration C but restored by the Monte Carlo sampling. In general,the user defines the size and the number of bins in the parameter file, each bin having a given
amount of sweeps. Within a sweep we run sequentially through the HS and Ising fields, fromtime slice 1 to time slice LTrotter and back. The results of each bin are written to a file andanalyzed at the end of the run.
To accomplish the reweighting of observables (see Sec. 2.1.3), for each configuration themeasured value of an observable is multiplied by the factors ZS and ZP:

ZS = sign(C) ,
ZP =

(62)

−S(C)

e


ℜ e−S(C).

(63)

They are computed from the Monte Carlo phase of a configuration,phase =e−S(C)e−S(C),(64)

which is provided by the main program. Note that each observable structure also includes theaverage sign [Eq. (25)].
3.2.1

Scalar observables

This data type is described in Table 5 and is useful to compute an array of scalar observables. Consider a variable Obs of type Obser_vec. At the beginning of each bin, a call to
23

SciPost Phys. 3, 013 (2017)

Obser_Vec_Init in the module observables_mod.f90 will set Obs%N=0, Obs%Ave_sign=0 and Obs%Obs_vec(:)=0. Each time the main program calls the routine Obser in the
Hamiltonian module, the counter Obs%N is incremented by one, the sign (see Eq. 23) isaccumulated in the variable Obs%Ave_sign, and the desired observables (multiplied by the−S(C)sign and ℜ ee−S(C) , see Sec. 2.1.2) are accumulated in the vector Obs%Obs_vec. At the end of][

VariableType

Description

Obs%N
Obs%Ave_sign
Obs%Obs_vec(:)

Int.
Real
Compl.

Number of measurements
Cumulated sign [Eq. (25)]
Cumulated vector ofsign(C)observables [Eq. (24)]〈〈Ô(:)〉〉C ℜ ee−S(C) sign (C)][Obs%File_Vec

Char.

Name of output file

Contribution ofconfiguration C−S(C)Table 5: Components of the Obser_vec type. The table lists the data included in a variable

Obs of type Obser_vec.
the bin, a call to Print_bin_Vec in module observables_mod.f90 will append the resultof the bin in the file File_Vec_scal. Note that this subroutine will automatically append thesuffix _scal to the the filename File_Vec. This suffix is important to allow automatic analysisof the data at the end of the run.
3.2.2

Equal time and time displaced correlation functions

VariableType

Description

Obs%N
Obs%Ave_sign
Obs%Obs_latt

Int.
Real
Compl.

Number of measurements
Cumulated sign [Eq. (25)]
Cumululated correlation

sign(C)

function [Eq. (24)]〈〈Ô~i,α (τ)Ô~j,β 〉〉C ℜ ee−S(C) sign(C)[](~i − ~j, τ, α, β)

Obs%Obs_latt0(α)

Compl.

Contribution ofconfiguration C

Cumulated expectation
value [Eq. (24)]

Obs%File_LattChar.

Name of output file

−S(C)

−S(C)

〈〈Ô~i,α 〉〉C ℜ ee−S(C) sign (C)[]Table 6: Components of the Obser_latt type. The table lists the data included in a variable
Obs of type Obser_latt.
This data type (see Table 6) is useful so as to deal with equal time as well as imaginary-time
displaced correlation functions of the form:SÔ,α,β (~k, τ) =1

X

Nunit cell ~ ~



~ ~ ~
e−k·(i− j ) 〈Ô~i,α (τ)Ô~j,β 〉 − 〈Ô~i,α 〉〈Ô~j,β 〉 .

(65)

i, j

Here, translation symmetry of the Bravais lattice is explicitly taken into account. The correla-

24

SciPost Phys. 3, 013 (2017)tion function splits in a correlated part S
S(back) ~(corr) ~(k):(k, τ) and a background part SÔ,α,βÔ,α,β1~ ~ ~(corr) ~(k, τ) =e−i k·(i− j ) 〈Ô~i,α (τ)Ô~j,β 〉 ,Ô,α,βNunit cell ~ ~
i, j
X

S

(66)

1
~ ~ ~
(back) ~(k) =e−i k·(i− j ) 〈Ô~i,α (τ)〉〈Ô~j,β 〉Ô,α,βNunit cell ~ ~i, jX= Nunit cell 〈Ôα 〉〈Ôβ 〉 δ(~k) ,(67)

where translation invariance in space and time has been exploited to obtain the last line. The
background part depends only on the expectation value 〈Ôα 〉, for which we use the following
estimatorX
1〈Ôα 〉 ≡(68)〈Ô~i,α 〉 .
Nunit cell ~iConsider a variable Obs of type Obser_latt. At the beginning of each bin a call to thesubroutine Obser_Latt_Init in the module observables_mod.f90 will initialize the elements of Obs to zero. Each time the main program calls the Obser or ObserT routines one−S(C)accumulates the quantity 〈〈Ô~i,α (τ)Ô~j,β 〉〉C ℜ ee−S(C) sign(C) in Obs%Obs_latt(~i − ~j, τ, α, β )
[
]
−S(C)
and 〈〈Ô~i,α 〉〉C ℜ ee−S(C) sign (C) in Obs%Obs_latt0(α). At the end of each bin, a call to[]
Print_bin_Latt in the module observables_mod.f90 will append the result of the binin the specified file Obs%File_Latt. Note that the routine Print_bin_Latt carries out theFourier transformation and prints the results in k-space. We have adopted the following naming conventions. For equal time observables, defined by having the second dimension of thearray Obs%Obs_latt(~i − ~j, τ, α, β ) set to unity, the routine Print_bin_Latt attaches thesuffix _eq to Obs%File_Latt. For time displaced correlation functions we use the suffix _tau.

3.3

File structure
Directory

Description

Prog/
Libraries/
Analysis/
Examples/Start/
Documentation/

Main program and subroutines
Collection of mathematical routines
Routines for error analysis
Example simulations for Hubbard-type models
Parameter files and scripts
Documentation of the QMC code.

Table 7: Overview of the directories.
The code package consists of the program directories Prog/, Libraries/ andAnalysis/. The example simulations corresponding to the walkthroughs of Sec. 4.1 - 4.4
are included in Examples/. The package content is summarized in Table 7.
3.3.1

Input files

The input files are listed in Table 8. To enable restarting a previous simulation (see Table 10)
or to use a given HS and Ising field configuration as input for a new simulation, the program
reads in the files confin_<threadnumber> in case they are present. It goes without sayingthat the dimensions of the thereby defined field configuration (number of threads, lattices size,
and number of time slices) have to match the corresponding values of the parameter file. The
25

SciPost Phys. 3, 013 (2017)

File

Description

parameters

This collects input data. We can set here the parameters forthe lattice, which model, variables of the QMC process, andthe error analysis.
List of integer numbers to initialize the random number
generator and to start a simulation from scratch.
(Optionally, a HS and Ising field configuration can beprovided as input.)

seeds(confin_<threadnumber> )Table 8: Overview of the input files in Start/ required for a simulation.
parameter file Start/parameters has the following form – using as an example the SU(2)symmetric Hubbard model on a square lattice (see Sec. 4.1 for a detailed walkthrough):!============================================================================! Variables for the Hubb program!---------------------------------------------------------------------------&VAR_lattice
L1 = 4! Length in direction a_1
L2 = 4! Length in direction a_2
Lattice_type = "Square"! a_1 = (1,0),a_2=(0,1), Norb=1, N_coord=2!Lattice_type ="Honeycomb"! a_1 = (1,0),a_2 =(1/2,sqrt(3)/2),Norb=2,N_coord=3
Model = "Hubbard_SU2"! Sets Nf=1, N_sun=2. HS field couples to the! density!Model = "Hubbard_Mz"! Sets Nf=2, N_sun=1. HS field couples to the! z-component of magnetization.
!Model="Hubbard_SU2_Ising"! Sets Nf_1, N_sun=2 and runs only for the square! lattice! Hubbard model coupled to transverse Ising field/&VAR_Hubbard! Variables for the Hubbard model
ham_T= 1.D0! Hopping parameter
ham_chem= 0.D0
! chemical potential
ham_U= 4.D0! Hubbard interaction
Beta= 5.D0! inverse temperaturedtau= 0.1D0! Thereby Ltrot=Beta/dtau/&VAR_Ising
Ham_xi = 1.d0
Ham_J = 0.2d0
Ham_h = 2.d0/

! Model parameters for the Ising code
! Only needed if Model="Hubbard_SU2_Ising"&VAR_QMC
Nwrap
= 10

! Variables for the QMC run! Stabilization. Green functions is computed from! scratch after each time interval Nwrap*Dtau! Number of sweeps! Number of bins
! 1 for calculation of time displ. Green functions;! 0 otherwise! Start measurements at time slice LOBS_ST! End
measurements at time slice LOBS_EN! Code will stop after CPU_MAX hours.
! If not specified, code will stop after Nbin bins.

NSweep
NBin
Ltau

= 500
= 2
= 1

LOBS_ST = 1
LOBS_EN = 50
CPU_MAX = 0.1/&VAR_errors! Variables for analysis programs

26

SciPost Phys. 3, 013 (2017)n_skip = 1
N_rebin = 1
N_Cov
= 0
N_auto= 100! Number of bins that will be skipped.
! Rebinning! If set to 1 covariance will be computed! for unequal time correlation functions.
! If set to >0 autocorrelation function will be! computed for scalar and equal time observables.

/

3.3.2

Output: Observables

File

Description

info

After completion of the simulation, this file documents parameters of the model, the QMC run and simulation metrics(precision, acceptance rate, wallclock time).
Results of equal time measurements of scalar observables.
The placeholder X stands for the observables Kin,Pot,Part,and Ener.
Results of equal time and time displaced measurements of correlation functions. The placeholder Y stands for Green,SpinZ,
SpinXY, and Den.
Output files for the HS and Ising field configuration.

X_scalY_eq,Y_tau

confout_<threadnumber>Table 9: Overview of the standard output files. See Sec. 3.2 for the definitions of observables
and correlation functions.
The standard output files are listed in Table 9. The output of the measured data is organizedin bins. One bin corresponds to the arithmetic average over a fixed number of individual
measurements which depends on the chosen measurement interval [LOBS_ST,LOBS_EN] on
the imaginary-time axis and on the number NSweep of Monte Carlo sweeps. If the user runs an
MPI parallelized version of the code, the average also extends over the number of MPI threads.
The formatting of the output for a single bin depends on the observable type, Obs_vec or
Obs_Latt:• Observables of type Obs_vec: For each additional bin, a single new line is added to the
output file. In case of an observable with N_size components, the formatting is

N_size+1 <measured value,1> ... <measured value,N_size><measured sign>
The counter variable N_size+1 refers to the number of measurements per line, including the phase measurement. This format is required by the error analysis routine (seeSec. 3.5). Scalar observables like kinetic energy, potential energy, total energy and particle number are treated as a vector of size N_size=1.
• Observables of type Obs_Latt: For each additional bin, a new data block is added to the
output file. The block consists of the expectation values [Eq. (68)] contributing to thebackground part [Eq. (67)] of the correlation function, and the correlated part [Eq. (66)]of the correlation function. For imaginary-time displaced correlation functions, the formatting of the block follows this scheme:<measured sign> <N_orbital> <N_unit_cell> <N_time_slices> <dtau>
do alpha = 1, N_orbital
27

SciPost Phys. 3, 013 (2017)〈Ôα 〉

enddo
do i = 1, N_unit_cell<reciprocal lattice vector k(i)>do tau = 1, N_time_slices
do alpha = 1, N_orbital
do beta = 1, N_orbital〈S(cor r)(k(i), τ)〉Ô,α,β

enddo
enddo
enddo
enddo
The same block structure is used for equal time correlation functions, except for theentries <N_time_slices> and <dtau> which are not present in the latter. Using thisstructure for the bins as input, the full correlation function SÔ,α,β (~k, τ) [Eq. (65)] is thencalculated by calling the error analysis routine (see Sec. 3.5).
3.3.3

Output: Precision

The finite temperature auxiliary field QMC algorithm is known to be numerically unstable,
as discussed in Sec. 2.3. The origin the numerical instabilities arises from the imaginarytime propagation which invariably leads to exponentially small and exponentially large scales.
Numerical stabilization of the code is delicate and has been pioneered in Ref. [2] for the finitetemperature algorithm and in Refs. [3, 4] for the zero temperature projective algorithm. As
shown in Ref. [7] scales can be omitted in the ground state algorithm – thus rendering it very
stable – but have to be taken into account in the finite-temperature code. Apart from runtime
information, the file info contains important information concerning the stability of the code.
It is important to know that numerical stabilization is delicate and there is no guarantee that
it will work for all models.
If the numerical stabilization turns out to be bad, one option is to reduce the value of theparameter Nwrap in the parameter file. For performing the stabilization of the involved matrix
multiplications we rely on routines from LAPACK. Hence it is very likely that your results maychange significantly if you switch the LAPACK implementation. In order to offer a simple
baseline to which people can quickly switch if they want to see whether their results dependon the library used for linear algebra routines we have included parts of the LAPACK-3.6.1
reference implementation from http://www.netlib.org/lapack/. You can switch to the QR
decomposition related routines from the LAPACK reference implementation by including theswitch -DQRREF into the PROGRAMCONFIGURATION string. To use these routines you needto link against a lapack library that implements at least the LAPACK-3.4.0 interface.5
To provide further flexibility, we have kept the history of different stabilization schemes.
Our default strategy is quick and generically works well but we have encountered some models
where it fails. If this applies to your model, you can use the switch -DSTAB2 (stabilization
scheme based on the QR decomposition, but not using the LAPACK reference implementation)or -DSTAB1 (stabilization scheme based on singular value decomposition) in the header ofthe file Makefile and recompile the code.
Typical values for the numerical precision can be found in the examples of Sec. 4 (seeSec. 4.1.3 and 4.2.3).
5

We have encountered some compiling issues with this flag. In particular the older intel ifort compiler version
10.1 fails for all optimization levels.

28

SciPost Phys. 3, 013 (2017)

3.4

Scripts

ScriptDescription

Sectionsetenv.shStart/analysis.shStart/out_to_in.sh

Exports the path variable.
Starts the error analysis.
Copies the output configurations of HS and Ising spinsto the respective input files.

3.6
3.5, 3.6
3.6

Table 10: Overview of the bash script files.

3.5

Analysis programs

ProgramDescription

cov_scal.f90

In combination with the script analysis.sh, the bin files with suffix _scal
are read in, and the corresponding files with suffix _scalJ are produced.
They contain the result of the Jackknife rebinning analysis (see Sec. 2.4).
In combination with the script analysis.sh, the bin files with suffix _eq
are read in, and the corresponding files with suffix _eqJR and _eqJK are
produced. They correspond to correlation functions in real and Fourier
space, respectively.
In combination with the script analysis.sh, the bin files X_tau are read in,and the directories X_kx_ky are produced for all kx and ky greater or equal
to zero. Here X is a place holder from Green, SpinXY, etc as specified in
Alloc_obs(Ltau) (See section 4.1.2). Each directory contains a fileg_kx_ky containing the time displaced correlation function traced over theorbitals. It also contains the covariance matrix if N_cov is set to unity inthe parameter file (see Sec. 3.3.1).
Equally, a directory X_R0 for the local time displaced correlation function
is generated.

cov_eq.f90

cov_tau.f90

Table 11: Overview of analysis programs that are called within the script analysis.sh.
Here we briefly discuss the analysis programs which read in bin/s and carry out the error
analysis. (See Sec. 2.4 for a more detailed discussion.) Error analysis is based on the central
limit theorem, which requires bins to be statistically independent, and also the existence of awell-defined variance for the observable under consideration. The former will be the case if
bins are longer than the autocorrelation time. The latter has to be checked by the user. In theparameter file listed in Sec. 3.3.1, the user can specify how many initial bins should be omitted(variable n_skip). This number should be at least comparable or lager than the autocorrelation time. The analysis of the autocorrelation time is triggered by specifying a positive value
for N_auto that is turned off be default (N_auto = 0). The rebinning variable N_rebin will
merge N_rebin bins into a single new bin. If the autocorrelation time is smaller than the effective bin size, the error should become independent of the bin size and thereby of the variable
N_rebin. Our analysis is based on the Jackknife resampling [57, 85], which includes proper
treatment of the sign. As listed in Table 11 we provide three analysis programs to account for
the three observable types. The programs can be found in the directory Analysis and are
executed by running the bash shell script analysis.sh. In the following, we describe theformatting of the output files mentioned in Table 13.
29

SciPost Phys. 3, 013 (2017)

File

Description

parameters

Contains also variables for the error analysis:n_skip, N_rebin, N_Cov and N_auto (see Sec. 3.3.1)
Monte Carlo bins (see Table 9)

X_scal, Y_eq, Y_tau

Table 12: Standard input files for the error analysis.
File

Description

X_scalJ

Jackknife mean and error of X, where X stands for Kin, Pot, Part,and Ener.
QMC-time resolved autocorrelation and rebinning analysis of X,
where X stands for Kin, Pot, Part, and Ener and N labels thecomponent if X is a vector.
Jackknife mean and error of Y, where Y stands for Green, SpinZ,
SpinXY, and Den.
The suffixes R and K refer to real and reciprocal space, respectively.
Time-resolved and spatially local Jackknife mean and error of Y,
where Y stands for Green, SpinZ, SpinXY, and Den.
QMC-time resolved autocorrelation and rebinning analysis of Y,
where Y stands for Green, SpinZ, SpinXY, and Den.
Time resolved and ~k-dependent Jackknife mean and error of Y,
where Y stands for Green, SpinZ, SpinXY, and Den.

X_scal_Auto_N

Y_eqJR and Y_eqJK

Y_R0/g_R0
Y_eq_Auto_Tr_kx_ky
Y_kx_ky/g_kx_ky

Table 13: Standard output files of the error analysis.
• For the scalar quantities X, the output files X_scalJ have the following formatting:

Effective number of bins, and bins:<N_bin - n_skip>

OBS :

1

<mean(X)><error(X)>

OBS :

2

<mean(sign)><error(sign)><N_bin>• For the autocorrelation analysis of scalar quantities X, the output files X_scal_Auto_N
have the following formatting:do i = 1, N_autotau(i)/n_rebin
enddo

Auto_X(tau)<error( X )>• For the equal time correlation functions Y, the formatting of the output files Y_eqJR and
Y_eqJK follows this structure:

do i = 1, N_unit_cell<k_x(i)><k_y(i)>do alpha = 1, N_orbital
do beta = 1, N_orbitalalpha beta Re<mean(Y)>error(Y)>enddo
30

Re<error(Y)>Im<mean(Y)>Im<

SciPost Phys. 3, 013 (2017)

enddo
enddo
where Re and Im refer to the real and imaginary part, respectively.
• For the autocorrelation analysis of equal time quantities Y, the output files
Y_eq_Auto_Tr_kx_ky have the following formatting:do i = 1, N_autotau(i)/n_rebin
enddo

Auto_Tr[Y](tau)<error( Tr[Y] )>• The imaginary-time displaced correlation functions Y are written to the output files
Y_R0/g_R0, when measured locally in space, and to the output files Y_kx_ky/g_kx_ky
when they are measured ~k-resolved. Both output files have the following formatting:do i = 0, Ltau
tau(i)<mean( Tr[Y] )>enddo<error( Tr[Y])>where Tr corresponds to the trace over the orbital degrees of freedom.

3.6

Running the code

In this section we describe the steps how to compile and run the code, as well as how to
perform the error analysis of the data.
3.6.1

Compilation

The environment variables and the directives to compile the code are set in the following
makefile Makefile:# -DMPI selects MPI.
# -DSTAB1
Alternative stabilization, using the singular value decomposition.
# -DSTAB2
Alternative stabilization, lapack QR with manual pivoting.
#
Packed form of QR factorization is not used.
# (no flag) Default stabilization, using lapack QR with pivoting.
#
Packed form of QR factorization is used.
# -DQRREF
Enables reference lapack implementation of QR decomposition.
# Recommendation: just use the -DMPI flag if you want to run in parallel or#leave it empty for serial jobs.
#
The default stabilization, no flag, is generically the best.
PROGRAMCONFIGURATION = -DMPIPROGRAMCONFIGURATION =f90 = gfortranexport f90
F90OPTFLAGS = -O3 -Wconversion -fcheck=all
F90OPTFLAGS = -O3
export F90OPTFLAGS
F90USEFULFLAGS = -cpp -std=f2003
F90USEFULFLAGS = -cpp
export F90USEFULFLAGS
FL = -c ${F90OPTFLAGS} ${PROGRAMCONFIGURATION}
export FL
DIR = ${CURDIR}
export DIR

31

SciPost Phys. 3, 013 (2017)

Libs = ${DIR}/Libraries/
export Libs
LIB_BLAS_LAPACK = -llapack -lblas
export LIB_BLAS_LAPACK
all: lib ana programlib:
cd Libraries && $(MAKE)
ana:
cd Analysis && $(MAKE)
program:cd Prog && $(MAKE)

clean: cleanall
cleanall: cleanprog cleanlib cleanana
cleanprog:
cd Prog && $(MAKE) clean
cleanlib:
cd Libraries && $(MAKE) clean
cleanana:
cd Analysis && $(MAKE) clean
help:@echo "The following are some of the valid targets of this Makefile"@echo "all, program, lib, ana, clean, cleanall, cleanprog, cleanlib,
cleanana"

In the above, the GNU Fortan compiler gfortran is set.6 We provide a set of options for compilation of the QMC code. The present options are -DMPI, -DQRREF, -DSTAB1, and -DSTAB2.
They can be included in the string variable PROGRAMCONFIGURATION by the user, as shown
above. The program can be compiled and ran either in single-thread mode (default) or in
multi-threading mode (define -DMPI) using the MPI standard for parallelization. The remaining three compiler options select a particular stabilization scheme for the matrix multiplications (see Sec. 3.3.3). To compile the libraries, the analysis routines and the QMC program atonce, just execute the single command:make
To clean up all directories and remove the object files and executables, execute the command
make clean. As can be seen in the above makefile, there exist also rules to compile/clean up
the library, the analysis routines and the QMC program separately.
3.6.2

Starting a simulation

To start a simulation from scratch, the following files have to be present: parameters andseeds. To run a single-thread simulation, for example by using the parameters of one of the
Hubbard models described in Sec. 4, issue the command./Prog/Examples.out
To restart the code using an existing simulation as a starting point, first run the script
out_to_in.sh to set the input configuration files.
6

A known issue with the alternative Intel Fortran compiler ifort is the handling of automatic, temporary arrays
which ifort allocates on the stack. For large system sizes and/or low temperatures this may lead to a runtime
error. One solution is to demand allocation of arrays above a certain size on the heap instead of the stack. Thisis accomplished by the ifort compiler flag -heap-arrays [n] where [n] is the minimal size (in kilobytes, forexample n=1024) of arrays that are allocated on the heap.

32

SciPost Phys. 3, 013 (2017)

3.6.3

Error analysis

Note that the error analysis script requires the presence of the environment variable DIR which
defines the path to the error analysis programs. So before starting the error analysis, one has
to make this variable available which is done by the script setenv.sh. The command issource ./setenv.sh
To perform an error analysis based on the Jackknife resampling method (Sec. 2.4.1) of theMonte Carlo bins for all observables run the script analysis.sh (see Sec. 3.5). In case thatthe parameter N_auto is set to a finite value the script will also trigger the computation of
autocorrelation functions (Sec. 2.4.2).

4

Examples

4.1

The SU(2)-Hubbard model on a square lattice

To implement a Hamiltonian, the user has to provide a module which specifies the lattice, themodel, as well as the observables they wish to compute. In this section, we describe the module
Hamiltonian_Examples.f90 which contains an implementation of the Hubbard model on
the square lattice. A sample run for this model can be found in Examples/Hubbard_SU2_
Square/. The input files are parameters and seeds (see Tab. 8). The output files are info,
confout, and files with suffixes _scal, _eq, and _tau that contain the raw measurements(see Tab. 9).
The Hamiltonian reads
2
 2
2 NunitX
Xcell

UX X ††
H =c xσ Tx, y c yσ +.
(69)
c c − 1/2
2 x σ=1 xσ xσ
σ=1 x, y=1
We can make contact with the general form of the Hamiltonian by setting: Nfl = 1,
U1
(ks)
Ncol ≡ N_SUN = 2, M T = 1, Tx(ks)y = T x, y , MV = Nunit cell , Uk = − 2 , Vx y = δ x, y δ x,k , αks = − 2
and M I = 0.
4.1.1

Setting the Hamiltonian: Ham_set

The main program will call the subroutine Ham_set in the module Hamiltonian_Hub.f90.
The latter subroutine defines the public variables

Type (Operator), dimension(:,:), allocatable
Type (Operator), dimension(:,:), allocatable
Integer, allocatable :: nsigma(:,:)
Integer:: Ndim, N_FL, N_SUN,:: Op_V:: Op_T
Ltrot

which specify the model. The array nsigma contains the HS field. The routine Ham_set
will first read the parameter file, then set the lattice, Call Ham_latt, set the hopping, Call
Ham_hop, and set the interaction, call Ham_V. The parameters are read in from the fileparameters, see Sec. 3.3.1.
The lattice: Call Ham_latt The choice Lattice_type = "Square" sets a~1 = (1, 0)
and a~2 = (0, 1) and for an L1 × L2 lattice ~L1 = L1 a~1 and ~L2 = L2 a~2 . The call to Call
Make_Lattice( L1, L2, a1, a2, Latt) will generate the lattice Latt of type Lattice.
For the Hubbard model on the square lattice, the number of orbitals per unit cell is given by
NORB=1 such that Ndim ≡ Nunit cell · NORB = Latt%N · NORB, since Nunit cell = Latt%N.
33

SciPost Phys. 3, 013 (2017)

The hopping term: Call Ham_hop The hopping matrix is implemented as follows. We
allocate an array of dimension 1 × 1 of type operator called Op_T and set the dimension forthe hopping matrix to N = Ndim . One allocates and initializes this type by a single call to thesubroutine Op_make:call Op_make(Op_T(1,1),Ndim)
Since the hopping does not break down into small blocks, we have P = 1 andDo i= 1,Ndim
Op_T(1,1)%P(i) = i
Enddo
We set the hopping matrix withDO I = 1, Latt%N
Ix = Latt%nnlist(I,1,0)
Iy = Latt%nnlist(I,0,1)
Op_T(1,1)%O(I ,Ix) = cmplx(-Ham_T,
0.d0,kind(0.D0))
Op_T(1,1)%O(Ix, I ) = cmplx(-Ham_T,
0.d0,kind(0.D0))
Op_T(1,1)%O(I ,Iy) = cmplx(-Ham_T,
0.d0,kind(0.D0))
Op_T(1,1)%O(Iy, I ) = cmplx(-Ham_T,
0.d0,kind(0.D0))
Op_T(1,1)%O(I ,I ) = cmplx(-Ham_chem,0.d0,kind(0.D0))
ENDDO
Here, the integer function j= Latt%nnlist(I,n,m) is defined in the lattice module andreturns the index of the lattice site ~I + n~
a1 + m~a2 . Note that periodic boundary conditions are already taken into account. The hopping parameter Ham_T as well as the chemical potential Ham_chem are read from the parameter file. To completely define the hopping we further set: Op_T(1,1)%g = -Dtau , Op_T(1,1)%alpha = cmplx(0.d0,0.d0,
kind(0.D0)) and call the routine Op_set(Op_T(1,1)) so as to generate the unitary transformation and eigenvalues as specified in Table 3. Recall that for the hopping, the variable
Op_set(Op_T(1,1))%type is not required. Note that although a checkerboard decomposition is not used here, it can be implemented by considering a larger number of sparse hopping
matrices.
The interaction term: Call Ham_V To implement this interaction, we allocate an array of
Operator type. The array is called Op_V and has dimensions Ndim × Nfl = Ndim × 1. We set
the dimension for the interaction term to N = 1, and allocate and initialize this array of typeOperator by repeatedly calling the subroutine Op_make:do i = 1,Ndim
call Op_make(Op_V(i,1),1)
enddo
For each lattice site i, the matrices P are of dimension 1×Ndim and have only one non-vanishingentry. Thereby we can set:
Do i = 1,Ndim
Op_V(i,1)%P(1)

= i

34

SciPost Phys. 3, 013 (2017)

Op_V(i,1)%O(1,1) = cmplx(1.d0,0.d0, kind(0.D0))
Op_V(i,1)%g= sqrt(cmplx(-dtau*ham_U/dble(N_SUN),0.D0,kind(0.D0)))
Op_V(i,1)%alpha = cmplx(-0.5d0,0.d0, kind(0.D0))
Op_V(i,1)%type= 2Call Op_set( Op_V(i,1) )Enddo

so as to completely define the interaction term.
4.1.2

Observables

At this point, all the information for the simulation to start has been provided. The code willsequentially go through the operator list Op_V and update the fields. Between time slices
LOBS_ST and LOBS_EN the main program will call the routine Obser(GR,Phase,Ntau)
which is provided by the user and handles equal time correlation functions. If Ltau=1 themain program will call the routine ObserT(NT, GT0,G0T,G00,GTT, PHASE) which isagain provided by the user and handles imaginary-time displaced correlation functions.
The user will have to implement the observables he/she wants to compute. Here we will
describe how to proceed.
Allocating space for the observables: Call Alloc_obs(Ltau)
observables, the user will have to declare the following:

For four scalar or vector

Allocate ( Obs_scal(4) )Do I = 1,Size(Obs_scal,1)
select case (I)case (1)N = 2; Filename ="Kin"case (2)N = 1; Filename ="Pot"case (3)N = 1; Filename ="Part"
case (4)N = 1, Filename ="Ener"case default
Write(6,*) ’ Error in Alloc_obs ’
end selectCall Obser_Vec_make(Obs_scal(I),N,Filename)
enddo
Here, Obs_scal(1) contains a vector of two observables so as to account for the x- andy-components of the kinetic energy for example.
For equal time correlation functions we allocate Obs_eq of type Obser_Latt. Here we
include the calculation of spin-spin and density-density correlation functions alongside equal
time Green functions.

Allocate ( Obs_eq(4) )Do I = 1,Size(Obs_eq,1)
select case (I)case (1)Ns = Latt%N; No = Norb;

Filename ="Green"
35

SciPost Phys. 3, 013 (2017)case (2)Ns = Latt%N; No = Norb; Filename ="SpinZ"
case (3)Ns = Latt%N; No = Norb; Filename ="SpinXY"
case (4)Ns = Latt%N; No = Norb; Filename ="Den"case default
Write(6,*) ’ Error in Alloc_obs ’
end selectNt = 1
Call Obser_Latt_make(Obs_eq(I),Ns,Nt,No,Filename)enddo
For the Hubbard model Norb = 1 and for equal time correlation functions Nt = 1. If Ltau
= 1 then the code will allocate space for time displaced quantities. The same structure as
for equal time correlation functions will be used albeit with Nt = Ltrot + 1. At the beginning of each bin, the main program will set the bin observables to zero by calling the routine
Init_obs(Ltau). The user does not have to edit this routine.
Measuring equal time observables: Obser(GR,Phase,Ntau)
tion,
GR(x,y, σ) = 〈c x,σ c †y,σ 〉,The equal time Green func(70)the phase factor phase [Eq. (64)], and time slice Ntau are provided by the main program.
Here, x and y label both unit cell as well as the orbital within the unit cell. For the Hubbard
model described here, x corresponds to the unit cell. The Green function does not depend on
the color index, and is diagonal in flavor. For the SU(2) symmetric implementation there is
only one flavor, σ = 1 and the Green function is independent on the spin index. This renders
the calculation of the observables particularly easy. P
An explicit calculation of the potential energy 〈U ~i n̂~i,↑ n̂~i,↓ 〉 reads
Obs_scal(2)%N= Obs_scal(2)%N + 1
Obs_scal(2)%Ave_sign = Obs_scal(2)%Ave_sign + Real(ZS,kind(0.d0))
Do i = 1,Ndim
Obs_scal(2)%Obs_vec(1)=Obs_scal(2)%Obs_vec(1)+(1-GR(i,i,1))**2*Ham_U*ZS*ZP
Enddo
−S(C)

Here ZS = sign(C) [see Eq. (23)], ZP = ℜ ee−S(C) [see Eq. (64)] and Ham_U corresponds to the[]
Hubbard-U term.
Equal time correlations are also computed in this routine. As an explicit example, we consider
the equal time density-density correlation:〈n~i,α n~j,β 〉 − 〈n~i,α 〉〈n~j,β 〉 .

(71)

For the calculation of such quantities, it is convenient to define:GRC(x,y,s) = δ x, y − GR(y,x,s) ,(72)

such that GRC(x,y,s) corresponds to 〈〈ĉ †x,s ĉ y,s 〉〉. In the program code, the calculation of theequal time density-density correlation function looks as follows:

36

SciPost Phys. 3, 013 (2017)

Obs_eq(4)%N= Obs_eq(4)%N + 1! Even if it is redundant, each observable! carries its own counter and sign.
Obs_eq(4)%Ave_sign = Obs_eq(4)%Ave_sign + Real(ZS,kind(0.d0))
Do I1 = 1,Ndim
I= List(I1,1)! = I1 (The Hubbard model on the squareno_I = List(I1,2)! = 1
lattice has one orbital per unit!
cell)
Do J1 = 1,Ndim
J= List(J1,1)
no_J = List(J1,2)
imj = latt%imj(I,J)
Obs_eq(4)%Obs_Latt(imj,1,no_I,no_J) = && Obs_eq(4)%Obs_Latt(imj,1,no_I,no_J) + &&(GRC(I1,I1,1) * GRC(J1,J1,1) * N_SUN * N_SUN+ &&GRC(I1,J1,1) * GR(I1,J1,1) * N_SUN) * ZP * ZS
Enddo
Obs_eq(4)%Obs_Latt0(no_I) = && Obs_eq(4)%Obs_Latt0(no_I)+GRC(I1,I1,1) * N_SUN * ZP * ZS
Enddo

Note that we consider the square lattice of the single site Hubbard model as a special caseof a multiorbital problem as described in Sec. 4.3.1 At the end of each bin the main program
will call the routine Pr_obs(LTAU). This routine will append the result of the bins in thespecified file, with appropriate suffix.
Measuring time displaced observables: ObserT(NT, GT0,G0T,G00,GTT, PHASE) This
subroutine is called by the main program at the beginning of each sweep, provided that LTAU
is set to unity. NT runs from 0 to Ltrot and denotes the imaginary time difference. For a given
time displacement, the main program provides:

GT0(x,y,s) = 〈〈ĉ x,s (N t∆τ)ĉ †y,s (0)〉〉 = 〈〈T ĉ x,s (N t∆τ)ĉ †y,s (0)〉〉G0T(x,y,s) = −〈〈ĉ †y,s (N t∆τ)ĉ x,s (0)〉〉 = 〈〈T ĉ x,s (0)ĉ †y,s (N t∆τ)〉〉G00(x,y,s) = 〈〈ĉ x,s (0)ĉ †y,s (0)〉〉GTT(x,y,s) = 〈〈ĉ x,s (N t∆τ)ĉ †y,s (N t∆τ)〉〉(73)

In the above we have omitted the color index since the Green functions are color independent.
The time displaced spin-spin correlations 4〈〈Ŝ~z (τ)Ŝ~z (0)〉〉 are thereby given by:ij

4〈〈Ŝ~z (τ)Ŝ~z (0)〉〉 = −2 G0T(J1,I1,1) GT0(I1,J1,1) .
i

j

(74)

Note that the above holds for the SU(2) HS transformation discussed in this chapter. The
handling of time displaced correlation functions is identical to that of equal time correlations.
4.1.3

Numerical precision

The directory Examples/Hubbard_SU2_Square contains an example simulation of the 4 × 4
Hubbard model at U/t = 4 and β t = 10. Information on the numerical stability is includedin the following lines of the corresponding file info:
Precision Green Mean, Max : 1.2918865817224671E-014
Precision Phase, Max: 5.0272908791449966E-012
Precision tau
Mean, Max : 8.4596701790588625E-015

37

4.0983018995027644E-011
3.5033530012121281E-011

SciPost Phys. 3, 013 (2017)

showing the mean and maximum difference between the wrapped and from scratched computed equal and time displaced Green functions [7]. A stable code should produce results
where the mean difference is smaller than the stochastic error. The above example shows avery stable simulation since the Green function is of order one.

4.2

The Mz -Hubbard model on a square lattice

The Hubbard Hamiltonian can equally be written as:
H =2 NunitX
Xcell
σ=1 x, y=1c †xσ Tx, y c yσ −

2U X †c x,↑ c x↑ − c †x,↓ c x↓ .
2 x

(75)

We can make contact with the general form of the Hamiltonian (see Eq. 2) by setting: Nfl = 2,U(k,s=1)
Ncol ≡ N_SUN = 1, M T = 1, Tx(ks)= δ x, y δ x,k ,y = T x, y , MV = Nunit cell , Uk = 2 , Vx y

Vx(k,s=2)= −δ x, y δ x,k , αks = 0 and M I = 0. The coupling of the HS fields to the z-component ofythe magnetization breaks the SU(2) spin symmetry. Nevertheless the z-component of the spin
remains a good quantum number such that the imaginary-time propagator – for a given HS
field – is block diagonal in this quantum number. This corresponds to the flavor index whichruns from one to two labelling spin up and spin down degrees of freedom. In the parameter
file listed in Sec. 3.3.1 setting the model variable to Hubbard_Mz will carry out the simulation in the above representation. With respect to the SU(2) case, the changes required in the
Hamiltonian_Examples.f90 module are minimal and essentially effect only the interaction term, and the calculation of observables. We note that in this formulation the hopping
matrix can be flavor dependent such that a Zeeman magnetic field can be introduced. If thechemical potential is set to zero, this will not generate a negative sign problem [37, 88, 89].
A sample run for this model can be found in Examples/Hubbard_Mz_Square/. The input
files are parameters and seeds (see Tab. 8). The output files are info, confout, and files
with suffixes _scal, _eq, and _tau that contain the raw measurements (see Tab. 9).
4.2.1

The interaction term: Call Ham_V

The interaction term is now given by:
Allocate(Op_V(Ndim,N_FL))
do nf = 1,N_FL
do i = 1, Ndim
Call Op_make(Op_V(i,nf),1)
enddo
enddo
Do nf = 1,N_FL
nc = 0
X = 1.d0
if (nf == 2) X = -1.d0
Do i = 1,Ndim
nc = nc + 1
Op_V(nc,nf)%P(1) = I
Op_V(nc,nf)%O(1,1) = cmplx(1.d0, 0.d0, kind(0.D0))
Op_V(nc,nf)%g= X*SQRT(CMPLX(DTAU*ham_U/2.d0, 0.D0, kind(0.D0)))
Op_V(nc,nf)%alpha = cmplx(0.d0, 0.d0, kind(0.D0))
Op_V(nc,nf)%type= 2Call Op_set( Op_V(nc,nf) )
Enddo
Enddo

In the above, one will see explicitly that there is a sign difference between the coupling of theHS field in the two flavor sectors.
38

SciPost Phys. 3, 013 (2017)

4.2.2

The measurements: Call Obser, Call ObserT

Since the spin up and spin down Green functions differ for a given HS configuration, the Wickdecomposition will take a different form. In particular, the equal time spin-spin correlation
functions 4〈〈Ŝ~z Ŝ~z 〉〉 calculated in the subroutine Obser will take the form:ij

4〈〈Ŝ zx Ŝ zy 〉〉 = GRC(x,y,1) * GR(x,y,1) + GRC(x,y,2) * GR(x,y,2) +(GRC(x,x,2) - GRC(x,x,1))*(GRC(y,y,2) - GRC(y,y,1))Here, GRC is defined in Eq. 72. Equivalent changes will have to be carried out for other equal
time and time displaced observables.
Apart from these modifications, the program will run in exactly the same manner as for
the SU(2) case.
4.2.3

Numerical precision

The directory Examples/Hubbard_Mz_Square contains an example simulation of the 4 × 4
Hubbard model at U/t = 4 and β t = 10. Information on the numerical stability is includedin the following lines of the corresponding file info:
Precision Green Mean, Max : 5.0823874429126405E-011
Precision Phase, Max: 0.0000000000000000
Precision tau
Mean, Max : 1.5929357848647394E-011

5.8621144596315844E-006
1.0985132530727526E-005

This is still an excellent precision but nevertheless choosing a HS field which couples tothe z-component of the magnetization apparently leads to numerical results that are a coupleof order of magnitudes less precise than a HS decomposition coupling to the charge (comparewith Sec. 4.1.3).

4.3

The SU(2)-Hubbard model on the honeycomb lattice

The Hamilton module Hamiltonian_Examples.f90 can also carry out simulations for the
Hubbard model on the Honeycomb lattice by setting in the parameter file Lattice_type="Honeycomb" (see Sec. 3.3.1). A sample run for this model can be found in Examples/
Hubbard_SU2_Honeycomb/. The input files are parameters and seeds (see Tab. 8). The
output files are info, confout, and files with suffixes _scal, _eq, and _tau that contain the
raw measurements (see Tab. 9).
4.3.1

Working with multi-orbital unit cells: Call Ham_Latt

This model is an example of a multi-orbital unit cell, and the aim of this section is to document
how to implement this in the code. The Honeycomb lattice is a triangular Bravais lattice with
two orbitals per unit cell. The routine Ham_Latt will set:

Norb
= 2
N_coord = 3
a1_p(1) = 1.D0
; a1_p(2) =a2_p(1) = 0.5D0 ; a2_p(2) =L1_p= dble(L1) * a1_p
L2_p
= dble(L2) * a2_p0.d0
sqrt(3.D0)/2.D0

and then call Make_Lattice(L1_p,L2_p,a1_p,a2_p,Latt) so as to generate the triangular lattice. The coordination number of this lattice is N_coord=3 and the number
39

SciPost Phys. 3, 013 (2017)of orbitals per unit cell corresponds to NORB=2. The total number of orbitals is thereby:
Ndim =Latt%N*NORB. To easily keep track of the orbital and unit cell, we define a super-index
as shown below:
Allocate (List(Ndim,2), Invlist(Latt%N,Norb))
nc = 0Do I = 1,Latt%N! Unit-cell index
Do no = 1,Norb! Orbital index
nc = nc + 1! Super-index labeling unit cell
! and orbital
List(nc,1) = I! Unit-cell of super index nc
List(nc,2) = no
! Orbital of super index nc
Invlist(I,no) = nc
! Super-index for given unit cell
! and orbital
Enddo
Enddo

With the above lists one can run through all the orbitals and at each time keep track of
the unit-cell and orbital index. We note that when translation symmetry is completely absent
one can work with a single unit cell, and the number of orbitals will then correspond to thenumber of lattice sites.
4.3.2

The hopping term: Call Ham_Hop

Some care has to be taken when setting the hopping matrix.
Hamiltonian_Examples.f90 we do this in the following way:
DO I = 1, Latt%N
do no = 1,Norb

In the Hamilton module! Loop over unit cell
! Runs over orbitals and! sets the chemical potential

I1 = invlist(I,no)
Op_T(nc,n)%O(I1 ,I1) = cmplx(-Ham_chem, 0.d0, kind(0.D0))
enddo
I1 = Invlist(I,1)! Orbital A of unit cell I
Do nc1 = 1,N_coord! Loop over coordination number
select case (nc1)case (1)J1 = invlist(I,2)! Orbital B of unit cell icase (2)J1=invlist(Latt%nnlist(I,1,-1),2) ! Orbital B of unit cell i+a_1-a_2
case (3)J1=invlist(Latt%nnlist(I,0,-1),2) ! Orbital B of unit cell i-a_2
case default
Write(6,*) ’ Error in Ham_Hop ’
end selectOp_T(nc,n)%O(I1,J1) = cmplx(-Ham_T,
0.d0, kind(0.D0))
Op_T(nc,n)%O(J1,I1) = cmplx(-Ham_T,
0.d0, kind(0.D0))
Enddo
Enddo

As apparent from the above, hopping matrix elements are non-zero only between the A and B
sublattices.
4.3.3

Observables: Call Obser, Call ObserTIn the multi-orbital case, the correlation functions have additional orbital indices. This isautomatically taken care of in the routines Call Obser and Call ObserT since we have40

SciPost Phys. 3, 013 (2017)

already considered the Hubbard model on the square lattice to correspond to a multi-orbital
unit cell albeit with the special choice of one orbital per unit cell.

4.4

The SU(2)-Hubbard model on a square lattice coupled to a transverse Isingfield

The model we consider here is very similar to the above, but has an additional coupling to atransverse field:
H =

2 X
X
σ=1 x, yUc †xσ Tx, y c yσ +
2

 2
X X
x

σ=1

2
c †xσ c xσ − 1/2



+ξ

X



Ẑ〈x, y〉 c †xσ c yσ + h.c.

σ,〈x, y〉−hX〈x, y〉X̂ 〈x, y〉 − JXẐ〈x, y〉 Ẑ〈x 0 , y 0 〉 (76)〈〈x, y〉〈x 0 , y 0 〉〉We can make contact with the general form of the Hamiltonian by setting: Nfl = 1,U
(ks)
Ncol ≡ N_SUN = 2, M T = 1, Tx(ks)y = T x, y , MV = Nunit cell ≡ Ndim , Uk = − 2 , Vx y = δ x, y δ x,k ,αks = − 12 and M I = 2Nunit cell . The last two terms of the above Hamiltonian describe a transverse Ising field model on the bonds of the square lattice. This type of Hamiltonian has recently
been extensively discussed [19, 22, 90]. Here we adopt the notation of Ref. [19]. Note that〈〈x, y〉〈x 0 , y 0 〉〉 denotes nearest neighbor bonds. The modifications required to generalize the
Hubbard model code to the above model are two-fold. First, one has to specify the function Real(Kind=8)functionS0(n,nt), and second, modify the interaction Call Ham_V.
A sample run for this model can be found in Examples/Hubbard_SU2_Ising_Square/.
The input files are parameters and seeds (see Tab. 8). The output files are info, confout,
and files with suffixes _scal, _eq, and _tau that contain the raw measurements (see Tab. 9).
4.4.1

The Ising term

Since the Ising field lives on bonds we have to provide a data structure defining this quantity.
A bond has an anchor site as well as an orientation. The routine Setup_Ising_action
initializes the arrays L_bond and L_bond_inv that contain this information.

nc = 0Do n_orientation = 1,N_coord
Do I = 1, Latt%N
nc = nc + 1
L_bond(I,n_orientation) = nc
L_bond_inv(nc,1) = I
L_bond_inv(nc,2) = n_orientation
enddo
enddo
The two legs of the bond are given by the anchor ~I and ~I + a~norientation .
4.4.2

The interaction term: Call Ham_V

The dimension of Op_V is now (M I + MV ) × Nfl = ((Ncoor d + 1)Ndim ) × 1 since each site has
Ncoor d = 2 bonds for the square lattice.
do i = 1,N_coord*Ndim
call Op_make(Op_V(i,1),2)! Runs over bonds for Ising inter.

41

SciPost Phys. 3, 013 (2017)

enddo
do i = N_coord*Ndim+1, (N_coord+1)*Ndim ! Runs over sites for Hubbard inter.
call Op_make(Op_V(i,1),1)
enddo

The first N_coord*Ndim operators run through the 2N bonds of the square lattice and are
given by:
Do nc = 1,Ndim*N_coord

I1 = L_bond_inv(nc,1)! Runs over bonds. Coordination number = 2.
! For the square lattice Ndim = Latt%N! Anchor of the bond
! L_bond_inv is setup in Setup_Ising_action== 1 ) I2 = Latt%nnlist(I1,1,0) ! Second site of== 2 ) I2 = Latt%nnlist(I1,0,1) ! the bondif ( L_bond_inv(nc,2)
if ( L_bond_inv(nc,2)
Op_V(nc,1)%P(1) = I1
Op_V(nc,1)%P(2) = I2
Op_V(nc,1)%O(1,2) = cmplx(1.d0 ,0.d0, kind(0.D0))
Op_V(nc,1)%O(2,1) = cmplx(1.d0 ,0.d0, kind(0.D0))
Op_V(nc,1)%g= cmplx(-dtau*Ham_xi,0.D0,kind(0.D0))
Op_V(nc,1)%alpha = cmplx(0d0,0.d0, kind(0.D0))
Op_V(nc,1)%type
= 1
Enddo

Here, ham_xi defines the coupling strength between the Ising and fermion degree of freedom.
As for the Hubbard case, the last Ndim operators read:nc = N_coord*Ndim
Do i = 1, Ndim
nc = nc + 1
Op_V(nc,1)%P(1) = i
Op_V(nc,1)%O(1,1)= cmplx(1.d0 ,0.d0, kind(0.D0))
Op_V(nc,1)%g= sqrt(cmplx(-dtau*ham_U/(DBLE(N_SUN)),0.D0,kind(0.D0)))
Op_V(nc,1)%alpha = cmplx(-0.5d0,0.d0, kind(0.D0))
Op_V(nc,1)%type = 2
Enddo

4.4.3

The function Real (Kind=8) function S0(n,nt)

As mentioned above, a configuration now includes both HS spins and Ising spins and is givenby
C = si,τ , l j,τ with i = 1 · · · M I , j = 1 · · · MV , τ = 1, L Tr ot t er .
(77)
This configuration is stored in the integer array nsigma(M_V + M_I, Ltrot). With theabove ordering of Hubbard and Ising interaction terms, and a for a given imaginary time,
the first 2*Ndim fields correspond to the Ising interaction and the next Ndim ones to the
Hubbard interaction. The first argument of the function S0, namely n, corresponds to theindex of the operator string Op_V(n,1). If Op_V(n,1)%type = 2 then S0(n,nt) returns 1. Note that type=2 refers to spins that stem from a HS transformation. If however
Op_V(n,1)%type = 1 then function S0 returns
e−S0,I s1,τ ,··· ,−sn,τ ,···sMI ,τ




(78)
e−S0,I s1,τ ,··· ,sn,τ ,···sMI ,τ
That is, if n ≤ 2 ∗ Ndim, S0(n,nt) returns the ratio of the new weight to the old weight of theIsing Hamiltonian upon flipping a single Ising spin sn,τ . Otherwise, S0(n,nt) returns unity.

42

SciPost Phys. 3, 013 (2017)

5

Miscellaneous

5.1

Other models

The aim of this section is to briefly mention a small selection of other models that can bestudied using the QMC code of the ALF project.
5.1.1

Kondo lattice model

Simulating the Kondo lattice with the QMC code of the ALF project requires rewriting of themodel along the lines of Refs. [20, 21, 91]. Adopting the notation of these articles, the Hamiltonian that one will simulate reads:
2
 J X XX 
2
U X f†††Hˆ = −tĉ~ ĉ~ + H.c. −+ĉ~ fˆ~ + fˆ~ ĉ~n̂~ − 1 .
(79)
i,σ j,σi,σ i,σi,σ i,σ
i
4
2
σ
~i
~
〈~i,~j〉,σ| i {z
}
|
{z
}
≡HˆU

≡Hˆt

This form is included in the general Hamiltonian (2) such that the above Hamiltonian can be
implemented in our program package. The relation to the Kondo lattice model follows fromexpanding the square of the hybridization to obtain:Xˆ f + η̂z,c · η̂z, f − η̂ x,c · η̂ x, f − η̂ y,c · η̂ y, f + Hˆ .
ˆ c · S~
Hˆ = Hˆt + J
S~
(80)
U
~
~
~
~
~
~
~
~
~i

i

i

i

i

i

i

i

i

where the η-operators relate to the spin-operators via a particle-hole transformation in one
spin sector:η̂~α = P̂ −1 Ŝ~α P̂ with P̂ −1 ĉ~ P̂ = (−1)i x +i y ĉ~† and P̂ −1 ĉ~ P̂ = ĉ~
i

i

i,↓

i,↑

i,↑

n̂

i,↓

(81)

f

Since the η̂ f - and Ŝ f -operators do not alter the parity [(−1) ~i ] of the f -sites,Hˆ , HˆU = 0.

(82)

Thereby, and for positive values of U , doubly occupied or empty f -sites – corresponding toeven parity sites – are suppressed by a Boltzmann factor e−β U/2 in comparison to odd parity
sites. Choosing β U adequately essentially allows to restrict the Hilbert space to odd parity
f -sites. In this Hilbert space η̂ x, f = η̂ y, f = η̂z, f = 0 such that the Hamiltonian (79) reduces
to the Kondo lattice model.
5.1.2 SU(N )-Hubbard-Heisenberg models
SU(2N )-Hubbard-Heisenberg [26, 27] models can be written as:
Hˆ = −tX〈~i,~j〉| U X
X
N 2
†
†
†~
~ĉ † ~ĉ + H.c. − J~
D̂~ ~ D̂~ ~ + D̂~ ~ D̂~ ~ +ĉ~ ĉ~ −~i ~j
i, j i, j
i, j i, j
i i
2N ~ ~N ~2〈i, j〉
i
{z
}
{z
} |
{z
} |
≡Hˆt

(83)

≡HˆU

≡HˆJ

Here, ~ĉ~† = (ĉ~† , ĉ~† , · · · , ĉ~† ) is an N -flavored spinor, and D̂~i,~j = ~ĉ~† ~ĉ~j . To use the QMC code
i

i,1

i,2

i,N

i

of the ALF package to simulate this model, one will rewrite the J-term as a sum of perfectsquares,2 
2J X †
HˆJ = −
D̂ ~ ~ + D̂〈~i,~j〉 − D̂†~ ~ − D̂〈~i,~j〉 ,(84)〈i, j〉〈i, j〉
4N ~ ~〈i, j〉

43

SciPost Phys. 3, 013 (2017)

so to manifestly bring it into the form of the general Hamiltonian(2). It is amusing to notethat setting the hopping t = 0, charge fluctuations will be suppressed by the Boltzmann factor2−β U/N ~ĉ~† ~ĉ~ − N2iiesince in this case ĤJ , HˆU = 0. This provides a route to use the auxiliary
field QMC algorithm to simulate – free of the sign problem – SU(2N )-Heisenberg models inthe self-adjoint antisymmetric representation 7 . For odd values of N recent progress in our
understanding of the origins of the sign problem [40] allows us to simulate a set of non-trivial
Hamiltonians [19, 92], without encountering the sign problem.

5.2

Performance, memory requirements and parallelization

As mentioned in the introduction, the auxiliary field QMC algorithm scales linearly in inverse
temperature β and cubic in the volume Ndim . Using fast updates, a single spin flip requires(Ndim )2 operations to update the Green function upon acceptance. As there are LTrotter × Ndim
spins to be visited, the total computational cost for one sweep is of the order of β(Ndim )3 . This
operation dominates the performance, see Fig. 2. A profiling analysis of our code shows that80-90% of the CPU time is spend in ZGEMM calls of the BLAS library provided in the MKL
package by Intel. Consequently, the single-core performance is next to optimal.

Wall- lo k time

10000

QMC simulation
t ∼ x3

1000

100

10
256

400

576

784 1024 1296

Ndim

Figure 2: Volume scaling behavior of the auxiliary field QMC code of the ALF project on SuperMUC (phase 2/Haswell nodes) at the LRZ in Munich. The number of sites Ndim corresponds
to the system volume. The plot confirms that the leading scaling order is due to matrix multiplications such that the runtime is dominated by calls to ZGEMM.
For the implementation which scales linearly in β, one has to store LTrotter /NWrap intermediate propagation matrices of dimension N × N . For large lattices and/or low temperatures
this dominates the total memory requirements that can exceed 2 GB memory for a sequential
version.
At the heart of Monte Carlo schemes lies a random walk through the given configuration
space. This is easily parallalized via MPI by associating one random walker to each MPI task.
For each task, we start from a random configuration and have to invest the autocorrelation
time Tauto to produce an equilibrated configuration. Additionally we can also profit from an
OpenMP parallelized version of the BLAS/LAPACK library for an additional speedup, which
also effects equilibration overhead NMPI × Tauto /NOMP , where NMPI is the number of cores and
NOMP the number of OpenMP threads. For a given number of independent measurements
7This corresponds to a Young tableau with single column and N /2 rows.

44

SciPost Phys. 3, 013 (2017)

Nmeas , we therefore need a wall-clock time given by


Tauto
Nmeas
T=
1+
.
NOMP
NMPI

(85)

As we typically have Nmeas /NMPI  1, the speedup is expected to be almost perfect, in accordance with the performance test results for the auxiliary field QMC code on SuperMUC [seeFig. 3 (a)].
For many problem sizes, 2 GB memory per MPI task (random walker) suffices such that we
typically start as many MPI tasks as there are physical cores per node. Due to the large amountof CPU time spent in MKL routines, we do not profit from the hyper-threading option. For large
systems, the memory requirement increases and this is tackled by increasing the amount of
OpenMP threads to decrease the stress on the memory system and to simultaneously reduce
the equilibration overhead [see Fig. 3 (b)]. For the displayed speedup, it was crucial to pin
the MPI tasks as well as the OpenMP threads in a pattern which keeps the threads as compact
as possible to profit from a shared cache. This also explains the drop in efficiency from 14 to
28 threads where the OpenMP threads are spread over both sockets.
We store the field configurations of the random walker as checkpoints, such that a longsimulation can be easily split into several short simulations. This procedure allows us to takeadvantage of chained jobs using the dependency chains provided by the batch system.
(b)

QMC simulation
ideal ∼ x

Performan e (1/T )Performan e (1/T )(a) 100

10

1
28

56

112

224

448

QMC simulation
ideal ∼ x

101
1

896 1792

Number of ores

2

4

7

14

28

Number of OpenMP threads

Figure 3: MPI (a) and OpenMP (b) scaling behavior of the auxiliary field QMC code of the ALF
project on SuperMUC (phase 2/Haswell nodes) at the LRZ in Munich. The MPI performancedata was normalized to 28 cores and was obtained using a problem size of Ndim = 400. Thisis a medium to small system size that is the least favorable in terms of MPI synchronization
effects. The OpenMP performance data was obtained using a problem size of Ndim = 1296.
Employing 2 and 4 OpenMP threads introduces some synchronization/management overhead
such that the per-core performance is slightly reduced, compared to the single thread efficiency.
Further increasing the amount of threads to 7 and 14 keeps the efficiency constant. The dropin performance of the 28 thread configuration is due to the architecture as the threads arenow spread over both sockets of the node. To obtain the above results, it was crucial to pin
the processes in a fashion that keeps the OpenMP threads as compact as possible.

6

Conclusions and Future Releases

In its present form, the auxiliary field QMC code of the ALF project allows to simulate a large
class of non-trivial models, both efficiently and at minimal programming cost. There are many

45

SciPost Phys. 3, 013 (2017)

possible extensions which deserve to be considered in future releases. The model Hamiltonians we have presented so far are imaginary-time independent. This however can be easily
generalized to imaginary-time dependent model Hamiltonians thus allowing, for example, toaccess entanglement properties of interacting fermionic systems [33–35, 93]. Generalizations
to include global moves are equally desirable. This is a prerequisite to play with recent ideas ofself-learning algorithms [94] so as to possibly avoid the issue of critical slowing down. Paralleltempering schemes are equally desirable, so as to possibly alleviate long autocorrelation times.
Most of the above has already been tested and will be available in the next major release of
the ALF package.
On the longer term, we foresee further possible developments. At present, the QMC codeof this package is restricted to discrete HS fields such that implementations of the long-range
Coulomb repulsion – as introduced in [28,65,66] – are not yet included. Extensions to continuous HS fields are certainly possible, but require an efficient upgrading scheme such as hybrid
molecular dynamics [45]. An implementation of the ground state projective QMC method
within ALF is equally desirable. Efforts in the above directions will be pursued on the longerterm.
As it stands, programming a new model certainly requires some detailed knowledge of thealgorithm. To facilitate access we hope to maintain an increasing number of model Hamiltonians in the ALF repository. A further step is to aim at cross compatibility with other major
projects, especially the ALPS [58] project.

Acknowledgments
We are very grateful to S. Beyl, M. Hohenadler, F. Parisen Toldin, M. Raczkowski, T. Sato,J. Schwab, Z. Wang, and M. Weber for constant support during the development of this project.
We equally thank G. Hager, M. Wittmann, and G. Wellein for useful discussions and support.
FFA would also like to thank T. Lang and Z. Y. Meng for developments of the auxiliary fieldcode as well as T. Grover for many discussions.
Funding information MB thanks the Bavarian Competence Network for Technical and Scientific High Performance Computing (KONWIHR) for financial support. FG and JH thankthe SFB-1170 for financial support under projects Z03 and C01. FFA thanks the DFG-funded
FOR1807 and FOR1346 for partial financial support. Part of the optimization of the code wascarried out during the Porting and Tuning Workshop 2016 offered by the Forschungszentrum
Jülich. Calculations to extensively test this package were carried out both on SuperMUC at
the Leibniz Supercomputing Centre and on JURECA [95] at the Jülich Supercomputing Centre.
We thank both institutions for generous allocation of computing time.

A

LicenseThe ALF code is provided as an open source software such that it is available to all and we hopethat it will be useful. If you benefit from this code we ask that you acknowledge the ALF collaboration as mentioned on our homepage alf.physik.uni-wuerzburg.de. The Git repository at
alf.physik.uni-wuerzburg.de gives us the tools to create a small but vibrant community aroundthe code and provides a suitable entry point for future contributors and future developments.
The homepage is also the place where the original source files can be found. With the coming public release it was necessary to add copyright headers to our source files. The Creative
Commons licenses are a good way to share our documentation and it is also well accepted by46

SciPost Phys. 3, 013 (2017)

publishers. Therefore this documentation is licensed to you under a CC-BY-SA license. Thismeans you can share it and redistribute it as long as you cite the original source and license
your changes under the same license. The details are in the file license.CCBYSA that you shouldhave received with this documentation. The source code itself is licensed under a GPL license to
keep the source as well as any future work in the community. To express our desire for a proper
attribution we decided to make this a visible part of the license. To that end we have exercised
the rights of section 7 of GPL version 3 and have amended the license terms with an additional
paragraph that expresses our wish that if an author has benefitted from this code that he/sheshould consider giving back a citation as specified on alf.physik.uni-wuerzburg.de. This is not
something that is meant to restrict your freedom of use, but something that we strongly expect
to be good scientific conduct. The original GPL license can be found in the file license.GPL andthe additional terms can be found in license.additional. In favour to our users, the ALF code
contains part of the lapack implementation version 3.6.1 from http://www.netlib.org/lapack.
Lapack is licensed under the modified BSD license whose full text can be found in license.BSD.
With that being said, we hope that the ALF code will prove to you to be a suitable and highperformance tool that enables you to perform quantum Monte Carlo studies of solid statemodels of unprecedented complexity.
The ALF project’s contributors.

References[1] R. Blankenbecler, D. J. Scalapino and R. L. Sugar, Monte Carlo calculations of coupledboson-fermion systems, Phys. Rev. D 24, 2278 (1981), doi:10.1103/PhysRevD.24.2278.
[2] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis and R. T. Scalettar,
Numerical study of the two-dimensional Hubbard model, Phys. Rev. B 40, 506-516 (1989),
doi:10.1103/PhysRevB.40.506.
[3] G. Sugiyama and S. E. Koonin, Auxiliary field Monte-Carlo for quantum many-body groundstates, Ann. Phys. 168, 1 (1986), doi:10.1016/0003-4916(86)90107-7.
[4] S. Sorella, S. Baroni, R. Car and M. Parrinello, A Novel Technique for the Simulation of Interacting Fermion Systems, Europhys. Lett. 8, 663 (1989), doi:10.1209/02955075/8/7/014.
[5] S. Duane, A. D. Kennedy, B. J. Pendleton and D. Roweth, Hybrid Monte Carlo, Phys. Lett.
B 195, 216 (1987), doi:10.1016/0370-2693(87)91197-X.
[6] A. D. Kennedy, Algorithms for Dynamical Fermions, arXiv:hep-lat/0607038v2.
[7] F. F. Assaad and H. G. Evertz, World-line and determinantal quantum Monte Carlo methods for spins, phonons and electrons, in H. Fehske, R. Schneider and A. Weiße (eds.),
Computational Many-Particle Physics, vol. 739 of Lecture Notes in Physics, pp. 277–356.
Springer, Berlin Heidelberg, ISBN 978-3-540-74685-0, (2008) doi:10.1007/978-3-54074686-7_10.
[8] J. Gubernatis, N. Kawashima and P. Werner, Quantum Monte Carlo Methods, Cambridge University Press, Cambridge, ISBN 9780511902581 (2016),
doi:10.1017/CBO9780511902581.

47

SciPost Phys. 3, 013 (2017)[9] D. J. Scalapino, Numerical Studies of the 2D Hubbard Model, pp. 495–526, Springer
New York, New York, NY, ISBN 978-0-387-68734-6, doi:10.1007/978-0-387-68734-6_13(2007).
[10] J. P. F. LeBlanc et al., Solutions of the two-dimensional Hubbard model: Benchmarks
and results from a wide range of numerical algorithms, Phys. Rev. X 5, 041041 (2015),
doi:10.1103/PhysRevX.5.041041.
[11] M. Hohenadler, T. C. Lang and F. F. Assaad, Correlation effects in quantum spinHall insulators: A quantum Monte Carlo study, Phys. Rev. Lett. 106, 100403 (2011),
doi:10.1103/PhysRevLett.106.100403.
[12] D. Zheng, G.-M. Zhang and C. Wu, Particle-hole symmetry and interaction
effects in the Kane-Mele-Hubbard model, Phys. Rev. B 84, 205121 (2011),
doi:10.1103/PhysRevB.84.205121.
[13] F. F. Assaad and I. F. Herbut, Pinning the order: The nature of quantum criticality in the Hubbard model on honeycomb lattice, Phys. Rev. X 3, 031010 (2013),
doi:10.1103/PhysRevX.3.031010.
[14] F. Parisen Toldin, M. Hohenadler, F. F. Assaad and I. F. Herbut, Fermionic quantum criticality in honeycomb and π-flux Hubbard models: Finite-size scaling of renormalizationgroup-invariant observables from quantum Monte Carlo, Phys. Rev. B 91, 165108 (2015),
doi:10.1103/PhysRevB.91.165108.
[15] Y. Otsuka, S. Yunoki and S. Sorella, Universal quantum criticality in the metal-insulatortransition of two-dimensional interacting Dirac electrons, Phys. Rev. X 6, 011029 (2016),
doi:10.1103/PhysRevX.6.011029.
[16] S. Chandrasekharan and A. Li, Quantum
mensional lattice Gross-Neveu models, Phys.
doi:10.1103/PhysRevD.88.021701.

critical behavior in three diRev. D 88, 021701 (2013),[17] V. Ayyar and S. Chandrasekharan, Massive fermions without fermion bilinear condensates,
Phys. Rev. D 91, 065035 (2015), doi:10.1103/PhysRevD.91.065035.
[18] Z.-X. Li, Y.-F. Jiang, S.-K. Jian and H. Yao, Fermion-induced quantum critical points: type-II
Landau-forbidden transitions, arXiv:1512.07908.
[19] F. F. Assaad and T. Grover, Simple fermionic model of deconfined phases and phase transitions, Phys. Rev. X 6, 041049 (2016), doi:10.1103/PhysRevX.6.041049.
[20] F. F. Assaad, Quantum Monte Carlo simulations of the half-filled two-dimensional Kondo
lattice model., Phys. Rev. Lett. 83, 796 (1999), doi:10.1103/physrevlett.83.796.
[21] S. Capponi and F. F. Assaad, Spin and charge dynamics of the ferromagnetic and antiferromagnetic two-dimensional half-filled Kondo lattice model, Phys. Rev. B 63, 155114 (2001),
doi:10.1103/physrevb.63.155114.
[22] Y. Schattner, S. Lederer, S. A. Kivelson and E. Berg, Ising nematic quantum critical point in a metal: A Monte Carlo study, Phys. Rev. X 6, 031028 (2016),
doi:10.1103/PhysRevX.6.031028.
[23] X. Y. Xu, K. Sun, Y. Schattner, E. Berg and Z. Y. Meng, Non-Fermi-liquid at (2+1)d ferromagnetic quantum critical point, arXiv:1612.06075.

48

SciPost Phys. 3, 013 (2017)[24] E. Berg, M. A. Metlitski and S. Sachdev, Sign-problem–free quantum Monte
Carlo of the onset of antiferromagnetism in metals, Science 338, 1606 (2012),
doi:10.1126/science.1227769.
[25] H.-K. Tang, X. Yang, J. Sun and H.-Q. Lin, Berezinskii-Kosterlitz-Thoules phase transition of
spin-orbit coupled Fermi gas in optical lattice, EPL 107, 40003 (2014), doi:10.1209/02955075/107/40003.
[26] F. F. Assaad, Phase diagram of the half-filled two-dimensional SU(n) HubbardHeisenberg model: A quantum Monte Carlo study, Phys. Rev. B 71, 075103 (2005),
doi:10.1103/PhysRevB.71.075103.
[27] T. C. Lang, Z. Y. Meng, A. Muramatsu, S. Wessel and F. F. Assaad, Dimerized solids andresonating plaquette order in SU(n)-Dirac fermions, Phys. Rev. Lett. 111, 066401 (2013),
doi:10.1103/PhysRevLett.111.066401.
[28] M. Hohenadler, F. Parisen Toldin, I. F. Herbut and F. F. Assaad, Phase diagram of the Kane-Mele-Coulomb model, Phys. Rev. B 90, 085146 (2014),
doi:10.1103/PhysRevB.90.085146.
[29] H.-K. Tang, E. Laksono, J. N. B. Rodrigues, P. Sengupta, F. F. Assaad and S. Adam,
Interaction-driven metal-insulator transition in strained graphene, Phys. Rev. Lett. 115,
186602 (2015), doi:10.1103/PhysRevLett.115.186602.
[30] M. Rigol, A. Muramatsu, G. G. Batrouni and R. T. Scalettar, Local quantum criticality in confined fermions on optical lattices, Phys. Rev. Lett. 91, 130403 (2003),
doi:10.1103/PhysRevLett.91.130403.
[31] D. Lee, Lattice simulations for few- and many-body systems, Prog. Part. Nucl. Phys. 63,
117 (2009), doi:10.1016/j.ppnp.2008.12.001.
[32] T. Grover, Entanglement of interacting fermions in quantum Monte Carlo calculations, Phys.
Rev. Lett. 111, 130402 (2013), doi:10.1103/PhysRevLett.111.130402.
[33] P. Broecker and S. Trebst, Rényi entropies of interacting fermions from determinantal
quantum Monte Carlo simulations, J. Stat. Mech. P08015 (2014), doi:10.1088/17425468/2014/08/p08015.
[34] F. F. Assaad, T. C. Lang and F. Parisen Toldin, Entanglement spectra of interacting fermions in quantum Monte Carlo simulations, Phys. Rev. B 89, 125121 (2014),
doi:10.1103/PhysRevB.89.125121.
[35] F. F. Assaad, Stable quantum Monte Carlo simulations for entanglement spectra of interacting fermions, Phys. Rev. B 91, 125146 (2015), doi:10.1103/PhysRevB.91.125146.
[36] P. Broecker and S. Trebst, Numerical stabilization of entanglement computation in
auxiliary-field quantum Monte Carlo simulations of interacting many-fermion systems,
Phys. Rev. E 94, 063306 (2016), doi:10.1103/PhysRevE.94.063306.
[37] C. Wu and S.-C. Zhang, Sufficient condition for absence of the sign problem inthe fermionic quantum Monte Carlo algorithm, Phys. Rev. B 71, 155115 (2005),
doi:10.1103/PhysRevB.71.155115.
[38] E. F. Huffman and S. Chandrasekharan, Solution to sign problems in halffilled spin-polarized electronic systems, Phys. Rev. B 89, 111101 (2014),
doi:10.1103/PhysRevB.89.111101.
49

SciPost Phys. 3, 013 (2017)[39] Z.-X. Li, Y.-F. Jiang and H. Yao, Solving the fermion sign problem in quantum Monte
Carlo simulations by Majorana representation, Phys. Rev. B 91, 241117 (2015),
doi:10.1103/PhysRevB.91.241117.
[40] Z. C. Wei, C. Wu, Y. Li, S. Zhang and T. Xiang, Majorana positivity and the fermion
sign problem of quantum Monte Carlo simulations, Phys. Rev. Lett. 116, 250601 (2016),
doi:10.1103/PhysRevLett.116.250601.
[41] Z.-X. Li, Y.-F. Jiang and H. Yao, Majorana-time-reversal symmetries: A fundamental principle for sign-problem-free quantum Monte Carlo simulations, Phys. Rev. Lett. 117, 267002(2016), doi:10.1103/PhysRevLett.117.267002.
[42] J. Hubbard, Calculation of partition functions, Phys. Rev. Lett. 3, 77 (1959),
doi:10.1103/PhysRevLett.3.77.
[43] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino and R. L. Sugar,
Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301
(1990), doi:10.1103/PhysRevB.41.9301.
[44] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations
to fermionic quantum Monte Carlo simulations, Phys. Rev. Lett. 94, 170201 (2005),
doi:10.1103/PhysRevLett.94.170201.
[45] S. Duane and J. B. Kogut, Hybrid stochastic differential equations applied to quantum
chromodynamics, Phys. Rev. Lett. 55, 2774 (1985), doi:10.1103/PhysRevLett.55.2774.
[46] J. E. Hirsch, Discrete Hubbard-Stratonovich transformation for fermion lattice models,
Phys. Rev. B 28, 4059 (1983), doi:10.1103/PhysRevB.28.4059.
[47] A. D. Sokal, Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms, in C. DeWitt-Morette, P. Cartier and A. Folacci (eds.), Functional Integration,
NATO ASI Series (Series B: Physics) 361, Springer, Boston, MA, doi:10.1007/978-14899-0319-8_6.
[48] H. G. Evertz, G. Lana and M. Marcu, Cluster algorithm for vertex models, Phys. Rev. Lett.
70, 875 (1993), doi:10.1103/PhysRevLett.70.875.
[49] A. W. Sandvik, Stochastic series expansion method with operator-loop update, Phys. Rev. B
59, R14157 (1999), doi:10.1103/PhysRevB.59.R14157.
[50] O. F. Syljuåsen and A. W. Sandvik, Quantum Monte Carlo with directed loops, Phys. Rev. E
66, 046701 (2002), doi:10.1103/PhysRevE.66.046701.
[51] J. E. Hirsch and R. M. Fye, Monte Carlo method for magnetic impurities in metals, Phys.
Rev. Lett. 56, 2521 (1986), doi:10.1103/PhysRevLett.56.2521.
[52] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer and P. Werner, Continuoustime Monte Carlo methods for quantum impurity models, Rev. Mod. Phys. 83, 349 (2011),
doi:10.1103/RevModPhys.83.349.
[53] F. F. Assaad, Continuous-time QMC Solvers for Electronic Systems in Fermionic and BosonicBaths, in E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein (eds.), DMFT at 25:
Infinite Dimensions: Lecture Notes of the Autumn School on Correlated Electrons 4, Verlag
des Forschungszentrum Jülich, Jülich, ISBN 978-3-89336-953-9 (2014).

50

SciPost Phys. 3, 013 (2017)[54] F. F. Assaad and T. C. Lang, Diagrammatic determinantal quantum Monte Carlo methods: Projective schemes and applications to the Hubbard-Holstein model, Phys. Rev. B 76,
035116 (2007), doi:10.1103/PhysRevB.76.035116.
[55] R. T. Scalettar, D. J. Scalapino and R. L. Sugar, New algorithm for the numerical simulationof fermions, Phys. Rev. B 34, 7911 (1986), doi:10.1103/PhysRevB.34.7911.
[56] S. Durr et al., Ab Initio Determination of Light Hadron Masses, Science 322, 1224 (2008),
doi:10.1126/science.1163233.
[57] F. F. Assaad, Quantum Monte Carlo methods on lattices: The determinantal method, in
J. Grotendorst, D. Marx and A. Muramatsu (eds.), Lecture notes of the Winter School on
Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms. 10,
pp. 99–155. Publication Series of the John von Neumann Institute for Computing, Jülich
(2002).
[58] B. Bauer et al., The ALPS project release 2.0: open source software for strongly correlated
systems, J. Stat. Mech. P05001 (2011), doi:10.1088/1742-5468/2011/05/p05001.
[59] Quest: Quantum electron simulation toolbox, http://quest.ucdavis.edu, accessed: 201706-14.
[60] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio and P. Seth, TRIQS:
A toolbox for research on interacting quantum systems, Comput. Phys. Commun. 196, 398(2015), doi:10.1016/j.cpc.2015.04.023.
[61] w2dynamics,

https://w2dynamics.physik.uni-wuerzburg.de/W2Dynamics/

W2Dynamics, accessed: 2017-06-14.
[62] L. Huang, Y. Wang, Z. Y. Meng, L. Du, P. Werner and X. Dai, iQIST: An open source
continuous-time quantum Monte Carlo impurity solver toolkit, Comput. Phys. Commun.
195, 140 (2015), doi:10.1016/j.cpc.2015.04.020.
[63] Y. Motome and M. Imada, A Quantum Monte Carlo Method and Its Applications to MultiOrbital Hubbard Models, J. Phys. Soc. Jpn. 66, 1872 (1997), doi:10.1143/JPSJ.66.1872.
[64] F. F. Assaad, M. Imada and D. J. Scalapino, Charge and spin structures of a d x 2 − y2 superconductor in the proximity of an antiferromagnetic Mott insulator, Phys. Rev. B 56, 15001(1997), doi:10.1103/PhysRevB.56.15001.
[65] R. Brower, C. Rebbi and D. Schaich, Hybrid Monte Carlo simulation on the graphene hexagonal lattice, arXiv:1204.5424.
[66] M. V. Ulybyshev, P. V. Buividovich, M. I. Katsnelson and M. I. Polikarpov, Monte Carlostudy of the semimetal-insulator phase transition in monolayer graphene with arealistic interelectron interaction potential, Phys. Rev. Lett. 111, 056801 (2013),
doi:10.1103/PhysRevLett.111.056801.
[67] M. Roger, J. H. Hetherington and J. M. Delrieu, Magnetism in solid 3 he, Rev. Mod. Phys.
55, 1 (1983), doi:10.1103/RevModPhys.55.1.
[68] J. Werner and F. F. Assaad, Ring-exchange periodic Anderson model for 3 He bilayers, Phys.
Rev. B 90, 205122 (2014), doi:10.1103/PhysRevB.90.205122.
[69] P. Werner and A. J. Millis, Efficient dynamical mean field simulation of the Holstein-Hubbard
model, Phys. Rev. Lett. 99, 146404 (2007), doi:10.1103/PhysRevLett.99.146404.
51

SciPost Phys. 3, 013 (2017)[70] R. M. Fye, New results on Trotter-like approximations, Phys. Rev. B 33, 6271 (1986),
doi:10.1103/PhysRevB.33.6271.
[71] M. Iazzi and M. Troyer, Efficient continuous-time quantum Monte Carlo algorithm for fermionic lattice models, Phys. Rev. B 91, 241118 (2015),
doi:10.1103/PhysRevB.91.241118.
[72] S. M. A. Rombouts, K. Heyde and N. Jachowicz, Quantum Monte Carlo method
for fermions, free of discretization errors, Phys. Rev. Lett. 82, 4155 (1999),
doi:10.1103/PhysRevLett.82.4155.
[73] E. Gull, P. Werner, O. Parcollet and M. Troyer, Continuous-time auxiliary-field Monte
Carlo for quantum impurity models, EPL 82, 57003 (2008), doi:10.1209/02955075/82/57003.
[74] S. Rombouts, K. Heyde and N. Jachowicz, A discrete Hubbard-Stratonovich decomposition for general, fermionic two-body interactions, Phys. Lett. A 242, 271 (1998),
doi:10.1016/S0375-9601(98)00197-2.
[75] D. Rost, E. V. Gorelik, F. Assaad and N. Blümer, Momentum-dependent pseudogaps
in the half-filled two-dimensional Hubbard model, Phys. Rev. B 86, 155109 (2012),
doi:10.1103/PhysRevB.86.155109.
[76] D. Rost, F. F. Assaad and N. Blümer, Quasi-continuous-time impurity solver for the dynamical mean-field theory with linear scaling in the inverse temperature, Phys. Rev. E 87,053305 (2013), doi:10.1103/PhysRevE.87.053305.
[77] N. Blümer, Multigrid Hirsch-Fye quantum Monte Carlo method for dynamical mean-field
theory, arXiv:0801.1222.
[78] Z. Bai, C. Lee, R.-C. Li and S. Xu, Stable solutions of linear systems involving long chain of matrix multiplications, Lin. Algebra Appl. 435, 659 (2011),
doi:10.1016/j.laa.2010.06.023.
[79] A. van der Sluis, Condition numbers and equilibration of matrices, Numer. Math. 14, 14(1969), doi:10.1007/BF02165096.
[80] J. W. Negele and H. Orland, Quantum many-particle systems, Frontiers in physics,
Addison-Wesley, Redwood City, California. (1988).
[81] W. Krauth, Statistical Mechanics: Algorithms and Computations, Oxford University Press,
ISBN 9780198515364 (2006).
[82] C. J. Geyer, Practical Markov chain Monte Carlo, Statistical Science 7(4), 473 (1992).
[83] R. M. Neal, Probabilistic inference using Markov chain Monte Carlo methods, (1993).
[84] M. Bercx, J. S. Hofmann, F. F. Assaad and T. C. Lang, Spontaneous particle-hole symmetry breaking of correlated fermions on the Lieb lattice, Phys. Rev. B 95, 035108 (2017),
doi:10.1103/PhysRevB.95.035108.
[85] B. Efron and C. Stein, The jackknife estimate of variance, Ann. Statist. 9, 586 (1981),
doi:10.1214/aos/1176345462.
[86] S. Chakravarty, B. I. Halperin and D. R. Nelson, Low-temperature behavior of
two-dimensional quantum antiferromagnets, Phys. Rev. Lett. 60, 1057 (1988),
doi:10.1103/PhysRevLett.60.1057.
52

SciPost Phys. 3, 013 (2017)[87] M. B. Thompson, A Comparison of Methods for Computing Autocorrelation Time,
arXiv:1011.0175.
[88] I. Milat, F. F. Assaad and M. Sigrist, Field induced magnetic ordering transition in Kondo
insulators., Eur. Phys. J. B 38, 571 (2004), doi:10.1140/epjb/e2004-00154-5.
[89] M. Bercx, T. C. Lang and F. F. Assaad, Magnetic field induced semimetal-to-cantedantiferromagnet transition on the honeycomb lattice, Phys. Rev. B 80, 045412 (2009),
doi:10.1103/PhysRevB.80.045412.
[90] X. Y. Xu, K. S. D. Beach, K. Sun, F. F. Assaad and Z. Y. Meng, Topological phase transitions with SO(4) symmetry in (2+1)d interacting Dirac fermions, Phys. Rev. B 95, 085110(2017), doi:10.1103/PhysRevB.95.085110.
[91] K. S. D. Beach, P. A. Lee and P. Monthoux, Field-induced antiferromagnetism in the Kondo
insulator, Phys. Rev. Lett. 92, 026401 (2004), doi:10.1103/physrevlett.92.026401.
[92] Z.-X. Li, Y.-F. Jiang and H. Yao, Fermion-sign-free Majorana-quantum-Monte-Carlo studies
of quantum critical phenomena of Dirac fermions in two dimensions, New J. Phys. 17,085003 (2015), doi:10.1088/1367-2630/17/8/085003.
[93] F. F. Assaad, Interface superconductivity: Get it strained, Nat. Phys. 10, 905 (2014),
doi:10.1038/nphys3161.
[94] X. Y. Xu, Y. Qi, J. Liu, L. Fu and Z. Y. Meng, Self-learning determinantal quantum Monte
Carlo method, Phys. Rev. B 96, 041119 (2017), doi:10.1103/PhysRevB.96.041119.
[95] Jülich Supercomputing Centre, JURECA: General-purpose supercomputer at Jülich Supercomputing Centre, JLSRF 2, A62 (2016), doi:10.17815/jlsrf-2-121.

53