DCA++: A software framework to solve correlated electron problems with   modern quantum cluster methods

Urs R. H\"ahner; Gonzalo Alvarez; Thomas A. Maier; Raffaele Solc\`a,; Peter Staar; Michael S. Summers; Thomas C. Schulthess

arXiv:1901.01047·physics.comp-ph·October 30, 2019·Comput. Phys. Commun.

DCA++: A software framework to solve correlated electron problems with modern quantum cluster methods

Urs R. H\"ahner, Gonzalo Alvarez, Thomas A. Maier, Raffaele Solc\`a,, Peter Staar, Michael S. Summers, Thomas C. Schulthess

PDF

TL;DR

DCA++ is an open-source, high-performance software framework that implements advanced quantum cluster algorithms to solve strongly correlated electron problems, scalable across various modern computing architectures.

Contribution

It introduces a scalable, portable implementation of the dynamical cluster approximation and its extension, with sustainable development practices for quantum many-body simulations.

Findings

01

Achieves efficient scalability from workstations to supercomputers

02

Provides portable performance on CPU-GPU hybrid architectures

03

Supports insights into high-temperature superconductivity

Abstract

We present the first open release of the DCA++ project, a high-performance research software framework to solve quantum many-body problems with cutting edge quantum cluster algorithms. DCA++ implements the dynamical cluster approximation (DCA) and its DCA $^{+}$ extension with a continuous self-energy. The algorithms capture nonlocal correlations in strongly correlated electron systems, thereby giving insight into high- $T_{c}$ superconductivity. The code's scalability allows efficient usage of systems at all scales, from workstations to leadership computers. With regard to the increasing heterogeneity of modern computing machines, DCA++ provides portable performance on conventional and emerging new architectures, such as hybrid CPU-GPU, sustaining multiple petaflops on ORNL's Titan and CSCS' Piz Daint supercomputers. Moreover, we show how sustainable and scalable development of the code base…

Figures19

Click any figure to enlarge with its caption.

Tables1

Table 1. Table 1 : Relevant CMake options for building DCA++.

Option	Default value	Description
CMAKE_BUILD_TYPE		The build type. Options are: [empty], Debug, Release, RelWithDebInfo, MinSizeRel. According to the build type CMake sets different optimization and debug compiler flags. Release is recommended for production runs.
FFTW_INCLUDE_DIR		Path to fftw3.h.
FFTW_LIBRARY		The FFTW3(-compatible) library. If the -D flag is used, multiple libraries have to be separated by semicolons and enclosed by double quotes.
DCA_HAVE_LAPACK		Set to TRUE if you use a compiler wrapper that automatically links against BLAS and LAPACK libraries, to prevent CMake from searching for these libraries.
DCA_WITH_MPI	ON	Enable MPI. If MPI is enabled, set the environment variable CXX to the MPI compiler wrapper before running CMake.
DCA_WITH_CUDA	OFF	Enable GPU support.
DCA_WITH_PINNED_HOST_MEMORY	OFF	Enable pinned host memory. (advanced)
CUDA_GPU_ARCH	sm_60	Name of the real architecture to build for.
CUDA_TOOLKIT_ROOT_DIR		Path to the CUDA Toolkit. Determined by CMake. Set it manually, if CMake cannot find it.
MAGMA_DIR		Path to the MAGMA installation directory. Hint for CMake to find MAGMA.
DCA_BUILD_DCA	ON	Build main_dca.
DCA_BUILD_ANALYSIS	ON	Build main_analysis.
DCA_BUILD_CLUSTER_SOLVER_CHECK	OFF	Build cluster_solver_check, an application that tests a QMC solver against exact diagonalization.
DCA_WITH_TESTS_FAST	OFF	Build DCA++’s fast tests.
DCA_WITH_TESTS_EXTENSIVE	OFF	Build DCA++’s extensive tests.
DCA_WITH_TESTS_PERFORMANCE	OFF	Build DCA++’s performance tests. (Only in Release mode.)
TEST_RUNNER		Command for executing (MPI) programs. Required if MPI is enabled or on systems that use a command for launching executables (e.g. aprun or srun).
MPIEXEC_NUMPROC_FLAG	-n	Flag used by TEST_RUNNER to specify the number of processes.
MPIEXEC_PREFLAGS		Flags to pass to TEST_RUNNER directly before the executable to run.
DCA_CLUSTER_SOLVER	CT-AUX	The cluster solver for the DCA/DCA⁺ calculation. Options are: CT-AUX, SS-CT-HYB.
DCA_WITH_THREADED_SOLVER	ON	Use multiple walker and accumulator threads in the cluster solver.
DCA_LATTICE	square	Lattice type, options are: bilayer, square, triangular.
DCA_POINT_GROUP	D4	Point group symmetry, options are: C6, D4.
DCA_MODEL	tight-binding	Model type, options are: tight-binding.
DCA_RNG	std::mt19937_64	Random number generator, options are: std::mt19937_64, std::ranlux48, custom. See Random number generator section in the Wiki for details.
DCA_RNG_CLASS		Class name including namespaces of the custom random number generator.
DCA_RNG_HEADER		Header file of the custom random number generator.
DCA_RNG_LIBRARY		Custom random number generator library.
DCA_PROFILER	None	Profiler type, options are: None, Counting, PAPI.
DCA_WITH_AUTOTUNING	OFF	Enable auto-tuning. Requires a profiler type other than None. (advanced)
DCA_WITH_GNUPLOT	OFF	Enable Gnuplot.
DCA_WITH_SINGLE_PRECISION_MEASUREMENTS	OFF	Measure in single precision.
DCA_WITH_SINGLE_PRECISION_COARSEGRAINING	OFF	Coarsegrain in single precision. (advanced)
DCA_WITH_QMC_BIT	OFF	Enable QMC solver built-in tests. (advanced)

Equations99

