Exact Diagonalization library for quantum electron models

TL;DR

EDLib is a C++ library that efficiently performs exact diagonalization for quantum electron models, enabling calculation of Green's functions and spin correlations for complex systems.

Contribution

It introduces a flexible C++ template library with multiple Hamiltonian storage options for solving quantum electron models using exact diagonalization.

Findings

Supports single-band Hubbard and multi-orbital Anderson models

Provides computation of Green's functions and spin correlations

Offers three Hamiltonian storage types for efficiency

Abstract

We present an exact diagonalization C++ template library (EDLib) for solving quantum electron models, including single-band finite Hubbard cluster and multi-orbital impurity Anderson model. The observables that can be computed using EDLib are single particle Green's functions and spin-spin correlation functions. This code provides three different types of Hamiltonian matrix storage that can be chosen based on the model.

Figures

Tables

Table 1. Table 1: Wall-clock time used per iteration as a function of the number of MPI processes on various matrix sizes.

Processes	41409225	64128064
36	0.91 sec	2.47 sec
72	0.55 sec	1.41 sec

Table 2. Table 2: Input parameters description

Parameter name	Description
NSITES	Number of sites
NSPINS	Number of spins
INPUT_FILE	HDF5 input file
Storage parameters
storage.MAX_SIZE	Number of eigenvalues to find
storage.MAX_DIM	Number of eigenvalues to find
storage.EIGENVALUES_ONLY	Compute only eigenvalues
spinstorage.ORBITAL_NUMBER	Number of orbitals with interaction
ARPACK parameters
arpack.SECTOR	Read symmetry sectors from file
arpack.NEV	Number of eigenvalues to find
arpack.NCV	Number of convergent values
Lanczos parameters
lanc.NOMEGA	Number of Matsubara frequencies
lanc.NLANC	100, ”Number of Lanczos iterations
lanc.BETA	10.0, ”Inverse temperature
lanc.BOLTZMANN_CUTOFF	Cutoff for Boltsman factor
single impurity Anderson Model
siam.NORBITALS	Number of impurity orbitals

Table 3. Table 3: Descitized bath parameters for orbitals of different symmetries.

Orbital	$V_k$, eV	${ϵ}_k$, eV
$xy$, $x^2-y^2$	0.56434; 0.68392; 0.29519	-2.37325; -0.87328; 2.01265
$xz$, $yz$	0.81892; 0.99136	-3.15496; -1.69066
$3z^2-r^2$	0.77347; 0.79785	-5.59842; -2.95325

Table 4. Table 4: The resulting electronic configuration for exact diagonalization study of Co adatom adsorbed on Pt(111).

$\Delta E$, eV	n↑	n↓	Sector size	Major contribution to g.s.
0.0	12	15	841568	$\downarrow$ $\uparrow \downarrow$ $\downarrow$ $\uparrow \downarrow$ $\downarrow$
0.0	13	14	1618400	$\uparrow$ $\uparrow \downarrow$ $\downarrow$ $\uparrow \downarrow$ $\downarrow$ + $\downarrow$ $\uparrow \downarrow$ $\downarrow$ $\uparrow \downarrow$ $\uparrow$ + $\downarrow$ $\uparrow \downarrow$ $\uparrow$ $\uparrow \downarrow$ $\downarrow$
0.0	14	13	1618400	$\downarrow$ $\uparrow \downarrow$ $\uparrow$ $\uparrow \downarrow$ $\uparrow$ + $\uparrow$ $\uparrow \downarrow$ $\uparrow$ $\uparrow \downarrow$ $\downarrow$ + $\uparrow$ $\uparrow \downarrow$ $\downarrow$ $\uparrow \downarrow$ $\uparrow$
0.0	15	12	841568	$\uparrow$ $\uparrow \downarrow$ $\uparrow$ $\uparrow \downarrow$ $\uparrow$
13.471	9	8	590976100
13.471	8	9	590976100

Equations

$$\mathcal{H}={\mathcal{H}}_{loc}+{\mathcal{H}}_{hop}.$$

$$H=\sum _{p\sigma }{ϵ}_p{c}_{p\sigma }^{+}{c}_{p\sigma }+$$

$$+$$

