Skip to content

User Guide

A step-by-step guide to using XDiag

Installation

Julia

Enter the package mode using ] in the Julia REPL and type: 

add XDiag

That's it!

C++

For using the C++ version, we need to compile the library first. To do so, we first download the code using git, 

cd /path/to/where/xdiag/should/be
git clone https://github.com/awietek/xdiag.git

and then compile the library using CMake,

cd xdiag
cmake -S . -B build
cmake --build build
cmake --install build

The resulting library is now installed at /path/to/where/xdiag/should/be/install. There are various options when compiling, including optimizations which can be used. For more details on the compilation process, we refer to the Compilation guide.

First steps

Writing code

Let us set up our first program using the xdiag library. 

#include <xdiag/all.hpp>

using namespace xdiag;

int main() try {
  say_hello();
} catch (Error e) {
  error_trace(e);
}

using XDiag
say_hello()

The function say_hello() prints out a welcome message, which also contains information which exact XDiag version is used. In Julia this is all there is to it.

What is maybe a bit unfamiliar is the try / catch block in C++. XDiag implements a traceback mechanism for runtime errors, which is activated by the error_trace function. While not stricly necessary here, it is a good practice to make use of this.

Compilation

In C++, now that the application program is written, we next need to set up the compilation instructions using CMake. To do so we create a second file called CMakeLists.txt in the same directory.

cmake_minimum_required(VERSION 3.19)

project(hello_world)

find_package(xdiag REQUIRED HINTS "/path/to/where/xdiag/should/be/install")
add_executable(main main.cpp)
target_link_libraries(main PRIVATE xdiag::xdiag)

You should replace /path/to/where/xdiag/should/be/install with the appropriate directory where your XDiag library is installed after compilation. We then compile the application code,

cmake -S . -B build
cmake --build build
and finally run our first XDiag application.

./build/main

Hilbert spaces

We are now ready to run our first actual calculation using XDiag. Our immediate goal will be to determine the ground state energy of a \(S=1/2\) Heisenberg model on a 1D chain lattice with periodic boundary conditions,

\[ H = J \sum_{\langle i, j\rangle} \mathbf{S}_i \cdot \mathbf{S}_j,\]

where \(\mathbf{S}_i = (S^x_i, S^y_i, S^z_i)\) denotes the vector of spin matrices at a given site \(i\). The notation \(\langle i, j\rangle\) refers to summatation over neighboring sites \(i\) and \(j\).

The first thing to define before any computation, is the Hilbert space our model will be defined on. For the present example, we use the Hilbert space class Spinhalf. Further possible Hilbert spaces include tJ and Electron, see Blocks. We consider a chain lattice with \(N=8\) sites and create a \(S=1/2\) Hilbert space:

using namespace xdiag;
int N = 8;
auto hspace = Spinhalf(N);

N = 8
hspace = Spinhalf(8)

We would like to know which spin configurations, the Hilbert space is made up of. To do so, we can iterate over the Hilbert space and print out the spin configurations. 

for (auto spins : hspace) {
  Log("{}", to_string(spins));
}

for spins in hspace
    println(to_string(spins))
end

This produces an output similar to the following:

↓↓↓↓↓↓↓↓
↓↓↓↓↓↓↓↑
↓↓↓↓↓↓↑↓
↓↓↓↓↓↓↑↑
↓↓↓↓↓↑↓↓
↓↓↓↓↓↑↓↑
↓↓↓↓↓↑↑↓
...

Here we already see several things at work. XDiag features a convenient way to write logs in C++, with the Log class. The first argument to Log() is a format string. In C++ we use the fmt library, to be able to write structured format and format our output. The second argument turns our spins into a string. spins is of type ProductState, whose configuration on each site can be individually addressed.

Further, we notice that all \(2^N\) spin configurations are included in this Hilbert space. However, the Heisenberg model conserves the total \(S^z = \sum_i S^z_i\), and thus we could limit ourselves to a block of the Hilbert space, which only contains configurations of a certain magnetization: 

int nup = 4;
auto block = Spinhalf(N, nup);
for (auto spins : block) {
  Log("{}", to_string(spins));
}

nup = 4;
block = Spinhalf(N, nup);
for spins in block
    println(to_string(spins))
end

This produces an output similar to: 

↓↓↓↓↑↑↑↑
↓↓↓↑↓↑↑↑
↓↓↓↑↑↓↑↑
↓↓↓↑↑↑↓↑
↓↓↓↑↑↑↑↓
↓↓↑↓↓↑↑↑
↓↓↑↓↑↓↑↑
...

We see that the block now only contains configurations, where the total number of spins pointing up is 4. We can now quickly check the dimension of the Hilbert spaces, and confirm that the dimension of the block is reduced from \(2^8=256\) to \(\begin{pmatrix} 8 \\ 4 \end{pmatrix} = 70\),