\phi_{\bf K}({\bf k})=\left\{\begin{array}[]{ll}1,&{\bf k}\in{\bf K}^{\mathrm{th}}\,\mathrm{patch},\\ 0,&\mathrm{otherwise}.\end{array}\right.

\phi_{\bf K}({\bf k})=\left\{\begin{array}[]{ll}1,&{\bf k}\in{\bf K}^{\mathrm{th}}\,\mathrm{patch},\\ 0,&\mathrm{otherwise}.\end{array}\right.

\frac{N _{c}}{V _{BZ}} \int_{BZ} d k ϕ_{K} (k) ϕ_{K^{'}} (k) = δ_{K, K^{'}},

\frac{N _{c}}{V _{BZ}} \int_{BZ} d k ϕ_{K} (k) ϕ_{K^{'}} (k) = δ_{K, K^{'}},

ϕ_{K} (k) = ϕ (k - K),

ϕ_{K} (k) = ϕ (k - K),

ϕ (k - K) = ϕ (K - k) .

ϕ (k - K) = ϕ (K - k) .

Σ (k = K + \tilde{k}, i ω_{n}) \approx Σ_{c} (K, i ω_{n}) .

Σ (k = K + \tilde{k}, i ω_{n}) \approx Σ_{c} (K, i ω_{n}) .

Σ_{DCA} (k, i ω_{n}) \equiv K \sum ϕ_{K} (k) Σ_{c} (K, i ω_{n}) .

Σ_{DCA} (k, i ω_{n}) \equiv K \sum ϕ_{K} (k) Σ_{c} (K, i ω_{n}) .

G_{c} (K, i ω_{n}) \to \overset{ˉ}{G} (K, i ω_{n}) = \frac{N _{c}}{V _{BZ}} \int_{BZ} d k ϕ_{K} (k) [G_{0}^{- 1} (k, i ω_{n}) - Σ_{DCA} (k, i ω_{n})]^{- 1} .

G_{c} (K, i ω_{n}) \to \overset{ˉ}{G} (K, i ω_{n}) = \frac{N _{c}}{V _{BZ}} \int_{BZ} d k ϕ_{K} (k) [G_{0}^{- 1} (k, i ω_{n}) - Σ_{DCA} (k, i ω_{n})]^{- 1} .

G_{0} (k, i ω_{n}) = [i ω_{n} - ε_{k} + μ]^{- 1} .

G_{0} (k, i ω_{n}) = [i ω_{n} - ε_{k} + μ]^{- 1} .

G_{0} (K, i ω_{n}) = [\overset{ˉ}{G}^{- 1} (K, i ω_{n}) + Σ_{c} (K, i ω_{n})]^{- 1} .

G_{0} (K, i ω_{n}) = [\overset{ˉ}{G}^{- 1} (K, i ω_{n}) + Σ_{c} (K, i ω_{n})]^{- 1} .

{G_{0} (K, i ω_{n}), H_{int}} ⟶ QMC Σ_{c} (K, i ω_{n}) .

{G_{0} (K, i ω_{n}), H_{int}} ⟶ QMC Σ_{c} (K, i ω_{n}) .

∣ Σ_{c}^{(n)} - Σ_{c}^{(n - 1)} ∣ < ε_{QMC} .

∣ Σ_{c}^{(n)} - Σ_{c}^{(n - 1)} ∣ < ε_{QMC} .

Σ_{c}^{(n)} = α Σ_{QMC} + (1 - α) Σ_{c}^{(n - 1)}, 0 < α \leq 1 .

Σ_{c}^{(n)} = α Σ_{QMC} + (1 - α) Σ_{c}^{(n - 1)}, 0 < α \leq 1 .

Σ_{c} (K, i ω_{n}) = \frac{N _{c}}{V _{BZ}} \int_{BZ} d k ϕ_{K} (k) Σ_{DCA^{+}} (k, i ω_{n}) .

Σ_{c} (K, i ω_{n}) = \frac{N _{c}}{V _{BZ}} \int_{BZ} d k ϕ_{K} (k) Σ_{DCA^{+}} (k, i ω_{n}) .

\displaystyle\big{\{}{\bf K},\Sigma_{c}({\bf K},{i\omega_{n}})\big{\}}

\displaystyle\big{\{}{\bf K},\Sigma_{c}({\bf K},{i\omega_{n}})\big{\}}

with \tilde{Σ} (K, i ω_{n})

\tilde{Σ} (k, i ω_{n}) = \frac{N _{c}}{V _{BZ}} \int_{BZ} d k^{'} ϕ (k - k^{'}) Σ_{DCA^{+}} (k^{'}, i ω_{n}),

\tilde{Σ} (k, i ω_{n}) = \frac{N _{c}}{V _{BZ}} \int_{BZ} d k^{'} ϕ (k - k^{'}) Σ_{DCA^{+}} (k^{'}, i ω_{n}),

\overset{ˉ}{G} (K, i ω_{n})

\overset{ˉ}{G} (K, i ω_{n})

\overset{ˉ}{Σ} (K, i ω_{n})

G_{0} (K, i ω_{n}) = [\overset{ˉ}{G}^{- 1} (K, i ω_{n}) + \overset{ˉ}{Σ} (K, i ω_{n})]^{- 1} .

G_{0} (K, i ω_{n}) = [\overset{ˉ}{G}^{- 1} (K, i ω_{n}) + \overset{ˉ}{Σ} (K, i ω_{n})]^{- 1} .

∣ Σ_{c}^{(n)} - \overset{ˉ}{Σ}^{(n - 1)} ∣ < ε_{QMC} .

∣ Σ_{c}^{(n)} - \overset{ˉ}{Σ}^{(n - 1)} ∣ < ε_{QMC} .

H

H

H_{0}

H_{int}

Z = e^{- K} Tr [e^{- β H_{0}} T_{τ} e^{- \int_{0}^{β} d τ (H_{int} - K / β)}] .

Z = e^{- K} Tr [e^{- β H_{0}} T_{τ} e^{- \int_{0}^{β} d τ (H_{int} - K / β)}] .

Z = e^{- K} k = 0 \sum \infty (\frac{K}{β})^{k} \int_{0}^{β} d τ_{1} \dots \int_{t_{k - 1}}^{β} d τ_{k} Tr [e^{- (β - τ_{k} + τ_{1}) H_{0}} (1 - \frac{β H _{int} ( τ _{k} )}{K}) \dots e^{- (τ_{2} - τ_{1}) H_{0}} (1 - \frac{β H _{int} ( τ _{1} )}{K})] .

Z = e^{- K} k = 0 \sum \infty (\frac{K}{β})^{k} \int_{0}^{β} d τ_{1} \dots \int_{t_{k - 1}}^{β} d τ_{k} Tr [e^{- (β - τ_{k} + τ_{1}) H_{0}} (1 - \frac{β H _{int} ( τ _{k} )}{K}) \dots e^{- (τ_{2} - τ_{1}) H_{0}} (1 - \frac{β H _{int} ( τ _{1} )}{K})] .

1 - \frac{β H _{int}}{K}

1 - \frac{β H _{int}}{K}

cosh (γ_{μν})

Z

Z

Z_{k} ({μ, ν, τ, s}_{k})

⟨ G (k, i ω_{n})⟩ = G_{0} (k, i ω_{n}) - \frac{1}{β} G_{0} (k, i ω_{n}) ⟨ M (k, i ω_{n})⟩ G_{0} (k, i ω_{n}) .

⟨ G (k, i ω_{n})⟩ = G_{0} (k, i ω_{n}) - \frac{1}{β} G_{0} (k, i ω_{n}) ⟨ M (k, i ω_{n})⟩ G_{0} (k, i ω_{n}) .

Σ (k, i ω_{n}) = G_{0}^{- 1} (k, i ω_{n}) - ⟨ G (k, i ω_{n}) ⟩^{- 1} .

Σ (k, i ω_{n}) = G_{0}^{- 1} (k, i ω_{n}) - ⟨ G (k, i ω_{n}) ⟩^{- 1} .

f_{ω} φ_{ω} = \int_{0}^{β} d τ e^{iω τ} [\int_{0}^{β} d τ^{'} φ (τ - τ^{'}) f (τ^{'})] .

f_{ω} φ_{ω} = \int_{0}^{β} d τ e^{iω τ} [\int_{0}^{β} d τ^{'} φ (τ - τ^{'}) f (τ^{'})] .

f_{ω} \leftarrow \frac{1}{φ _{ω}} FFT [{\overset{ˉ}{f_{l}} = i \sum φ (τ_{i} - \frac{l β}{m N _{ω}}) f_{i} ∣ l \in 0, 1, \dots, m N_{ω} - 1}] .

f_{ω} \leftarrow \frac{1}{φ _{ω}} FFT [{\overset{ˉ}{f_{l}} = i \sum φ (τ_{i} - \frac{l β}{m N _{ω}}) f_{i} ∣ l \in 0, 1, \dots, m N_{ω} - 1}] .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

DCA++: A software framework to solve correlated electron problems with modern quantum cluster methods

Urs R. Hähner

Gonzalo Alvarez

Thomas A. Maier

Raffaele Solcà

Peter Staar

Michael S. Summers

Thomas C. Schulthess

Institute for Theoretical Physics, ETH Zurich, 8093 Zurich, Switzerland

Computational Science and Engineering Division, Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA

Swiss National Supercomputing Centre, ETH Zurich, 6900 Lugano, Switzerland

IBM Research – Zurich, 8803 Rüschlikon, Switzerland

Abstract

We present the first open release of the DCA++ project, a high-performance research software framework to solve quantum many-body problems with cutting edge quantum cluster algorithms. DCA++ implements the dynamical cluster approximation (DCA) and its DCA+ extension with a continuous self-energy. The algorithms capture nonlocal correlations in strongly correlated electron systems, thereby giving insight into high- $T_{c}$ superconductivity. The code’s scalability allows efficient usage of systems at all scales, from workstations to leadership computers. With regard to the increasing heterogeneity of modern computing machines, DCA++ provides portable performance on conventional and emerging new architectures, such as hybrid CPU-GPU, sustaining multiple petaflops on ORNL’s Titan and CSCS’ Piz Daint supercomputers. Moreover, we show how sustainable and scalable development of the code base has been achieved by adopting standard techniques of the software industry. These include employing a distributed version control system, applying test-driven development and following continuous integration.

keywords:

Strongly correlated electron systems; Quantum cluster algorithms; Dynamical cluster approximation; Continuous-time quantum Monte Carlo; Extreme-scale computing; Sustainable software development

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: DCA++

Licensing provisions: BSD-3-Clause

Programming language: C++14 and CUDA

Nature of problem: Understanding the fascinating physics of strongly correlated electron systems requires the development of sophisticated algorithms and their implementation on leadership computing systems.

Solution method: The DCA++ code provides a highly scalable and efficient implementation of the dynamical cluster approximation (DCA) and its DCA+ extension.

1 Introduction

Computer simulations have become indispensable across all areas of science and physics in particular. In condensed matter physics significant computational effort is devoted to solving the quantum many-body problem of interacting electrons. Strong correlations between electrons are the origin of a variety of fascinating phenomena, most prominently high-temperature superconductivity. The enormous complexity of the quantum many-body problem - the Hilbert space grows exponentially with the system size - and lack of approximations and simplifications in the physically relevant regime pose a great challenge to date.

Progress in solving the many-electron problem has been made by introducing models that comprise the fundamental mechanisms that lead to the observed phenomena. The two-dimensional Hubbard model, for instance, is believed to capture the important physics in the superconducting cuprates [1, 2]. The Hubbard Hamiltonian describes fermions on a lattice, where the particles are allowed to hop between lattice sites and interact through Coulomb repulsion when they occupy the same site. Despite the simple structure, there exists no exact solution except for the one-dimensional case.

To get insight into the physics described by the Hubbard model, one has to resort to numerical methods. Exact diagonalization (ED) and a variety of quantum Monte Carlo (QMC) algorithms solve the model on a finite-size lattice. While ED is restricted to small systems due to the exponential scaling of the problem size with the number of lattice sites, the negative sign problem prevents QMC calculations on large lattices or at low temperatures. Finite-size effects arise from the truncation of the infinite lattice to a finite number of sites.

Mean-field theories choose a different approach and are formulated in the thermodynamic limit, that is on the infinite lattice. In dynamical mean-field theory (DMFT) [3] the infinite lattice problem becomes tractable by reducing it to the self-consistent solution of an effective impurity model. Underlying the single-site approximation is the assumption that the self-energy is local in real space, which holds in the limit of infinite dimensions. However, spatially nonlocal correlations are often essential to describe phase transitions, which is the reason why DMFT has been extended by a number of quantum cluster algorithms [4]. The dynamical cluster approximation (DCA) [4, 5, 6] maps the bulk lattice problem to a finite-size periodic cluster that is self-consistently embedded in a dynamical mean-field. The effective cluster problem is best solved by means of continuous-time QMC methods [7]. Compared to finite-size QMC calculations, DCA’s mean-field approach not only reduces the fermion sign problem [6], but also generally leads to results that converge faster to the thermodynamic limit [6, 8].

The DCA++ project provides an efficient and highly scalable C++ implementation of the DCA algorithm and the DCA*+* extension that introduces a continuous lattice self-energy [9]. The code has been developed by a team at ORNL and ETH Zurich, which, besides modern software design, has the focus on developing and implementing sophisticated algorithms [10, 11] and exploiting the benefits of emerging architectures, such as hybrid CPU-GPU systems. In addition to setting benchmarks for extreme scale scientific applications [12, 13], the DCA++ code has provided the numerical evidence for a series of scientific publications in the field of strongly correlated electron systems. Studies of the pairing mechanism in unconventional superconductors [14] and the influence of near neighbor Coulomb repulsion on $d$ -wave superconductivity [15], for example, were based on DCA calculations with DCA++.

In software-driven research, scientific productivity is strongly coupled to software productivity. Hence, the scientific output of a research group can be held by challenges such as new computer architectures, advanced algorithms or changing teams of developers. To prevent this productivity collapse, software development needs to be sustainable and scalable by producing comprehensible, maintainable, and extensible code. At the same time, it is essential to release changes to the code base, from new features to bug fixes, rapidly to the users, an ability that usually falls under the term continuous delivery. Last but not least, the scientific standard demands correctness, credibility and reproducibility of numerical results in published work. To fulfill these requirements in a challenging environment formed by complex algorithms, performance sensitive codes, the diversity of architectures, and the multidisciplinary of teams, the DCA++ project employs well-proven tools and successful techniques of the software industry [16]. While adopting these methods can require an effort, we believe that they represent a substantial factor for a research code to become a long-lived software project.

This paper accompanies the first open release of the DCA++ code and is structured as follows: Section 2 reviews the underlying DCA and DCA ${}^{+}\>$ algorithms, briefly discusses the continuous-time auxiliary-field and continuous-time hybridization expansion QMC algorithms, and outlines the computation of two-particle correlation functions. Section 3 gives an overview of the main components of the DCA++ code, addresses prerequisites and the CMake-based building routine, and explains how to run a DCA or DCA ${}^{+}\>$ calculation. Section 4 explores the software development practice that is followed in the DCA++ project, but whose applicability extends to many other collaborative scientific software efforts. Section 5 provides two typical use cases supplemented by performance analyses. The paper concludes with an outline of future development plans and a summary of ways to contribute back to the DCA++ project.

2 Formalism and methods

2.1 The dynamical cluster approximation

Like finite-size methods, the DCA algorithm replaces the infinite real-space lattice by a finite-size cluster to reduce the complexity of the problem. The reduction of the infinite lattice to a cluster of $N_{c}$ sites corresponds to a discretization of the first Brillouin zone into a set of $N_{c}$ cluster momenta ${\bf K}$ (see Fig. 1). In contrast to finite-size methods, which solve the cluster in isolation, the DCA embeds the cluster in a mean field to retain information about the degrees of freedom of the lattice that are not contained in the cluster. This is achieved through a coarse-graining procedure in momentum space, where these degrees of freedom are averaged out. For this purpose, the first Brillouin zone (BZ) is divided into $N_{c}$ patches (see Fig. 1). Each patch is centered around a unique cluster momentum ${\bf K}$ and represented by a patch function $\phi_{\bf K}({\bf k})$ , the characteristic function of the patch,

[TABLE]

In addition to the orthogonality condition,

[TABLE]

we require the patches to have equal size and shape,

[TABLE]

and to exhibit inversion symmetry,

[TABLE]

Despite these restrictions, the choice of patches is far from unique and it is even possible to create interleaved shapes [17]. In the standard coarse-graining the patches are defined as the Brillouin zones of the superlattice as illustrated in Fig. 1.

The underlying assumption of the DCA is that the lattice self-energy $\Sigma({\bf k},{i\omega_{n}})$ is a slowly varying function in momentum space and in the vicinity of a cluster momentum ${\bf K}$ well approximated by the self-energy of the cluster $\Sigma_{c}({\bf K},{i\omega_{n}})$ ,

[TABLE]

Based on this assumption, we can approximate the lattice self-energy $\Sigma({\bf k},{i\omega_{n}})$ by a piecewise constant continuation of the cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ by means of the patch functions $\phi_{\bf K}({\bf k})$ ,

[TABLE]

Unlike finite-size methods, the DCA, using above approximation of the lattice self-energy, replaces the Green’s function of the cluster $G_{c}({\bf K},{i\omega_{n}})$ with the coarse-grained average of the lattice Green’s function,

[TABLE]

$G_{0}({\bf k},{i\omega_{n}})$ is the non-interacting lattice Green’s function given by

[TABLE]

To formulate the effective cluster problem of the DCA, we need the corresponding bare propagator of the cluster $\mathcal{G}_{0}({\bf K},{i\omega_{n}})$ . We obtain it by reversing the Dyson equation in the cluster-exclusion step,

[TABLE]

$\mathcal{G}_{0}({\bf K},{i\omega_{n}})$ is fundamentally different from the free lattice propagator evaluated at the cluster momenta ${\bf K}$ , $G_{0}({\bf K},{i\omega_{n}})$ . It incorporates the coupling of the cluster to the mean-field that represents the remaining degrees of freedom of the lattice. The effective cluster problem is completed by the interacting part of the Hamiltonian $H_{\mathrm{int}}$ and best solved by QMC methods to calculate a new cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ ,

[TABLE]

This closes the DCA self-consistency loop.

To start the loop, one has to provide an initial cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ , which in practice is either zero or the result of a previous calculation for a slightly higher temperature (see Section 3.3). Steps (6), (7), (9) and (10) are then iterated until convergence is reached, where the criterion for convergence is dictated by the stochastic error of the QMC solver $\varepsilon_{\mathrm{QMC}}$ ,

[TABLE]

To stabilize convergence of the loop, the self-energy produced by the QMC solver $\Sigma_{\mathrm{QMC}}({\bf K},{i\omega_{n}})$ is often mixed with the previous value of the cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ ,

[TABLE]

The complete DCA algorithm is summarized in Fig. 2.

Relying on QMC methods, the DCA algorithm inherits their negative sign problem. Nevertheless, the sign problem turns out to be less severe than in finite-size lattice QMC simulations, an effect the dynamical mean-field has been credited with [6].

Like finite-size methods, the DCA yield the exact result of the infinite lattice problem as $N_{c}\rightarrow\infty$ . In other words, if clusters large enough to extrapolate to infinite cluster size are accessible, we can obtain exact results for the thermodynamic limit by finite-size scaling. In doing so, convergence with respect to cluster size indicates that the current size of the cluster captures all relevant nonlocal correlations. In the opposite limit, i.e. for $N_{c}=1$ , the DCA recovers DMFT.

2.2 The DCA ${}^{+}\>$ algorithm

The cluster shape dependence represents one of the major drawbacks of the DCA algorithm (see Fig. 4). The expansion of the lattice self-energy $\Sigma_{\mathrm{DCA}}({\bf k},{i\omega_{n}})$ , Eq. (6), evidently depends on the shape of the coarse-graining patches and is therefore sensitive to the definition of the cluster. A way to cure this problem is to abandon the representation of the lattice self-energy in terms of step functions and find a continuous description. Simple interpolation schemes of the cluster self-energy, however, have been shown to lead to causality violations [5]. Instead, we start from a coarse-graining equation for the self-energy, equivalent to the expression for the Green’s function in Eq. (7). Inverting Eq. (6) with the help of the orthogonality condition of the patch functions, Eq. (2), we obtain the relation,

[TABLE]

This equation implies that we need to generate a lattice self-energy $\Sigma_{\mathrm{DCA}^{+}}({\bf k},{i\omega_{n}})$ whose coarse-grained average equals the cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ . Note that the piecewise constant form of the standard DCA algorithm, Eq. (6), trivially satisfies this condition. Finding the lattice self-energy $\Sigma_{\mathrm{DCA^{+}}}({\bf k},{i\omega_{n}})$ in DCA ${}^{+}\>$ is a two step procedure. First, we interpolate the cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ between the cluster momenta ${\bf K}$ ,

[TABLE]

Next, we extend the coarse-graining condition, Eq. (13), to the interpolated lattice self-energy $\tilde{\Sigma}({\bf k},{i\omega_{n}})$ ,

[TABLE]

where we used the translation property, Eq. (3), and inversion symmetry, Eq. (4), of the patch function $\phi_{\bf K}({\bf k})$ . According to Eq. (15), $\tilde{\Sigma}({\bf k},{i\omega_{n}})$ is related to $\Sigma_{\mathrm{DCA^{+}}}({\bf k},{i\omega_{n}})$ by a convolution with the patch function $\phi({\bf k}-{\bf k}^{\prime})$ . Consequently, the lattice self-energy $\Sigma_{\mathrm{DCA^{+}}}({\bf k},{i\omega_{n}})$ is determined by solving the deconvolution problem. We want to point out that the convergence of both the interpolation and the deconvolution depends on the size of the real-space cluster. More precisely, if the assumption already the DCA algorithm is based on is violated, i.e. the self-energy of the lattice is longer-ranged, the lattice mapping problem, $\Sigma_{c}({\bf K},{i\omega_{n}})\rightarrow\Sigma_{\mathrm{DCA^{+}}}({\bf k},{i\omega_{n}})$ , is ill-defined and convergence poor [9]. This sensitivity to the behavior of the real-space self-energy defines the current limits of the DCA ${}^{+}\>$ algorithm.

From $\Sigma_{\mathrm{DCA}^{+}}({\bf k},{i\omega_{n}})$ we again calculate the coarse-grained Green’s function $\bar{G}({\bf K},{i\omega_{n}})$ and this time also a coarse-grained self-energy $\bar{\Sigma}({\bf K},{i\omega_{n}})$ , which should be close to the cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ ,

[TABLE]

The coarse-grained self-energy $\bar{\Sigma}({\bf K},{i\omega_{n}})$ enters the calculation of the bare Green’s function of the effective cluster problem in DCA ${}^{+}\>$ ,

[TABLE]

As in the DCA algorithm, the QMC cluster solver closes the self-consistency loop by generating a new cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ . The full DCA ${}^{+}\>$ algorithm is illustrated in Fig. 3, in which we also point out the differences to the standard DCA loop.

Distinguishing the coarse-grained self-energy $\bar{\Sigma}({\bf K},{i\omega_{n}})$ from the cluster self-energy $\Sigma_{c}({\bf K},{i\omega_{n}})$ , the DCA ${}^{+}\>$ algorithm requires a modified convergence criterion. Self-consistency is reached when the QMC solver produces a cluster self-energy $\Sigma_{c}^{(n)}$ that agrees, within the statistical Monte Carlo error $\varepsilon_{\mathrm{QMC}}$ , with the coarse-grained self-energy $\bar{\Sigma}^{(n-1)}$ used to set up the cluster problem,

[TABLE]

The benefits of the DCA+ algorithm in terms of full momentum resolution and independence of the results on the cluster shape are illustrated in Fig. 4, where the momentum dependence of the DCA and DCA+ lattice self-energy is shown for the two different 16-site clusters of Fig. 1. The DCA self-energy is characterized by jump discontinuities between the coarse-graining patches and shows a strong cluster shape dependence in the region between ( $\pi,0$ ) and ( $0,\pi$ ). In contrast, the DCA+ algorithm provides a self-energy with smooth momentum dependence and results for the two clusters that agree well. In addition to curing the cluster shape dependence, the DCA ${}^{+}\>$ algorithm further reduces the negative sign problem of the underlying QMC method. This was attributed to the removal of the apparent discontinuities in the DCA self-energy, which had caused artificial long-range correlations [9].

2.3 Continuous-time auxiliary-field quantum Monte Carlo

Developed by Gull et al. [18] specifically for quantum impurity problems, the continuous-time auxiliary field QMC algorithm (CT-AUX) is the state of the art method to solve the effective cluster problem in DCA and DCA+ calculations. CT-AUX can treat Hamiltonians containing all types of density-density interactions, including spatially non-local Coulomb repulsion and multiple orbital models, but the algorithm is most useful for large cluster calculations of the single-band Hubbard model. Instead of describing the algorithm at full length, we want to focus on the relevant parts for a DCA++ user and refer the interested reader to the original paper as well as the detailed review on continuous-time QMC methods given in Ref. [7].

We consider a generalized Hubbard model given by the Hamiltonian

[TABLE]

The indices $i$ and $j$ run over sites and orbitals to allow for the most general type of electron hopping. $\mu$ and $\nu$ represent combined indices of spin, orbital and site, and the sum $\sum_{\mu>\nu}$ runs over all pairs of correlated spin-orbitals ( $U_{\mu\nu}\neq 0$ ), whose number we denote by $N_{\mathrm{corr}}$ . Note that for the traditional single-band Hubbard model with only on-site Coulomb repulsion, the number of correlated orbitals is just the number of sites, $N_{\mathrm{corr}}=N_{c}$ . Based on this decomposition of the Hamiltonian $H$ into a non-interacting part $H_{0}$ and the interaction term $H_{\mathrm{int}}$ , we can formulate the partition function $Z=\mathrm{Tr}[e^{-\beta H}]$ in the corresponding interaction representation,

[TABLE]

The time-ordered exponential is now expanded with the parameter $K$ , introduced in Eq. (22), controlling the expansion order,

[TABLE]

The interaction terms $\left(1-\beta H_{\mathrm{int}}(\tau)/K\right)$ are decoupled with the help of auxiliary spins $s_{i}$ [19],

[TABLE]

As a result, the partition function $Z$ becomes a high-dimensional integral,

[TABLE]

which can be evaluated by a Markov Chain Monte Carlo procedure.

Following the Metropolis algorithm, we sample the configuration space $\mathcal{C}$ , in which an element $c\in\mathcal{C}$ consists of a set of $k$ vertices, $\{v=(\mu_{i},\mu_{j},\tau,s)\}_{k}$ . Each vertex possesses an auxiliary spin $s=\pm 1$ and represents an interaction between two correlated spin-orbitals, $\mu_{i}=\left\{\sigma_{i},\alpha_{i},{\bf r}_{i}\right\}$ and $\mu_{j}=\left\{\sigma_{j},\alpha_{j},{\bf r}_{j}\right\}$ , acting at a time $\tau\in[0,\beta)$ . A configuration $c\in\mathcal{C}$ in the Markov chain is updated by inserting a new vertex, removing an existing vertex, or flipping an auxiliary spin. These single spin updates correspond to rank-1 operations on the $N_{\sigma}$ -matrices defined in Eq. (27). The dimension of the $N_{\uparrow}(N_{\downarrow})$ -matrix is given by the number of $\uparrow(\downarrow)$ -electron spins contained in the current configuration. If interactions only occur between spins of opposite type, both $N_{\uparrow}$ and $N_{\downarrow}$ are $k\times k$ matrices. Fig. 5 shows a typical configuration of auxiliary spins on the imaginary time axis and the three types of single spin updates.

While we sample configurations on the imaginary time axis, measuring the Green’s function is most accurately performed in discrete Matsubara frequency space. Time measurements would lead to binning errors as the vertices can occupy arbitrary positions in the imaginary time interval. From a vertex configuration $\{v\}_{k}$ and its respective $N_{\sigma}$ -matrices we first compute a function $M$ , whose Fourier coefficients are then accumulated. The Green’s function is linearly related to $M$ ,

[TABLE]

The self-energy $\Sigma({\bf k},{i\omega_{n}})$ is finally gained from the measured Green’s function $\langle G({\bf k},{i\omega_{n}})\rangle$ and the bare propagator of the cluster impurity problem $\mathcal{G}_{0}({\bf k},{i\omega_{n}})$ using the Dyson equation,

[TABLE]

In addition to the single-particle Green’s function $G({\bf k},{i\omega_{n}})$ , also the two-particle Green’s function can be obtained within the CT-AUX algorithm. The measurement of the latter involves significantly higher computational cost, though.

The CT-AUX algorithm scales as the cube of the linear size of the $N_{\sigma}$ -matrices, i.e. $\mathcal{O}(\langle k\rangle^{3})$ . The average expansion order $\langle k\rangle$ , in turn, grows linearly with the expansion parameter $K$ , the inverse temperature $\beta$ , the strength the interactions $U_{\mu\nu}$ , and the number of interaction terms $N_{\mathrm{corr}}$ . In simulations of the Hubbard model away from half-filling, CT-AUX exhibits a negative sign problem. The signal-to-noise ratio decreases exponentially when the system size or the strength of the Coulomb interaction is increased, or the temperature lowered. Therefore, the behavior of both the average expansion order and the sign problem makes CT-AUX calculations of large systems, with strong coupling and at low temperatures computationally most demanding.

Submatrix updates

An important improvement in terms of computational efficiency of the CT-AUX algorithm has been achieved with the introduction of submatrix updates [20]. Combining $k_{s}$ successive single spin updates into a rank- $k_{s}$ operation, submatrix updates substantially improve data locality and thus benefit from the cache hierarchy of modern processors.

The formalism of submatrix updates assumes that individual spins are changed only once within a submatrix update and spins that do not satisfy this condition require special treatment at the cost of a small performance loss. One can either update these spins with the help of Bennett’s algorithm [21] or abort the current submatrix update when a spin is selected twice and propose the same spin again at the beginning of the next update. For systems with large expansion orders, for which the probability of touching the same spin twice within a submatrix update is small, ignoring proposals of already chosen spins is a valid option, too. This approximation avoids the performance loss, but violates the detailed balance condition.

Non-equidistant FFT measurements

A remarkable speed-up in the measurement of single- as well as the two-particle functions has been accomplished by developing optimized non-equispaced fast Fourier transform (NFFT) methods [11]. In the context of CT-AUX, the methods are used to measure and accumulate the $M$ function, but their applicability extends to other continuous-time QMC methods as well.

Since vertex configurations are not sampled on an equispaced grid but at random time locations, we cannot directly apply the fast Fourier transform (FFT) algorithm on the raw samples. To circumvent this problem, the NFFT algorithm [22] makes use of the convolution theorem,

[TABLE]

Given samples $f_{i}$ at random times $\tau_{i}\in[0,\beta)$ , we first transform the data-set onto an equispaced imaginary time grid by convoluting it with a localized, $\beta$ -periodic kernel $\varphi$ . We are then able to use the FFT algorithm to compute the Fourier transform of the convoluted data. According to Eq. (30), we obtain the Fourier components $f_{\omega}$ of the original data after normalizing with the Fourier coefficients of the kernel $\varphi_{\omega}$ , which finalizes the NFFT algorithm,

[TABLE]

In Eq. (31), $N_{\omega}$ denotes the desired number of positive Matsubara frequencies and $m$ is an oversampling factor.

The linearity of the Fourier transform allows us to delay it, meaning we can accumulate the equispaced data points $\bar{f}_{l}$ and perform a single FFT on the accumulated data at the end of the Monte Carlo integration. The algorithm can be further accelerated by approximating the kernel $\varphi$ with a low order polynomial at the cost of a controllable interpolation error.

While the delayed-NFFT scheme cannot be applied to the measurement of the two-particle Green’s function as Fourier transformation and accumulation are not interchangeable there, a specifically tuned 2D NFFT algorithm can still lead to a considerable speed-up compared to the standard 2D non-equidistant discrete Fourier transform [11].

2.4 Continuous-time hybridization expansion quantum Monte Carlo

With the continuous-time hybridization expansion QMC algorithm (CT-HYB) [23, 24], DCA++ provides an alternative continuous-time impurity solver. While CT-AUX is the favorable algorithm for large cluster DCA and DCA ${}^{+}\>$ simulations, CT-HYB is the method of choice for single-site multi-orbital impurity models in the strong coupling regime ( $U/t\approx 10-\infty$ ). The CT-HYB algorithm is capable of treating general interactions, but provides a particularly efficient segment formulation for density-density interactions. We will restrict ourselves to this formulation and again refer readers to the original papers and Ref. [7] for a complete discussion.

We consider a quantum impurity model, which describes a small system, the impurity, embedded in an infinite non-interacting bath,

[TABLE]

According to its name, the CT-HYB algorithm is based on an expansion of the partition function $Z$ in the hybridization $H_{\mathrm{hyb}}$ , which couples the impurity to the bath. In analogy to the procedure in CT-AUX, the partition function is transformed into a high-dimensional integral, which is sampled by means of a Markov chain Monte Carlo simulation. In the segment picture the configuration space $\mathcal{C}$ is represented by multiple imaginary time lines. Each line spans the interval $[0,\beta)$ and corresponds to a flavor (spin, orbital, site) of the impurity. A configuration $c\in\mathcal{C}$ consists of a set of segments on these time lines during which the corresponding spin-orbitals are occupied. For the Markov chain to be ergodic, two types of updates are needed: segment insertions and removals. Additional updates in form of insertion and removal of “anti-segments” and segment shifts, however, can improve the sampling efficiency. Segment swaps, which exchange the time lines of two flavors, may be necessary in certain cases. A typical configuration and examples of segment updates are illustrated in Figs. 7 and 7, respectively.

The CT-HYB algorithm accumulates the Matsubara coefficients of the Green’s function $G({i\omega_{n}})$ . In addition, it has proven useful to measure $G({i\omega_{n}})\Sigma({i\omega_{n}})$ . Although it is not necessary to accumulate this product, since the self-energy can be obtained from the Dyson equation (29), the approach leads to much smaller statistical errors [25].

For a diagonal hybridization term $H_{\mathrm{hyb}}$ , the CT-HYB algorithm in the segment formulation is sign problem free and scales as $\mathcal{O}(N\beta^{3})$ , where $N$ is the number of spin-degenerate orbitals.

2.5 Calculation of two-particle correlation functions

A possible approach to study phase transitions, for example to the superconducting state, is the extension of the DCA and DCA*+* frameworks to the two-particle level. In this section we want to briefly outline this formalism for the case of the particle-particle channel and refer to the article by Jarrel et al. [6] for a more detailed derivation in the DCA framework, to the work of Staar, Maier, and Schulthess [26] for the extension to DCA+, and to the didactic summary by Maier [27], which we follow here.

DCA

We start from the Bethe-Salpeter equation (BSE), which in the particle-particle channel (abbreviated with “pp”) is given by

[TABLE]

It relates the two-particle Green’s function $G_{2}$ to the irreducible vertex function of this channel $\Gamma^{pp}$ , and can be regarded as the two-particle level analog of the Dyson equation. In Eq. (33) we introduced a combined index for momentum and Matsubara frequency, $k\equiv({\bf k},{i\omega_{n}})$ , with the corresponding summation, $\frac{T}{N}\sum_{k}\equiv\frac{T}{V_{\mathrm{BZ}}}\int d{\bf k}\sum_{i\omega_{n}}$ . The variable $q\equiv({\bf q},{i\nu_{m}})$ represents the transferred momentum and transferred (bosonic) frequency. Applying the underlying assumption of the DCA at the two-particle level, we expand the irreducible vertex function of the cluster onto the coarse-graining patches to obtain a piecewise constant approximation of the lattice irreducible vertex function,

[TABLE]

where $\Gamma^{pp}(k,k^{\prime},q)\equiv\Gamma^{pp}(k,q-k,q-k^{\prime},k^{\prime})$ was used to simplify the notation. The cluster vertex function $\Gamma_{c}^{pp}(K,K^{\prime},q)$ , in turn, satisfies the cluster version of the Bethe-Salpeter equation (33),

[TABLE]

with the non-interacting two-particle Green’s function $G^{0,pp}_{2,c}(K,K^{\prime},q)\equiv G_{c,\uparrow}(K)G_{c,\downarrow}(q-K)\delta_{K,K^{\prime}}$ . Solving for the irreducible vertex function yields (matrix notation in $K,K^{\prime}$ for fixed $q$ ),

[TABLE]

The single- and two-particle Green’s functions are obtained from the cluster solver.

Possible phase transitions can be discovered from the eigenvalues and left and right eigenvectors of the Bethe-Salpeter kernel,

[TABLE]

where we dropped the dependence of $\lambda_{\alpha}$ , $g_{\alpha}(k)$ and $\tilde{g}_{\alpha}(k)$ on the transferred momentum and frequency $q$ , as well as on the considered channel. This can be seen when the two-particle Green’s function is expressed in terms of these eigenvalues and eigenvectors,

[TABLE]

The eigenvalues $\lambda_{\alpha}$ indicate an instability of the system: Eq. (39) diverges when the leading eigenvalue approaches 1. The eigenvectors $g_{\alpha}(k)$ and $\tilde{g}_{\alpha}(k)$ describe the momentum and frequency structure of the instability.

By using the approximation of the lattice irreducible vertex function by step functions, Eq. (34), and restricting the momentum resolution of the right eigenvectors $g_{\alpha}({\bf k},{i\omega_{n}})$ to the cluster momenta ${\bf K}$ , we can transform Eq. (38) into a coarse-grained eigenvalue equation,

[TABLE]

with

[TABLE]

Eq. (40) is easily solvable by standard diagonalization techniques.

DCA ${}^{+}\>$

To restore the momentum resolution of the lattice, we need to find a lattice irreducible vertex function $\Gamma^{pp}(k,k^{\prime},q)$ that is continuous in ${\bf k}$ and ${\bf k}^{\prime}$ . In analogy to the derivation of DCA ${}^{+}\>$ in Section 2.2, we start from a coarse-graining equation for the vertex function,

[TABLE]

which requires the coarse-grained average of the lattice function to match the corresponding cluster quantity. Consistent with the computation of the continuous lattice self-energy $\Sigma_{\mathrm{DCA}^{+}}({\bf k},{i\omega_{n}})$ , we obtain $\Gamma^{pp}_{\mathrm{DCA}^{+}}(k,k^{\prime},q)$ from an interpolation of the cluster vertex function,

[TABLE]

followed by a deconvolution,

[TABLE]

From $\Gamma_{\mathrm{DCA}^{+}}^{pp}(k,k^{\prime},q)$ we can compute the Bethe-Salpeter eigenvectors and eigenvalues with full momentum resolution (see Eq. (38)).

In some cases it has proven useful to project the Bethe-Salpeter kernel onto a set of crystal harmonics corresponding to localized lattice vectors. This folding not only leads to much smaller matrices to diagonalize, but more importantly reduces noise by projecting onto the physically relevant subspace.

3 The DCA++ code

3.1 Main components

The main purpose of this section is to give an overview of the key parts and main capabilities of the DCA++ code. A visual summary of the code base is shown in Fig. 8.

We have chosen C++ as the primary implementation language for its generic programming model, which, besides easy extensibility, allows us to hide architectural details. The C++14 standard, which DCA++ adopts, provides many modern language features for effective and safe programming, while still having a wide compiler support. On top of C++, we use CUDA to utilize GPU resources in performance critical parts of the code.

Cluster mapping

Distinguishing DCA and DCA+ from finite-size methods, the coarse-graining procedure is a central element of the code. In addition to the traditional coarse-graining, which partitions the first Brillouin zone according to the Brillouin zones of the superlattice, the interlaced coarse-graining developed by Staar et al. [17] is implemented with adjustable number of periods, that is variable degree of interleaving. Besides the coarse-graining routines, the cluster mapping module contains methods to update the chemical potential and execute the cluster-exclusion step. In LDA+DMFT calculations the double-counting problem is addressed by the correspondent submodule with a constant correction term [28], specified as an input parameter.

Lattice mapping

The framework of DCA ${}^{+}\>$ requires a non-trivial lattice mapping to generate a lattice self-energy and irreducible vertex function with full momentum resolution. After interpolation of the cluster functions onto the fine lattice grid, the lattice mapping becomes a deconvolution problem. In the case of the self-energy, we employ the iterative Richardson-Lucy deconvolution algorithm as explained in [9]. For the deconvolution of the irreducible vertex function we have chosen a simpler approach. After Fourier transformation of the convolution equation (44) to real space, the deconvolution problem presents itself as a simple division.

Cluster solver

Solving the effective cluster problem and computing a new estimate of the cluster self-energy is the most time consuming part in a DCA or DCA ${}^{+}\>$ calculation. For this reason, lots of effort has been put in the development of a cutting edge implementation of the CT-AUX algorithm including submatrix updates and efficient NFFT measurements. With the aim of featuring optimal methods for the entire spectrum of quantum impurity models, we provide a complementary cluster solver to CT-AUX in the form of the segment formulation of CT-HYB restricted to single-site problems (SS-CT-HYB). Besides the two continuous-time QMC methods, DCA++ implements two other cluster solvers, each bound to a special purpose. A high temperature series expansion (HTS) solver may help improve convergence of the lattice mapping by computing a fourth order perturbation expansion of the self-energy with respect to the on-site Coulomb interaction $U$ . To validate QMC methods and new models, an ED solver can be employed to compute exact results on small lattices in finite-size calculations (no mean-field).

Models

DCA++ is able to simulate standard, extended and frustrated Hubbard models on various types of lattices. 2D square and triangular lattices, and a bilayer square lattice corresponding to a two-orbital model are currently supported, but new lattices and models can easily be added upon need. In particular, the code’s generic structure based on C++ templates permits to easily treat 1D and 3D lattices, too. In the framework of LDA+DMFT real materials such as CuO2 and NiO can be studied in ab initio calculations.

Input/output

Running simulations requires reading input files and writing out results. The input and output (I/O) module supports two data formats, HDF5 [29] and JSON [30]. While JSON’s human-readability makes it very suitable for specifying input parameters, we recommend using the binary HDF5 format to write output. HDF5 features high speed in reading and writing data, and requires less storage size compared to other data formats.

Parallelization

The computationally most expensive part of a DCA or DCA ${}^{+}\>$ calculation is the QMC step. Fortunately, Monte Carlo simulations are embarrassingly parallel, which we exploit on distributed multi-core machines with a two level (MPI + threading) parallelization scheme, illustrated in Fig. 9.

Building DCA++ with MPI enabled allows to run the DCA or DCA ${}^{+}\>$ calculation with multiple MPI processes. Each MPI process performs an independent Monte Carlo integration and results are accumulated at the end in few collective reductions. In addition, a threaded Monte Carlo wrapper allows each MPI process to spawn multiple Monte Carlo walker and accumulator threads. Each Monte Carlo walker processes its own Markov chain and, as soon as it has executed the required number of updates per measurement, sends its current configuration over to a queue of idle Monte Carlo accumulators. A free accumulator then takes measurements on the sampled configuration. The total number of measurements during the Monte Carlo step of a single DCA or DCA ${}^{+}\>$ iteration is given by

[TABLE]

GPU

In addition to conventional multi-core architectures, DCA++ ports to the emerging hybrid CPU-GPU systems. Comprising large matrix-matrix multiplications, the CT-AUX cluster solver with submatrix updates particularly benefits from the computing power of the GPU. If DCA++ is built with GPU support enabled, all Monte Carlo walkers of a process are moved to its dedicated GPU. This leads to a considerable speed-up compared to the multi-core implementation [13].

Applications

main_dca is the central application of DCA++. An input parameter allows to choose at runtime whether the DCA or DCA ${}^{+}\>$ algorithm should be executed. One can turn off the mean-field and do a finite-size QMC calculation, too. Depending on input parameters, either a fixed number of iterations are carried out, or the DCA or DCA ${}^{+}\>$ algorithm is iterated until the desired accuracy is reached. The main output of the application is the cluster’s single-particle Green’s function and self-energy. The lattice self-energy can be written out, too, when DCA ${}^{+}\>$ is used. If required, the application additionally measures and stores the two-particle Green’s function of the cluster.

The output of main_dca is post processed by the application main_analysis. This application computes the eigenvectors and eigenvalues of the Bethe-Salpeter equation, which provide information about nature and extent of instabilities in the system and possible phase transitions. Like main_dca, main_analysis can be run in DCA or DCA ${}^{+}\>$ mode. While in DCA mode the momentum resolution of the BSE eigenvectors is restricted to the cluster momenta, the BSE solver produces eigenvectors with continuous momentum dependence when DCA ${}^{+}\>$ is used.

3.2 Building

Prerequisites

Building DCA++ requires a set of tools and libraries:

•

a C++14 compiler (tested: clang $\geq$ 3.5, gcc $\geq$ 4.9),

•

the CMake build system [31] version 3.3 or higher,

•

the HDF5 library with the C++ interface (tested: 1.8.13 and 1.10),

•

FFTW3 or an FFT library with the FFTW3 interface (e.g. Intel MKL),

•

BLAS and LAPACK libraries, and

•

MPI, if requested.

To enable the GPU support, one needs to provide

•

the CUDA Toolkit (tested: 7.0, 7.5 and 8.0), and

•

the MAGMA library (tested: 2.0.0 with CUDA 7.0 and 2.2.0 with CUDA 7.5 and 8.0).

Further optional requirements are

•

the EasyBuild framework [32],

•

Python, NumPy, SciPy, Matplotlib, and the h5py package to use the provided python scripts.

Building steps

The DCA++ project is maintained in a public GitHub repository at https://github.com/CompFUSE/DCA. Cloning the repository allows to obtain the latest version of the master branch. Alternatively, one can download a specific version from the release page. The recommended way of configuring and building DCA++ across all platforms is to use the CMake build system. The building procedure starts with creating a clean build directory. The cmake command followed by the path to the source and appropriate options generates the build files. Finally, by evoking make one can compile the applications and tests, if enabled. The building steps are summarized in Listing 1. More details on the building procedure can be found on the correspondent Wiki page in the GitHub repository, while a complete list of all relevant CMake options with descriptions and their default value is provided in A.

3.3 Running

Input files

There are two types of input files. All applications read simulation parameters, e.g. the DCA cluster, the temperature or output filenames, from a JSON-formatted input file, in which parameters are thematically grouped and sub-grouped in JSON objects. B provides descriptions of all input parameters including their type and default value.

In addition to providing simulation parameters, data needs to be transferred between individual runs. A main_dca run uses the output of a previous calculation to initialize the self-energy with a non-zero value (see below) and the application main_analysis requires a main_dca output file containing the cluster single-particle and two-particle Green’s functions, as well as the self-energy. While the file format for data can either be HDF5 (recommended) or JSON, it has to match between consecutive runs.

Parallelization

On modern systems with multi-core nodes we usually run up to two MPI tasks per node and one Monte Carlo walker or accumulator thread per core. The optimal distribution of Monte Carlo walker and accumulator threads, on the other hand, is strongly problem dependent. Note that memory restrictions on the GPU might limit the number of walker threads a process can spawn. The environment variable OMP_NUM_THREADS can be used to enable threading in BLAS, LAPACK and FFTW routines. One should be aware, however, that some of the calls to these routines are already within parallel regions.

Cooldown

A full DCA or DCA ${}^{+}\>$ calculation is usually performed by slowly cooling down the system. In a series of main_dca runs the temperature is gradually lowered to study its dependence on fluctuations leading to phase transitions. Sometimes this procedure is also necessary to guarantee convergence of the DCA and DCA+ loop at low temperatures.

The calculation is started at a high temperature ( $T\approx 1.0t$ ), at which a vanishing self-energy is a good initial guess, and the system then slowly cooled down. Each following main_dca run is initialized with the result of the previous temperature. While for the first temperature up to eight iterations are required until convergence is reached, three to six iterations are usually enough for the subsequent runs, given the temperature steps are sufficiently small. Lowering the temperature is usually done with decreasing step size.

In the beginning only single-particle (sp) functions, i.e. the cluster Green’s function and self-energy, are computed. When we reach the relevant temperature regime in which the system exhibits an instability, we additionally measure the required cluster two-particle Green’s function in a separate two-particle (tp) run, which only performs a single DCA or DCA ${}^{+}\>$ iteration. Subsequent to the main_dca runs, we separately analyze the output of all tp-runs with main_analysis. Listing 2 provides a schedule of jobs for a typical cooldown procedure.

We provide a python script that generates directories, input files and batch scripts for a cooldown. By default, this script and the provided template input files are configured for DCA or DCA ${}^{+}\>$ calculation of the 2D single-band Hubbard model with on-site Coulomb interaction $U$ and fixed electron filling $d$ , but it can easily be adjusted for other problems. Usage of the script is explained in the Wiki.

4 Sustainable and scalable software development

The DCA++ project follows several well-proven methods in software engineering to develop and maintain a high-performance research code in a multidisciplinary team. This section is devoted to sharing these methods and the tools that allow to implement them. In addition to the dedicated article [16], we want to call the interested reader’s attention to the “How to” documents of the IDEAS Scientific Software Productivity Project [33], which collect best practices in scientific software development, and which we partially follow here.

4.1 Software testing and test-driven development

Systematic testing is an essential element in software development to produce correct code. In addition to detecting errors, software tests provide an important form of documentation by describing correct usage and defining requirements and limits with edge cases.

Software tests can be divided according to their level into unit, integration, and system-level tests. While unit tests operate at the smallest code scale testing single functions and classes, integration tests check the interaction between these units, and system-level tests validate and verify the behavior of the full code. Orthogonal to this classification, we can categorize tests by their aspect into verification, validation, no-change, and performance tests [33]. Verification tests determine whether the code satisfies specific requirements or conditions. They answer the question “Are we building the product right?” [34]. Verification testing at the unit level is the tool kit of test-driven development (see below). Validation tests, in contrast, are motivated by the question “Are we building the right product?” [34]. These tests evaluate whether the software satisfies the intended purpose. No-change tests define the current behavior of the software and are indispensable for refactoring legacy code. Passing this type of software tests requires reproducing predefined results close to machine precision. Performance tests are used to compare alternative algorithms, to monitor the code’s performance, or to benchmark computing systems.

All of our software tests are automated and part of the regression test suite. The regression test suite is divided with respect to required computational resources into fast tests and extensive tests, while performance tests are kept separated. Fast tests include unit tests and lightweight integration tests. Since they build and run fast, they can be executed frequently in the code development to detect problems early. Extensive tests contain computationally more expensive integration tests and system-level tests, and are usually executed before a piece of code gets back into the main trunk (see workflow below). They provide additional evidence that the code’s behavior has not been unintentionally altered.

Test-driven development (TDD) [35] is a software development technique in which automated tests steer the design process. Based on the two rules of only writing new code if there is a failing test and removing duplication, TDD consistently follows a three-step cycle:

Create an automated test that fails. 2. 2.

Make the test pass quickly by any means necessary. 3. 3.

Refactor to eliminate duplication introduced by getting the test to pass.

As any functional code is protected by at least one unit test, TDD prevents regression and delivers high confidence in the correctness of the software. Writing tests first has the additional benefit that it encourages the programmer to think about the interface of the code first. This way, TDD fosters designs that are high in cohesion and low in coupling. The rhythm of TDD allows to work in small, safe steps and to focus on one task at a time. The tests provide immediate feedback.

In the DCA++ project, we follow TDD whenever the existing code allows it and refactor code to make it more testable. Where unit tests cannot be put in place, system-level no-change tests act as a global safety net. Moreover, we require defects to be reported in the form of minimal, failing tests. Our goal is a growing test suite and an increasing percentage of test driven code.

We implement tests using the Google C++ Testing Framework (Google Test) [36]. Google Test builds on the xUnit architecture and provides benefits both at writing and executing tests. It supports various fatal and non-fatal assertions, including death tests, and allows tests to share resources that are expensive to set up and not modified by the tests. Moreover, tests can be value- or type-parameterized to cover a wide range of cases without duplicating test code.

4.2 Software project management and collaborative development

Tools

Successfully maintaining any software project begins with keeping the code base under version control. By storing old versions of the software and monitoring changes to the code base, a version control system crucially contributes to reproducibility of simulation results and safe code development. Version control systems that support branching allow the isolated implementation of features in dedicated branches, while the master branch, which contains the production code, is always operational. We use the popular Git version control system [37], in which each developer clones their own local repository. Besides being distributed, Git offers fast and easy branching capabilities, which makes it well suited for collaborative software development. Our Git repository is hosted on the GitHub platform, which, in addition to the repository, provides a set of tools for managing a software project enabling us to implement the workflow explained below.

While GitHub helps us to maintain the DCA++ source, we employ the EasyBuild framework to manage builds of DCA++ on production systems. All dependencies of the DCA++ code, such as external libraries and compilers, but also the CMake build system, are specified in a build recipe. Bundled in the so-called “easyconfig” file, the dependencies can be installed with a one-line command into a single module111More commonly, EasyBuild is used to not only build the dependencies but the full application or library.. This module can be loaded and used by other collaborators to easily build the DCA++ code. We currently provide easyconfig files for CSCS’ Piz Daint (hybrid Cray XC50 / XC40) and ORNL’s Titan (Cray XK7) supercomputers.

To automate the quality control of the code, we run a Jenkins server [38] hooked up to DCA++’s GitHub repository. Triggered by changes to the repository, Jenkins builds all applications and tests, and runs the test suite. A test report allows to immediately detect any degradation and to document the status of the code for later reference.

Workflow

It is our goal to continuously and frequently deliver updates to the DCA++ code without compromising its quality. To achieve this, we organize the code development by tracking each task, bug fix or other enhancement as an issue. GitHub’s issue tracker offers to categorize issues with labels and to assign them to a team member. By associating issues with a milestone, the process to a specific goal can be monitored. The implementation of a single issue should not exceed few weeks and more complex features, which we call epics, should be broken down into smaller, independent tasks. To visualize the status of all issues and to prioritize work, we use a Kanban board, a popular tool in agile development. The board is organized in columns with each column representing a stage in the lifetime of an issue, i.e. Backlog, In progress, Review and Done. A typical Kanban board for a software project is shown in Fig. 10.

By developing each issue in a dedicated branch, we guarantee that the master branch is kept stable at all times. As part of the quality control, the finished work has to be submitted in form of a pull request. An experienced developer then reviews the code and requests changes if necessary. If the reviewer approves the new code and Jenkins confirms that all tests pass, it will be merged into the master branch and become production code. Fig. 11 illustrates this issue branch workflow.

For a successful realization of the workflow it is important to have early discussions about design decisions and implementation details. This assures that the development of a new issue starts immediately in the right direction. Collaborative development in large teams once more requires keeping individual issues smalls. As part of continuous integration, we demand developers to rapidly check in their changes via pull requests to master. This allows them to get feedback frequently and prevents a scenario, called merge hell, in which it involves considerable effort to integrate a far diverged branch back into the main line. Short issue cycles also mean that new features and bug fixes are quickly available to team members and users. This can, in particular, avoid duplication of code work.

5 Examples

The following two standard use cases study superconductivity in the 2D Hubbard model. The lattice model is described by the Hamiltonian

[TABLE]

where $c_{i\sigma}^{\dagger}$ ( $c_{i\sigma}$ ) creates (annihilates) an electron with spin $\sigma$ on lattice site $i$ , $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ is the corresponding number operator and $U$ denotes the on-site Coulomb repulsion. In the form above, where $\langle i,j\rangle$ indicates that the first sum only runs over pairs of nearest neighbor sites, the model is restricted to nearest neighbor hopping with isotropic amplitude $t$ . The simulations are performed on a square lattice.

In the first example, we determine the superconducting transition temperature $T_{c}$ for a small $2\times 2$ cluster employing the DCA method. The second example provides a large-scale DCA+ calculation in which we converge $T_{c}$ with respect to the cluster size $N_{c}$ . In both cases we use the CT-AUX QMC algorithm to solve the effective cluster problem.

5.1 Computing the superconducting transition temperature $T_{c}$ with DCA

As explained in section 2.5, we can analyze the transition to the superconducting phase by computing the leading eigenvalue of the Bethe-Salpeter equation in the particle-particle channel. The superconducting state of the 2D Hubbard model is characterized by $d$ -wave symmetry [26, 39]. In Fig. 12 we plot the leading $d$ -wave eigenvalue $\lambda_{d}$ as a function of temperature for strong coupling $U=8t$ , electron filling $\langle n\rangle=0.95$ and $N_{c}=4$ . As the temperature is lowered, the eigenvalue increases and the system approaches the phase instability. After fitting the data points with a function of the form

[TABLE]

we can solve for the transition temperature $T_{c}$ , defined as the temperature for which $\lambda_{d}(T)$ crosses one, i.e. $\lambda_{d}(T_{c})=1$ .

We performed the computation on Piz Daint’s multi-core partition, a Cray XC40 system equipped with two 18-core Intel Xeon processors (2.1 GHz) per node. Using a single compute node, results for 1 000 000 Monte Carlo measurements per iteration were obtained in 3.5 hours. The time-to-solution (TTS) can be greatly reduced by exploiting the code’s parallel performance and increasing the number of nodes. To quantify this, we measured the strong scaling efficiency $E_{s}(n)$ as a function of the number of nodes $n$ . It is defined as

[TABLE]

Results are shown in the left plot of Fig. 13. Even though the computational demand of the problem is small, DCA++ allows to efficiently use up to 200 compute nodes, for which the TTS reduces to less than two minutes.

We also examined the weak scaling behavior, where the problem size (number of Monte Carlo measurements) per node is kept constant. Weak scaling efficiency $E_{w}(n)$ as a function of the number of nodes $n$ is computed as

[TABLE]

The right plot of Fig. 13 presents the results for 1 000 000 Monte Carlo measurements per iteration and node. While overall, we observe weak scaling efficiency close to the ideal case of 100%, the 10 node result actually shows a deviation to “better than ideal”. This can be explained by the fact that we only kept the problem size with respect to the number of Monte Carlo measurements constant per node. Other parts of the code, in particular the coarse-graining step, benefited from a speed-up when the number of nodes was increased.

5.2 Cluster size scaling of $T_{c}$ with DCA+

The reduced fermionic sign problem in the DCA+ algorithm makes clusters sizes accessible that have been out of reach for standard DCA or finite-size QMC calculations. This allows to enter a regime of asymptotic convergence, in which the transition temperature $T_{c}$ becomes independent of the cluster size. Here, we consider the weak-coupling case of $U=4t$ at 10% doping. Fig. 15 shows the temperature dependence of the leading $d$ -wave eigenvalue $\lambda_{d}(T)$ for cluster sizes between 24 and 56 and Fig. 15 plots the resulting values for $T_{c}$ versus cluster size $N_{c}$ . We observe convergence of the transition temperature $T_{c}$ for clusters larger than 32 sites.

Error bars in Figs. 15 and 15 represent, with the exception of $N_{c}=36$ , the statistical error of the QMC algorithm. They tend to grow with the size of the cluster. This originates from the behavior of the fermionic sign problem, which becomes more severe with increasing cluster size leading to a larger statistical error.

For the case of $N_{c}=36$ , we used four different cluster shapes to illustrate the error associated with the choice of cluster for given cluster size. While the cluster shape dependence is substantially reduced by the introduction of DCA ${}^{+}\>$ [9], the corresponding spread in results still exceeds the statistical uncertainty of the Monte Carlo algorithm. Ref. [26] extends this analysis to more cluster sizes.

Large expansion orders resulting in large matrix operations make this problem computationally demanding, but also well suited for Piz Daint’s Cray XC50 hybrid compute nodes. They provide, in addition to a 12-core Intel Xeon processor (2.6 GHz), an NVIDIA Tesla P100 GPU accelerator. This leadership computing system was used to examine the large-scale parallel performance of the DCA++ code. We measured the TTS of the QMC integration for the computing intensive case of $N_{c}=56$ and $T=0.02t$ . Due to the problem size, we chose a minimum of 500 computes nodes to limit the run time. For illustration, the full cooldown for $N_{c}=56$ takes about 5.2 hours on 500 nodes. We performed 1 000 000 Monte Carlo measurements for each node count to investigate the strong scaling behavior and used 2000 Monte Carlo measurements per node in the weak scaling analysis. Fig. 16 presents strong and weak scaling results for up to 5000 compute nodes. The observed strong scaling efficiency implies that we achieve a considerable speed-up over 500 nodes when using nearly the entire machine. As in the multi-core case, weak scaling efficiency excels reaching almost 100%.

We want to conclude this example by looking at the speed-up that a Cray XC50 hybrid node provides over a Cray XC40 multi-core node. As in the scaling analysis, we compare the TTS of the QMC integration for the 56-site cluster at $T=0.02t$ . The total number of Monte Carlo measurements per node was fixed to 2000. We can see in Fig. 17 that the speed-up in a single-particle run is close to threefold, while the two-particle run is accelerated by a factor of 2.5. The smaller speed-up when the computationally expensive cluster two-particle Green’s function is accumulated accounts to the fact that the hybrid implementation of the CT-AUX solver performs the measurements on the host CPU.

6 Final remarks

Future development

DCA++ is an active software project with constantly extended functionality. The next release will add full multi-orbital support and a state of the art implementation of the continuous-time interaction expansion QMC algorithm (CT-INT) [40] to allow for simulations of multi-orbital models of systems such as the iron-based superconductors.

Since Monte Carlo algorithms are based on random sampling, they are non-deterministic by nature. While implementations of these algorithms can still be tested at the unit level, it is often difficult to design conventional integration or system-level tests due to the random behavior. For this reason, we have developed a framework based on hypothesis testing, which is able to validate stochastic programs. This framework is currently being integrated into DCA++.

The workload balance in the CT-AUX solver on hybrid CPU-GPU systems, i.e. running the Monte Carlo walker on the GPU and accumulating measurements on the CPU, has proven to be efficient on current systems. However, to respond to the emerging GPU-dominated node configurations, the next release of DCA++ will provide the capability to move the CT-AUX accumulator to the GPU, too. In the long term, the growing diversity of HPC architectures will make it necessary to revisit the on-node parallelization strategy altogether. We plan to follow a task based parallelization model based on the HPX-3 system [41] to ensure portability and efficiency for the future.

Getting started

The DCA++ code is licensed under the BSD 3-Clause License. Hence, it is open source and a version of the code can easily be obtained from the public GitHub repository at https://github.com/CompFUSE/DCA. We advise prospective users to visit the Wiki section of the repository, which, in particular, provides a small tutorial. Questions about building or running DCA++ can be posted to the repository’s built-in issue tracking system.

Contributing back

There are various ways to contribute to the DCA++ project. Reports about bugs and other problems can be submitted via the issue tracking system. The issue tracker is also the preferred communication channel for any other request concerning the code, including the suggestion of new features. We welcome code contributions, too. They are best submitted by creating a pull request in the GitHub repository.

Citation guideline

If the DCA++ code contributes to your research, we kindly ask to acknowledge this in any resulting scientific publication by citing this paper.

Acknowledgments

The authors would like to acknowledge the contributions of Giovanni Balduzzi, Peter W. Doak, Mi Jiang, Andrei Plamada and Bart Ydens to the DCA++ project. We would like to thank John Biddiscombe for his contributions to refactoring the code base, Guilherme Peretti-Pezzi for his assistance in setting up continuous integration and his support with EasyBuild, and Philipp Werner for providing the CT-HYB QMC solver. U.R.H. acknowledges the support by the National Center of Competence in Research (NCCR) MARVEL, funded by the Swiss National Science Foundation. The work of T.A.M. was supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, Division of Materials Sciences and Engineering. This research used resources of the Oak Ridge Leadership Computing Facility (OLCF) awarded by the INCITE program, and of the Swiss National Supercomputing Center (CSCS). OLCF is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725.

Appendix A CMake options

The DCA++ build can be configured with several CMake options. These options can either be appended to the cmake command with the -D flag:

⬇

$ cmake ../dca_source -D<variable1>=<value1> -D<variable2>=<value2> ...

or specified using the ccmake interface, which is launched in the existing build directory (after cmake was called) with:

⬇

$ ccmake .

Table 1 lists some general and all DCA++ specific CMake options.

Appendix B Input parameters

This appendix section provides short descriptions of each input parameter using the format:

•

"parameter-name": type (default value)

Description.

In addition, Listing 3 shows a sample input file as used for the DCA+ calculation of Section 5.2.

Output parameters

Group "output":

•

"directory": string ("./")

Directory to write the output to.

•

"output-format": string ("HDF5")

File format of the output files. Options are: "HDF5", "JSON".

•

"filename-dca": string ("dca.hdf5")

Filename for the output of the application main_dca.

•

"filename-analysis": string ("analysis.hdf5")

Filename for the output of the application main_analysis.

•

"filename-ed": string ("ed.hdf5")

Filename for the ED output in the application cluster_solver_check.

•

"filename-qmc": string ("qmc.hdf5")

Filename for the QMC output in the application cluster_solver_check.

•

"filename-profiling": string ("profiling.json")

Filename for the profiling output. The file format is always JSON.

•

"dump-lattice-self-energy": boolean (false)

Write out the lattice self-energy $\Sigma({\bf k},{i\omega_{n}})$ in the application main_dca.

•

"dump-cluster-Greens-functions": boolean (false)

Write out various cluster Green’s functions in the application main_dca.

•

"dump-Gamma-lattice": boolean (false)

Write out the lattice irreducible vertex function $\Gamma(k,k^{\prime},q)$ in the application main_analysis.

•

"dump-chi-0-lattice": boolean (false)

Write out $\chi_{0}(k)$ in the application main_analysis.

Physics parameters

Group "physics":

•

"beta": double (1.)

Inverse temperature.

•

"density": double (1.)

Target electron density.

•

"chemical-potential": double (0.)

Initial value of the chemical potential. This value is overwritten, if an output file with an initial self-energy is read (see section DCA parameters).

•

"adjust-chemical-potential": boolean (true)

Adjust the chemical potential to obtain the specified density.

Model parameters

$\diamond$

Single-band Hubbard model

Used if square lattice or triangular lattice are selected.

Group "single-band-Hubbard-model":

$\bullet$

"t": double (0.)

Nearest neighbor hopping parameter.

$\bullet$

"t-prime": double (0.)

Next nearest neighbor hopping parameter.

$\bullet$

"U": double (0.)

On-site Coulomb repulsion.

$\bullet$

"V": double (0.)

Nearest neighbor coulomb repulsion for opposite spins.

$\bullet$

"V-prime": double (0.)

Nearest neighbor coulomb repulsion for same spins.

$\diamond$

Bilayer Hubbard model

Used if bilayer lattice is selected.

Group "bilayer-Hubbard-model":

$\bullet$

"t": double (0.)

Nearest neighbor hopping parameter.

$\bullet$

"t-prime": double (0.)

Next nearest neighbor hopping parameter.

$\bullet$

"t-perp": double (0.)

Hopping parameter between layers.

$\bullet$

"U": double (0.)

On-site Coulomb repulsion.

$\bullet$

"V": double (0.)

Coulomb repulsion between layers for opposite spins.

$\bullet$

"V-prime": double (0.)

Coulomb repulsion between layers for same spins.

$\diamond$

Material model

Group "material-model":

$\bullet$

"t_ij-filename": string ("t_ij.txt")

Name of the CSV file containing the hopping matrix $t_{ij}$ .

$\bullet$

"U_ij-filename": string ("U_ij.txt")

Name of the CSV file containing the Coulomb repulsion matrix $U_{ij}$ .

DCA parameters

Group "DCA":

•

"initial-self-energy": string ("zero")

Either the name of the file (including path) with the initial self-energy (usually the main_dca output file of the previous temperature) or "zero" indicating that the initial self-energy should be zero.

•

"iterations": integer (1)

Number of DCA/DCA+ iterations.

•

"accuracy": double (0.)

Stop the DCA/DCA+ loop if this accuracy has been reached.

•

"self-energy-mixing-factor": double (1.)

Parameter $\alpha$ of Eq. (12).

•

"interacting-orbitals": array of integers ([0])

Indices of orbitals that are treated interacting orbitals. This parameter must be consistent with the model that is used.

•

"do-finite-size-QMC": boolean (false)

Do a finite-size QMC calculation (no mean-field).

•

Subgroup "coarse-graining":

$\bullet$

"k-mesh-recursion": integer (0)

Number of recursion steps in the creation of the momentum space mesh. See Ref. [17] for details.

$\bullet$

"periods": integer (0)

Number of periods of the interlaced coarse-graining patches, i.e. degree of their interleaving. Restricted by k-mesh-recursion. See Ref. [17] for details.

$\bullet$

"quadrature-rule": integer (1)

Determines the quadrature rule used in the coarse-graining.

quadrature-rule $<$ 0: use a flat mesh.

quadrature-rule $>=$ 0: use the Grundmann-Moeller rule [42] of index = quadrature-rule.

A rule of index $s$ is exact for polynomials up to order $2s+1$ .

$\bullet$

"threads": integer (1)

Number of threads used in the coarse-graining.

$\bullet$

"tail-frequencies": integer (0)

Number of tail frequencies used for updating the chemical potential.

•

Subgroup "DCA+":

$\bullet$

"do-DCA+": boolean (false)

Use the DCA+ algorithm.

$\bullet$

"deconvolution-iterations": integer (16)

Maximum number of iterations in the deconvolution step.

$\bullet$

"deconvolution-tolerance": double (1.e-3)

Termination criteria for the deconvolution step.

$\bullet$

"HTS-approximation": boolean (false)

Use a high-temperature series approximation in the lattice mapping.

$\bullet$

"HTS-threads": integer (1)

Number of threads used in the HTS-solver.

Domains parameters

Group "domains":

•

Subgroup "real-space-grids:

$\bullet$

"cluster": array of arrays of integers (lattice basis, e.g. in 2D: [[1, 0], [0, 1]])

Real space DCA cluster. Given as coordinates with respect to the lattice basis.

To choose a cluster of a given size, consult the exhaustive list of 2D and 3D clusters in

tools/cluster_definitions.txt.

$\bullet$

"sp-host": array of arrays of integers (lattice basis, e.g. in 2D: [[1, 0], [0, 1]])

Real space host grid for single-particle functions, also called (sp-)lattice. Given as coordinates with respect to the lattice basis.

$\bullet$

"tp-host": array of arrays of integers (lattice basis, e.g. in 2D: [[1, 0], [0, 1]])

Real space host grid for two-particle functions, also called tp-lattice. Only used in DCA+.

Given as coordinates with respect to the lattice basis.

•

Subgroup "imaginary-time":

$\bullet$

"sp-time-intervals": integer (128)

Discretization of the imaginary time axis.

$\bullet$

"time-intervals-for-time-measurements": integer (1)

Discretization of the imaginary time axis for additional time measurements.

•

Subgroup "imaginary-frequency":

$\bullet$

"sp-fermionic-frequencies": integer (256)

Number of fermionic Matsubara frequencies for single-particle functions.

$\bullet$

"HTS-bosonic-frequencies": integer (0)

Number of bosonic Matsubara frequencies in the HTS-solver.

$\bullet$

"four-point-fermionic-frequencies": integer (1)

Number of fermionic Matsubara frequencies for four-point functions.

•

Subgroup "real-frequency":

$\bullet$

"min": double (-10.)

Minimum real frequency.

$\bullet$

"max": double (10.)

Maximum real frequency.

$\bullet$

"frequencies": integer (3)

Number of real frequencies.

$\bullet$

"imaginary-damping": double (0.01)

Small positive shift $\eta$ , used to shift the real frequencies $\omega$ away from the real axis, $\omega\rightarrow\omega+i\eta$ .

Monte Carlo integration parameters

Group "Monte-Carlo-integration":

•

"seed": integer or string (985456376)

Seed for the random number generator(s) used in the Monte Carlo integration. Instead of an integer, the string "random" can be passed to generate a random seed.

•

"warm-up-sweeps": integer (20)

Number of warm-up sweeps.

•

"sweeps-per-measurement": double (1.)

Number of sweeps per measurement.

•

"measurements-per-process-and-accumulator": integer (100)

Number of independent measurements each accumulator of each process performs.

•

Subgroup "threaded-solver":

$\bullet$

"walkers": integer (1)

Number of Monte Carlo walkers.

$\bullet$

"accumulators": integer (1)

Number of Monte Carlo accumulators.

Monte Carlo solver parameters

$\diamond$

CT-AUX

Used if CT-AUX is selected as the cluster solver.

Group "CT-AUX":

$\bullet$

"expansion-parameter-K": double (1.)

The perturbation order in the CT-AUX algorithm increases linearly with the expansion parameter $K$ . While $K$ is only subject to the restriction of being positive, values of the order of one have proven to be a good choice [18].

$\bullet$

"initial-configuration-size": integer (10)

The CT-AUX solver is initialized with initial-configuration-size random interacting vertices.

$\bullet$

"initial-matrix-size": integer (128)

Initial size of the CT-AUX matrices.

$\bullet$

"max-submatrix-size": integer (128)

Maximum number of single spin updates per submatrix update.

$\bullet$

"neglect-Bennett-updates": boolean (false)

Neglect the two types of Bennett updates:

Removal of an interacting spin that has already been proposed for removal. 2. 2.

Removal of a non-interacting spin that has been proposed for insertion.

Turning on this option leads to a larger number of delayed spins ( $\approx$ max-submatrix-size) at the cost of a systematic error due to the violation of detailed balance.

$\bullet$

"additional-time-measurements": boolean (false)

Do additional time measurements.

$\diamond$

SS-CT-HYB

Used if SS-CT-HYB is selected as the cluster solver.

Group "SS-CT-HYB":

$\bullet$

"self-energy-tail-cutoff": integer (0)

Cut-off parameter for the imaginary frequency tail of the self-energy.

$\bullet$

"steps-per-sweep": double (0.5)

$\bullet$

"shifts-per-sweep": double (0.5)

The fraction of insert/removal steps is given by steps-per-sweep/(steps-per-sweep + shifts-per-sweep), while the fraction of segments shifts is shifts-per-sweep/(steps-per-sweep + shifts-per-sweep).

Four-point parameters

Group "four-point":

•

"type": string ("NONE")

Four-point type, options are:

–

"NONE"

–

"PARTICLE_PARTICLE_UP_DOWN"

–

"PARTICLE_HOLE_CHARGE"

–

"PARTICLE_HOLE_MAGNETIC"

–

"PARTICLE_HOLE_TRANSVERSE"

•

"momentum-transfer": array of doubles (null vector with the dimension of the lattice, e.g. for a 2D lattice: [0., 0.])

Transferred momentum $\mathbf{q}$ . Must be an element of the reciprocal lattice of the DCA cluster ( $\left\lVert\mathbf{q}-\mathbf{K}_{\mathrm{DCA}}\right\rVert_{2}<10^{-3}$ ).

•

"frequency-transfer": integer (0)

Transferred frequency $\omega_{n}$ . Given as the index $n$ of the bosonic Matsubara frequency $\omega_{n}=2n\pi/\beta$ .

Analysis parameters

Group "analysis":

•

"symmetrize-Gamma": boolean (true)

Symmetrize the cluster and lattice irreducible vertex functions according to the symmetry group and diagrammatic symmetries.

•

"Gamma-deconvolution-cut-off": double (0.5)

Cut-off parameter for the deconvolution of cluster irreducible vertex function. The deconvolution is done by Fourier transforming to real space, dividing by the coarse-graining patch functions in real space, $\phi(r)$ , and then transforming back to reciprocal space. The real space patch functions are only inverted down to the value of this parameter,

[TABLE]

•

"project-onto-chrystal-harmonics": boolean (false)

Project the product $\Gamma(k,k^{\prime})\chi_{0}(k)$ onto the crystal harmonic functions and diagonalize in this subspace.

•

"projection-cut-off-radius": double (1.5)

For the projection use crystal-harmonic functions of lattice vectors $\mathbf{r}$ with $\left\lVert\mathbf{r}\right\rVert_{2}<$ projection-cut-off-radius.

ED-solver parameters

Only required if the exact diagonalization solver is used, e.g. in the application cluster_solver_check.

Group "ED":

•

"eigenvalue-cut-off": double (1.e-6)

Only keep energy eigenvalues $E$ for which $e^{-\beta E}>$ eigenvalue-cut-off.

Double-counting parameters

Only required for LDA+DMFT.

Group "double-counting":

•

"method": string ("none")

The double-counting method to use, options are:

–

"none": no double-counting correction

–

"constant-correction-without-U-correction"

–

"constant-correction-with-U-correction"

•

"correction": double (0.)

The value of the double-counting correction.

Bibliography42

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] P. W. Anderson, The Resonating Valence Bond State in La 2Cu O 4 and Superconductivity, Science 235 (4793) (1987) 1196–1198.
2[2] F. C. Zhang, T. M. Rice, Effective Hamiltonian for the superconducting Cu oxides, Physical Review B 37 (7) (1988) 3759–3761.
3[3] A. Georges, G. Kotliar, W. Krauth, M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Reviews of Modern Physics 68 (1) (1996) 13–125.
4[4] T. Maier, M. Jarrell, T. Pruschke, M. H. Hettler, Quantum cluster theories, Reviews of Modern Physics 77 (3) (2005) 1027–1080.
5[5] M. H. Hettler, M. Mukherjee, M. Jarrell, H. R. Krishnamurthy, Dynamical cluster approximation: Nonlocal dynamics of correlated electron systems, Physical Review B 61 (19) (2000) 12739–12756.
6[6] M. Jarrell, T. Maier, C. Huscroft, S. Moukouri, Quantum Monte Carlo algorithm for nonlocal corrections to the dynamical mean-field approximation, Physical Review B 64 (19) (2001) 195130.
7[7] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, P. Werner, Continuous-time Monte Carlo methods for quantum impurity models, Reviews of Modern Physics 83 (2) (2011) 349–404.
8[8] S. Moukouri, M. Jarrell, Absence of a Slater Transition in the Two-Dimensional Hubbard Model, Physical Review Letters 87 (16) (2001) 167010.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

DCA++: A software framework to solve correlated electron problems with modern quantum cluster methods

Abstract

keywords:

1 Introduction

2 Formalism and methods

2.1 The dynamical cluster approximation

2.2 The DCA+ {}^{+}\>+algorithm

2.3 Continuous-time auxiliary-field quantum Monte Carlo

Submatrix updates

Non-equidistant FFT measurements

2.4 Continuous-time hybridization expansion quantum Monte Carlo

2.5 Calculation of two-particle correlation functions

DCA

DCA+ {}^{+}\>+

3 The DCA++ code

3.1 Main components

Cluster mapping

Lattice mapping

Cluster solver

Models

Input/output

Parallelization

GPU

Applications

3.2 Building

Prerequisites

Building steps

3.3 Running

Input files

Parallelization

Cooldown

4 Sustainable and scalable software development

4.1 Software testing and test-driven development

4.2 Software project management and collaborative development

Tools

Workflow

5 Examples

5.1 Computing the superconducting transition temperature TcT_{c}Tc​ with DCA

5.2 Cluster size scaling of TcT_{c}Tc​ with DCA+

6 Final remarks

Future development

Getting started

Contributing back

Citation guideline

Acknowledgments

Appendix A CMake options

Appendix B Input parameters

Output parameters

Physics parameters

Model parameters

DCA parameters

Domains parameters

Monte Carlo integration parameters

Monte Carlo solver parameters

Four-point parameters

Analysis parameters

ED-solver parameters

Double-counting parameters

2.2 The DCA ${}^{+}\>$ algorithm

DCA ${}^{+}\>$

5.1 Computing the superconducting transition temperature $T_{c}$ with DCA

5.2 Cluster size scaling of $T_{c}$ with DCA+