$${\mathcal{H}}_{hop}={\mathcal{H}}_{\uparrow }\oplus {\mathcal{H}}_{\downarrow }={\mathcal{H}}_{\uparrow }\otimes {\mathcal{I}}_{\downarrow }+{\mathcal{I}}_{\uparrow }\otimes {\mathcal{H}}_{\downarrow },$$

$$H\left(\begin{array}{c}{x}_1\\ x_2\\ \dots \\ x_3\end{array}\right)=H_{loc}\left(\begin{array}{c}{x}_1\\ x_2\\ \dots \\ x_3\end{array}\right)+\left(\begin{array}{c}{H}_{\downarrow }{x}_1^{(\downarrow )}\\ H_{\downarrow }{x}_2^{(\downarrow )}\\ \dots \\ H_{\downarrow }{x}_3^{(\downarrow )}\end{array}\right)+H_{\uparrow }\left(\begin{array}{c}{x}_1^{(\downarrow )}\\ x_2^{(\downarrow )}\\ \dots \\ x_3^{(\downarrow )}\end{array}\right),$$

$${\hat{c}}_j|m⟩$$

$$s$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

*Program Title: * EDLib

Licensing provisions: MIT

Programming language: C++, MPI

External routines: ARPACK-NG, ALPSCore library[1]

*Nature of problem:

*Finite Hubbard and Anderson models play an essential role in the description of strongly correlated many-particle systems. These models consist of a small number of localized orbitals with Coulomb interaction between electrons and (in case of the Anderson model) non-interacting bath energy levels. The finite Hubbard cluster can be used to study molecular magnets, such as $Mn_{12}$, $Fe_{4}$, $Mn_{4}$, and $V_{15}$, which are currently of interest due to their potential for use in novel technologies such as molecular electronics, solar energy harvesting, thermoelectrics, sensing, and other applications [2, 3, 4]. The Anderson model can be used to study impurities adsorbed on surfaces [5] and appears as an impurity model in the Dynamic Mean Field Theory[6].

*Solution method:

*The OpenMP and MPI parallelized versions of the finite temperature Lanczos diagonalization method is used to diagonalize Hamiltonian matrix and to compute observables.

1 Introduction

Further progress in material science is connected with the development of appropriate theoretical concepts and methods to treat realistic modern materials and devices taking their atomic structure, chemical composition, electronic and magnetic properties fully into account. Two of the basic quantum models for systems with strong electron-electron correlations are the Hubbard model[7] and the Anderson impurity model[8], which can be used to simulate lattice problems or an impurity in metal respectively.

At the moment, there are a number of well-developed numerical techniques one can use to solve these quantum electron models. For instance, many interesting and promising results were obtained by using QMC-type methods such as continuous-time quantum Monte Carlo method [9]. Since the main computational task is a sampling of a complex integral, these methods are ideally suited for parallelization. However, there is a fundamental problem of the QMC solvers called the sign problem, which can occur for models with non-diagonal Coulomb interaction matrix, lattice problem away from half-filling or when the simulation temperature is rather low [10].

Alternatively, truncating the infinite Hilbert space by solving a finite lattice problem or by discretizing an infinite bath with a finite set of energy levels allows one to use exact diagonalization techniques to treat the Anderson Hamiltonian. Such a method allows to diagonalize the electronic Hamiltonian for different geometries of lattice cluster or with different forms of the on-site Coulomb matrix [11, 12]. Another advantage of the exact diagonalization method is that it provides a natural way to calculate real-frequency correlation functions such as one- and two-particle Green’s functions at finite temperatures.

In this work, we present the parallel Exact diagonalization library for solving the eigenvalue problem of the Hubbard model or Anderson impurity model on distributed-memory and shared memory computing systems.

2 Exact diagonalization of finite quantum electron models

The Hamiltonian of the many quantum electron model can be expressed as the sum of local (diagonal) term and non-diagonal hopping term as follows:

[TABLE]

For example in case of Hubbard model [7] $\mathcal{H}_{loc}=\sum\limits_{i}U_{i}n_{i\uparrow}n_{i\downarrow}-\sum\limits_{i\sigma}\mu_{i}n_{i\sigma}$ and $\mathcal{H}_{hop}=\sum\limits_{\langle i,j\rangle\sigma}t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma},$ where $U_{i}$ is Coulomb potential on site $i$; $\mu_{i}$ – chemical potential on site $i$; $t_{ij}$ – hopping integral between sites $i$, $j$. $c^{(\dagger)}_{i\sigma}$ – annihilation (creation) operator of electron with spin direction $\sigma$ on $i$-th site. $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ – occupation number, number of electrons on the site.