XDIAG_SHOW(size(hspace));
XDIAG_SHOW(size(block));

@show size(hspace)
@show size(block)

which should print: 

size(hspace):
256
size(block):
70

Here, we introduced another functionality of XDiag in C++, the XDIAG_SHOW macro which can be used for quick debug printing of XDiag objects.

Operators

Next, we define our Hamiltonian. We do so by using an OpSum object. These encode sums of operators of the form

\[ \mathcal{O} = \sum_i c_i \mathcal{O}_i. \]

Here, \(\mathcal{O}_i\) denote single operators described as Op objects and \(c_i\) are the coupling constants. These can either be a real or complex number, or a string which later needs to be replaced. THe Hamiltonian of the spin \(S=1/2\) Heisenberg chain is created like this:

auto ops = OpSum();
for (int i=0; i<N; ++i) {
  ops += "J" * Op("SdotS", {i, (i+1) % N});
}
ops["J"] = 1.0;

ops = OpSum()
for i in 1:N
    ops += "J" * Op("SdotS", [i, mod1(i+1, N)])
end
ops["J"] = 1.0

We first create an empty OpSum and then add additional terms to it. The first part of the product denotes the coupling constant, here given as a string. Alternatively, one could have directly used real / complex numbers here. The second part of the product is a single Op object. It is created with two inputs:

  1. The type, here SdotS denoting an operator of the form \(\mathbf{S}_i\cdot\mathbf{S}_j\). XDiag features a wide variety of operator types, see Operator types.
  2. An array defining which site the operator lives on. Notice, that in julia we start counting the sites from 1, while in C++ we start counting the sites from 0.

Ground states

Now our Hamiltonian is defined, we can directly compute the ground state and its energy using the function eig0:

auto [e0, psi0] = eig0(ops, block);
Log("e0: {:.12f}", e0);

e0, psi0 = eig0(ops, block)
println("e0:", e0)

Here, e0 is a double precision real number and psi0 is a State object. The eig0 uses an iterative algorithms to compute the ground state by applying the Hamiltonian in a matrix-free manner, and is thus very useful for large system sizes. It is based on a Lanczos algorithm, which can also be called directly using eigs_lanczos which offers more control on the behavior of the algorithm.

Alternatively, we can also perform a full diagonalization by computing the full Hamiltonian matrix using the matrix function.

arma::mat H = matrix(ops, block); 
arma::vec evals;
arma::mat evecs;
arma::eig_sym(evals, evecs, H);
Log("e0: {:.12f}, e1: {:.12f}", evals[0], evals[1]);

H = matrix(ops, block); 
(evals, evecs) = eigen(Symmetric(H))
println("e0:", evals[1], "e1:", evals[2]);

Notice, that we are using the Armadillo library with the arma namespace in C++. This library serves as the linear algebra backend of XDiag and can be used to perform further calculations. In Julia, the eigen and Symmetric functions are part of the LinearAlgebra standard library.

Measurements

Finally, we might be interested in measuring some observables of the ground state. We compute the static spin correlation \(S^z_0 S^z_j\):

for (int i=1; i<N; ++i) {
  auto op = Op("SzSz", {0, i});
  double corr = inner(op, psi0);
  Log("<Sz_0 Sz_{}> = {:.12f}", i, corr);
}

for i in 2:N
    op = Op("SzSz", [1, i])
    corr = inner(op, psi0)
    println("<Sz_0 Sz_$i> = ", corr)
end

Here, we are using the inner function to compute an expectation value of the form \(\langle \psi_0 |\mathcal{O}| \psi_0 \rangle\) with \(\mathcal{O} = S^z_0 S^z_j\).

Input / Output

Julia features a variety of packages facilitating input and output of data. For C++, XDiag provides convenient functionality for TOML and HDF5 files.

Reading from TOML

While defining a Hamiltonian or other operators in code as above can be done, it is often preferable to define operators in a file and read it in. In XDiag, we use the TOML language to define input. The Hamiltonian of the \(N=8\) site Heisenberg chain we created above can be written in a TOML file as a list like this:

Interactions = [
  ['J', 'SdotS', 0, 1],
  ['J', 'SdotS', 1, 2],
  ['J', 'SdotS', 2, 3],
  ['J', 'SdotS', 3, 4],
  ['J', 'SdotS', 4, 5],
  ['J', 'SdotS', 5, 6],
  ['J', 'SdotS', 6, 7],
  ['J', 'SdotS', 7, 0]
]

The first entry in every list element is the coupling constant J, the second enty is the type SdotS, and the following two entries are the sites of the operator. To read in such a Hamiltonian from a toml file we can use the FileToml object together with the read_opsum Function,