The first step of exact diagonalization algorithm is to represent a Hamiltonian operator as a matrix. Despite the fact that for most quantum electron models this matrix is very sparse (99% of matrix elements being zeroes) the dimension still grows exponentially $M=2^{2N_{s}}$ in occupation number space $\ket{n_{1\uparrow},\dots,n_{N_{s}\uparrow}|n_{1\downarrow},\dots,n_{N_{s}\downarrow}}$, where $N_{s}$ is the number of electron levels in the studied quantum electron model. The exponential growth of basis size puts serious restriction on lattice size.

Based on the particle and spin conservation properties of the Hubbard model Hamitonian the matrix assumes block-diagonal form and the blocks – so called sectors – of much smaller dimension $M_{n_{\uparrow}n_{\downarrow}}=C_{N_{s}}^{n_{\uparrow}}\cdot C_{N_{s}}^{n_{\downarrow}}$, where $C_{n}^{k}$ is the number of combinations of k from n elements. The sectors can be diagonalized separately as matrices of local Hamiltonians for fixed total occupation numbers $n_{\uparrow}=\sum_{i=1}^{N_{s}}\hat{n}_{i\uparrow}$, $n_{\downarrow}=\sum_{i=1}^{N_{s}}\hat{n}_{i\downarrow}$ (see Fig.(1)). The size difference for $N_{s}=15$ in case of half-filling is 1.6 orders of magnitude, which is considerable, yet the amount of data and matrix sparsity (99% of the elements are zeroes for both matrices) remain high enough to cause difficulties[13]. In the case of the Lanczos diagonalization algorithm the main procedure is the Krylov subspace construction that requires a matrix-vector product operation.

3 Exact diagonalization of single multi-orbital impurity Anderson Model

The multi-orbital impurity Anderson model can be written in the following general form:

[TABLE]

Here $\epsilon_{i}$ and $\epsilon_{p}$ are energies of the impurity and bath states, $d^{+}_{i\sigma}$ and $c^{+}_{p\sigma}$ are the creation operators for impurity and surface electrons, $V_{ip}$ is a hopping between impurity and surface states, $U_{ijkl}$ is the Coulomb matrix element and the impurity orbital index $i$ ($j$, $k$, $l$) runs over the $d-$ states. Depending on the problem we solve the bath can correspond to either an effective Weiss field (DMFT) or, for instance, metallic surface states (adatom on a substrate).

4 Storage formats

4.1 Spin-resolved Hamiltonian storage format

Since the hopping Hamiltonian does not contain hopping between different spins it can be decompose into two parts for each spins as follows:

[TABLE]

where $\mathcal{I}_{\sigma}$ is the identity matrix with the same dimension as $\mathcal{H}_{\sigma}$, and can be stored separately. Since dimension of $\mathcal{H}_{\sigma}$ is much smaller than the original Hamiltonian matrix, the only problem is to store the eigen-vectors since the Hilbert space still grows exponentially. To deal with this issue in this library we implement the distributed storage of the vector as will be described in the next subsection.

4.1.1 MPI parallelization

In this library for solving the eigenvalue problem we use a parallel version of the implicitly restarted Arnoldi algorithm library [14], which requires implementation of the matrix-vector products. In the case of a matrix decoupled into diagonal and two off-diagonal matrices parts (See Eq. 3) this operation can be performed by three separate operations:

[TABLE]

where $x^{(\downarrow)}_{i}$ is an $i-$sub-vector of initial vector $x$ with a dimension of $dim(H_{\downarrow})$. It is clear to see vector can be simply distributed along different processors by integer numbers of $x^{(\downarrow)}_{i}$ sub-vectors. The only operation that needs to perform inter-processor communication is the last term in the right part of the Eq.4. In this case we can overlap communications and computations by using one-sided MPI communications:

⬇

MPI_Win_fence(MPI_MODE_NOPRECEDE, _win);

// Initiate remote data transfer for up-spin term

MPI_Get(...)

// Compute diagonal contribution.

...

// Compute down-spin contribution

...

// Perform final synchronization

MPI_Win_fence(MPI_MODE_NOSUCCEED | MPI_MODE_NOPUT | MPI_MODE_NOSTORE, _win);

// Compute up-spin contribution and off-diagonal interactions contribution

4.1.2 Scaling properties

Table 1 shows how the program scales with the number of MPI processes the two problems with different matrixs size on the University of Michigan high-performance cluster. We see that the computation time behaves like $1/0.8N_{p}$. The principal reason for this behavior is that a larger number of processes leads to more inter-node communication compared to mostly intra-node communication for a small number of processes. One way to improve the scalability is to dedicate a single core for communications. Work in this direction is currently in progress.

4.2 Signs-only Compressed Row storage format

In the new sparse matrix format – ¡¡Signs Only Compressed Row Storage¿¿ (SOCRS) – we attempt to balance the time-efficiency of CRS with the size-efficiency of on-the-fly generation.

According to estimation, the off-diagonal elements of the Hamiltonian matrix stored in CRS format make the largest contribution to memory footprint. The absolute values of inter-site hoppings $t_{ij}$ are readily available as an adjacency matrix, normally read once on solver initialization, whereas sign depends (for Fermi-Dirac statistics) on the number of sites occupied to the left of the changed state:

[TABLE]

which means their on-the-fly calculation would considerably increase SpMV time. In this situation it is sensible to store only signs in sparse format, one sign per bit. The dense diagonal can be stored in a separate vector. Scanning the adjacency matrix makes the number of elements in the row known, thus making row_ptr of CRS redundant.

5 Program description

The EDLib library is designed to solve the exact diagonalization problem for electronic quantum Hamiltonians. For a large matrix cases we use parallelization by means of MPI or OpenMP (depending storage type). The program is written in standard C++11 and distributed as a template library. The program has been checked using GNU, Intel and Clang C++ compilers. The test run has also been checked on the University of Michigan high-performance high-performance computing (HPC) cluster.

5.1 Class diagram

The structure of the EDLib library can be represented by the UML diagram shown at Fig. (3). The main class of the presented library is the $Hamiltonian$ which is parametrized by type of $Model$ and type of $Storage$. The main method for exact diagonalization is $diag$ that mostly delegate the work to the specified $Storage$ class. The Green’s functions ($GreenFunction$ class) are computed by the Lanczos continued fraction method ($Lanczos$ class).

5.2 Description of the input data

Input data to the EDLib library take the form of (1) parameter file, (2) Model specific HDF5 file. Table 2 represents the complete list of the parameters that can be defined in the parameter file. The structure of HDF5 files is model specific and for its generation we provide Python scripts for each represented model.

6 Prerequisites and Installation

To build the EDLib library, any recent C++11 compiler can be used; the libraries have been tested with GCC [15] 4.6 and above, Intel C++ 15 and above and Clang [16] 3.4 and above.

The library depends on the following packages:

The CMake build system [17] of version 2.8.12 and above.

• 2.

The ALPSCore libraries [1] of version 0.54.0 and above.

• 3.

The ARPACK library [14] of opencollab arpack-ng [18] version 3.5.0 and above.

To use optional distributed or shared memory parallel capabilities (the support depends on chosen matrix storage format), an MPI implementation supporting standard 2.1 [19] and above, or OpenMP version 3.0 [20] must be enabled by USE_MPI or USE_OPENMP, respectively.

The installation of the EDLib library follows the standard procedure for any CMake-based package. The first step is to download the EDLib source code. Assuming that all above mentioned prerequisite software is installed, the installation consists of running CMake from a temporary build directory, as outlined in the shell session example with MPI support below:

⬇

git clone https://github.com/Q-Solvers/EDLib.git

mkdir EDLib-build&& cd EDLib-build

cmake \

-DCMAKE_INSTALL_PREFIX=${HOME}/local/EDLib/ \

-DALPSCore_DIR=${HOME}/local/ALPSCore/share/ALPSCore/ \

-DARPACK_DIR=${HOME}/local/arpack-ng/lib/ \

-DUSE_MPI=TRUE\

../EDLib

make

make test

The command at line 1 will download the latest source code from github; at line 5 the destination install directory of the EDLib libraries is set (${HOME}/local/EDLib in this example).

6.1 Citation policy and contributing

EDLib, as an application of ALPSCore, is an open source project and we encourage feedback and contributions from the user community. Issues should be reported exclusively via the GitHub website at https://github.com/Q-solvers/EDLib/issues. For contributions, we recommend to use the pull request system on the GitHub website. Before any major contribution, we recommend coordination with the main EDLib developers. We kindly request that the present paper be cited in any published work using the EDLib library as well as the ALPSCore library on which the presented library is based also be cited[1]. This helps the EDLib and ALPSCore developers to better keep track of projects using the library and provides them guidance for future developments.

7 Examples

To show the ability of the presented library we consider two problems. The first problem is the groundstate calculation of isolated 4x4 Hubbard cluster. And the second one is the ground state electronic configuration of the single $Co$ impurity adsorbed on the $Pt(111)$ surface.

7.1 Finite Hubbard cluster diagonalization

The Hubbard 4-site plaquette is a minimal and generic electronic-structure model of cuprate superconductors suggested in [21] – which demonstrates critical behavior for certain doping. Its properties have been studied in isolation, in the bath and in the Bethe lattice. We study an isolated system of four such plaquettes with periodic boundary conditions. The schematic representation of the cluster is presented in the Fig. (4). The solid and dashed lines correspond to nearest neighbour hopping $t$ and second nearest neighbour hopping $t^{\prime}$ respectively. For the present calculation the following parameters have been chosen: $t=1.0eV$, $t^{\prime}=-0.3eV$, $U=6.0eV$ and $\mu=0.54eV$. The resulting lowest energy is $E_{0}=-22.6421eV$.

7.2 Co adatom on the Pt(111) surface

The electronic and transport properties of the single transition metal adsorbed on the various type of surface play a crucial role in the proper description of giant magnetic anysotropy [22] or Kondo physics [23, 24, 25]. In this simulation we present results for electronic configuration of the ground state of the single cobalt adatom adsorbed on the Pt(111) surface by means of the Anderson impurity model. The model parameters are obtained from first principles calculation. The experimental value of lattice constant for the bulk fcc Pt is 3.92 Å [26, 27]. Since the relaxation for different stackings, $fcc$ and $hcp$, as shown in previous study [28] does not show much difference, we perform simulation for $hcp$ position of cobalt adatom as presented in Fig. 5.

For the Anderson Impurity model we choose 5 orbitals for $d-$states of cobalt adatom, two orbitals in the bath for each $xz$,$yz$ and $3z^{2}-r^{2}$ cobalt orbitals and three orbitals in the bath for each $xy$ and $x^{2}-y^{2}$ cobalt orbitals. Based on the spin symmetries and block-diagonal structure of the Hamiltonian matrix the dimension of the largest block is about $\approx 590\times 10^{6}$. For the present calculations we choose the following parameters: $U=6.6$ eV, $J_{H}=0.9$ eV, $\mu=44.44$ eV, and the bath parameters are presented in the Table 3. The interaction part of the Hamiltonian is expressed by using Slater integral representation of the full rotational invariant Coloumb interaction tensor [29] with $F^{0}=U$, $F^{2}=14J_{H}/(1+0.625)$ and $F^{4}=0.625F^{2}$ [30]. We perform diagonalization of each symmetry sector to find the electronic configuration of the ground state.

The simulation is performed on Edison Cray machine and takes about 1600 core-hours on 10 nodes with maximum memory requirement about 10 Gb per node. The resulting electronic configuration is presented in Table 4, in addition we present the lowest energy for the half-filled states.

8 Summary

We have presented the free software EDLib, an implementation of the Exact diagonalization solver for Anderson Impurity and finite Hubbard models. Further developments (e.g., support for complex Hamiltonians, other measures) are planned for a future release.

Acknowledgment

We thank E. Gull and A. I. Lichtenstein for useful discussions. SI was supported by the Simons collaboration on the many-electron problem and by Act 211 Government of the Russian Federation, contract No. 02.A03.21.0006, M.D. acknowledge support from Deutsche Forschungsgemeinschaft via Project SFB 925. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.