auto fl = FileToml(XDIAG_DIRECTORY "/examples/user_guide/spinhalf_chain.toml");
auto ops = read_opsum(fl, "Interactions");

fl = FileToml("spinhalf_chain.toml")
ops = read_opsum(fl, "Interactions")

XDiag also features the functions read_permutation_group and read_representation to conventiently read PermutationGroup and Representation objects used to describe symmetries from file, see further below.

Writing results to hdf5

After a simulation, we want to store results to file. A standard scientific data format is the hdf5 format. Julia supports input and output to hdf5 with the HDF5.jl package.

For C++ we provide a convenient way of writing results to hdf5 files. In general all numerical data, including scalar real/complex numbers as well as armadillo vectors and matrices can be easily written to file using the FileH5 object.

auto fl = FileH5(XDIAG_DIRECTORY "/examples/user_guide/output.h5", "w!");
fl["e0"] = e0;
fl["evals"] = evals;
fl["evecs"] = evecs;

The second argument "w!" specifies the access mode to the file. "w!" is the forced write mode, where preexisting files are overwritten.

Symmetries

Permutations

Key functionality of XDiag comprises symmetry-adapted calculations, in particular with permutation symmtries which can be translation symmetries on a lattice. The Heisenberg chain with periodic boundary conditions is invariant under the translation operator

\[ T: i \rightarrow i+1 \]

To represent such symmetries, we use the Permutation object. For example, on the \(8\)-site chain lattice, the symmetry which maps the sites,

\[ \{0,1,2,3,4,5,6,7\} \rightarrow \{1,2,3,4,5,6,7,0\} \]

can be represented via: 

auto T = Permutation({1, 2, 3, 4, 5, 6, 7, 0});

T = Permutation([2, 3, 4, 5, 6, 7, 8, 1])

Notice, that also here we start counting from 1 in Julia, and from 0 in C++. Permutation objects can be multiplied, inverted, and raised to a power. We can use them to define a PermutationGroup object in the following way:

auto group = PermutationGroup({pow(T, 0), pow(T, 1), pow(T, 2), pow(T, 3),
                               pow(T, 4), pow(T, 5), pow(T, 6), pow(T, 7)});

group = PermutationGroup([T^0, T^1, T^2, T^3, T^4, T^5, T^6, T^7])

Representations

One-dimensional representations of a PermutationGroup are described by the Representation object. As one-dimensional irreducible representations (irreps) are given by their characters \(\chi(g)\) which is simply a list of (complex) numbers for every symmetry element \(g\). A representation can be created by handing the group and the list of characters: 

auto irrep_k_0 = Representation(group, arma::vec{1.0, 1.0, 1.0, 1.0,
                         1.0, 1.0, 1.0, 1.0});
auto irrep_k_pi = Representation(group, arma::vec{1.0, -1.0, 1.0, -1.0,
                          1.0, -1.0, 1.0, -1.0});

irrep_k_0 = Representation(group, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
irrep_k_pi = Representation(group, [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0])

Upon creation of a representation, XDiag checks whether the group axioms as well as the homomorphism property of the characters is fulfilled, i.e.,

\[ f * g = h \Rightarrow \chi(f) \cdot \chi(g) = \chi(h).\]

Representations can then be used to create symmetry adapted blocks, which can then be diagonalized independently:

auto block_k_0 = Spinhalf(N, nup, irrep_k_0);
auto block_k_pi = Spinhalf(N, nup, irrep_k_pi);
double e0_k_0 = eigval0(ops, block_k_0);
double e0_k_pi = eigval0(ops, block_k_pi);
Log("e0: k=0: {:.12f}, k=pi: {:.12f}", e0_k_0, e0_k_pi);

block_k_0 = Spinhalf(N, nup, irrep_k_0)
block_k_pi = Spinhalf(N, nup, irrep_k_pi)
e0_k_0 = eigval0(ops, block_k_0)
e0_k_pi = eigval0(ops, block_k_pi)
println("e0: k=0: $e0_k_0, k=pi: $e0_k_pi")

Also, permutation groups and representations can be defined in a TOML file and read in. Please refer to the read_permutation_group and read_representation documentation on the specific format required.

Using symmetries, thus, can give enhanced physical insights to the system, but also reduces the computational costs significantly. A common analysis tool for understanding quantum many-body systems is the Tower of states analysis. For an introduction, see e.g.

Studying Continuous Symmetry Breaking using Energy Level Spectroscopy
Alexander Wietek, Michael Schuler, Andreas M. Läuchli
arXiv:1704.08622 [cond-mat.str-el]
DOI: 10.48550/arXiv.1704.08622

Further features

Some further features of XDiag include: