The Model class

The main class for setting up a model is the Model class, which acts as a container for model specific parameters such as the Hamiltonian, temperature, chemical potential, etc. For a simple example, lets consider a simple two level system described by the Hamiltonian

\(H = U_{0}c_{0}^{\dagger}c_{0} + U_{1}c_{1}^{\dagger}c_{1} + V\left(c_{0}^{\dagger}c_{1} + H.c.\right)\),

which is to be studied at T = 300, and zero chemical potential. Before setting up the model we write the Hamiltonian on the canonical form

\(H = \sum_{\mathbf{i}\mathbf{j}}a_{\mathbf{i}\mathbf{j}}c_{\mathbf{i}}^{\dagger}c_{\mathbf{j}}\),

where \(a_{00} = U_0\), \(a_{11} = U_1\), and \(a_{01} = a_{10} = V\). Once on this form, the model is ready to be fed into a Model object, which is achieved as follows

double U_0 = 1;
double U_1 = 2;
double V = 0.5;
Model model;
model << HoppingAmplitude(U_0, {0}, {0});
model << HoppingAmplitude(U_1, {1}, {1});
model << HoppingAmplitude(V, {0}, {1}) + HC;

We first make a note about the terminology used in TBTK. From the canonical form for the Hamiltonian, and if we write the time evolution of a state as

\(\Psi(t+dt) = (1 - iHdt)\Psi(t)\),

it is clear that \(a_{\mathbf{i}\mathbf{j}}\) is the amplitude associated with the process where an electron is removed from \(\mathbf{j}\) and inserted at \(\mathbf{i}\). That is, in tight-binding terminology, the amplitude associated with a particle hopping from \(\mathbf{j}\) to \(\mathbf{i}\). While the term hopping amplitude often is restricted to the case where the two indices actually differ from each other, in TBTK it is used to refer to any \(a_{\mathbf{i}\mathbf{j}}\). Moreover, the indices \(\mathbf{i}\) and \(\mathbf{j}\) are often referred to as to- and from-indices, respectively, and a HoppingAmplitude is created as

HoppingAmplitude(value, toIndex, fromIndex);

In the example above two such HoppingAmplitudes are created for the diagonal entries and added to the Model using the <<-operator. Similarly, in the second last line a HoppingAmplitude that corresponds to \(a_{01}\) is created and added to the Model. However, on close inspection we see that + HC is added at the end of the line, which means that also the Hermitian conjugate is added. The last line that reads


creates a mapping between the indices that have been added to the Model and a linear Hilbert space index. It is important that the this call is performed after the HoppingAmplitudes has been added to the Model and that no more HoppingAmplitudes are added after this point.

Physical indices and Hilbert space indices

Although the above example is too simple to make the point since the index structure is trivial, we know from the Indices chapter that TBTK allows also much more complicated index structures to be used. One of the core features of TBTK is embodied in the call to model.construct(). To understand why, we note that although the indices in the problem above have a linear structure, this is rarely the case in general. In other problems, such as for example a two-dimensional lattice the index structure is not linear but is easily linearized by a defining a mapping such as \(h = L_y x + y\). The indices \((x, y)\) or {x, y} in TBTK notation we will call physical indices, while linearized indices such as \(h\) are referred to as Hilbert space indices.

So why do we need linearization? One answer is that algorithms almost universally are going to be most efficient in a linear basis. Moreover, algorithms can also be made much more general if the linearization procedure is not part of the algorithm itself. If mappings such as \(h = L_y x + y\) are hard coded into the Solver, it is essentially a single purpose code, possibly with a few knobs to turn to adjust the system size and parameter values, but the code is locked to a very specific type of problem. On the contrary, a Solver that accepts any Model that can provide itself on a linearized form can be reused for widely different types of problems.

So why physical indices? Is it not better if the developer performs the linearization from the start to improve performance? These questions would have been relevant if performance was an issue. The mapping from physical indices to Hilbert space and back again is certainly not a computationally cheap task if implemented through a naive approach such as a lookup table. In fact, this would have been prohibitively expensive. One of the most (or simply the most) important core data structures of TBTK is a tree structure in which the HoppingAmplitudes are stored. It allows for conversion back and forth between physical indices and Hilbert space indices with an overhead cost that is negligible in most parts of the code. Taking the responsibility to linearize indices of from the developer is a significant abstraction that allows mental capacity to be focused on physics instead of numerics.

The one place where the overhead really is important is in the Solvers where most of the computational time usually is spent. Solvers therefore usually internally convert the Model to some much more optimal format using the linearization provided by the Model before starting the main calculation. This is also the reason why, as described in the PropertyExtractor chapter, it is recommended to not use the Solvers immediately to extract properties, but to instead wrap them in PropertyExtractors. The PropertyExtractors provides the interface through which properties can be extracted from the Solvers using physical indices. In short, the distinction between physical indices and Hilbert space indices allows application developers to focus on the physics of the particular problem, while simultaneously allowing Solver developers to focus on numerical details and general properties rather than system specific details.

Example: A slightly more complex Model

To appreciate the ease with which Models can be setup in TBTK, we now consider a slightly more complex example. Consider a magnetic impurity on top of a two-dimensional substrate of size 51x51 described by

\(H = H_{S} + H_{Imp} + H_{Int}\)


\begin{eqnarray*} H_{S} &=& U_S\sum_{\mathbf{i}\sigma}c_{\mathbf{i}\sigma}^{\dagger}c_{\mathbf{i}\sigma} - t\sum_{\langle\mathbf{i}\mathbf{j}\rangle\sigma}c_{\mathbf{i}\sigma}^{\dagger}c_{\mathbf{j}\sigma},\\ H_{Imp} &=& (U_{Imp} - J)d_{\uparrow}^{\dagger}d_{\uparrow} + (U_{Imp} + J)d_{\downarrow}^{\dagger}d_{\downarrow},\\ H_{Int} &=& \delta\sum_{\sigma}c_{(25,25)\sigma}^{\dagger}d_{\sigma} + H.c. \end{eqnarray*}

Here \(\mathbf{i}\) is a two-dimensional index, \(\langle\mathbf{i}\mathbf{j}\rangle\) indicates summation over nearest neighbors, \(\sigma\) is a spin index, and \(c_{\mathbf{i}}\) and \(d_{\sigma}\) are operators on the substrate and impurity, respectively. The parameters \(U_S\) and \(U_{Imp}\) are on-site energies on the substrate and impurity, respectively, while \(t\) is a nearest neighbor hopping amplitude, \(J\) is a Zeeman term, and \(\delta\) is the coupling strength between the substrate and impurity.

We first note that an appropriate index structure is {0, x, y, s} for the substrate and {1, s} for the impurity. Using this index structure we next tabulate the hopping parameters on the canonical form given at the beginning of this chapter.

Hopping amplitude Value To Index From Index
\(a_{(0,x,y,\sigma),(0,x,y,\sigma)}\) \(U_S\) {0, x, y, s} {0, x, y, s}
\(a_{(0,x+1,y,\sigma),(0,x,y,\sigma)}\) \(-t\) {0, x+1, y, s} {0, x, y, s}
\(a_{(0,x,y,\sigma),(0,x+1,y,\sigma)}\) \(-t\) {0, x, y, s} {0, x+1, y, s}
\(a_{(0,x,y+1,\sigma),(0,x,y,\sigma)}\) \(-t\) {0, x, y+1, s} {0, x, y, s}
\(a_{(0,x,y,\sigma),(0,x,y+1,\sigma)}\) \(-t\) {0, x, y, s} {0, x, y+1, s}
\(a_{(1,\sigma),(1,\sigma)}\) \(U_{Imp}\) {1, s} {1, s}
\(a_{(1,\uparrow),(1,\uparrow)}\) \(-J\) {1, 0} {1, 0}
\(a_{(1,\downarrow),(1,\downarrow)}\) \(J\) {1, 1} {1, 1}
\(a_{(0,25,25,\sigma),(1,\sigma)}\) \(\delta\) {0, 25, 25, s} {1, s}
\(a_{(1,\sigma),(0,25,25,\sigma)}\) \(\delta\) {1, s} {0, 25, 25, s}

Symbolic subindices should be understood to imply that the values are valid for all possible values of the corresponding subindices. We also note that hopping amplitudes that appear multiple times should be understood to add the final value. For example does \(a_{(1,\uparrow),(1,\uparrow)}\) appear twice (sixth and seventh row) and should be understood to add to \(U_{Imp} - J\). While the first column is the analytical representation of the symbol for the hopping amplitudes, the third and fourth column is the corresponding numerical representation. In particular, we note that we use 0 to mean up spin and 1 to mean down spin. Next we note that the table can be reduced if we take into account that row 2 and 3, 4 and 5, and 9 and 10 are each others Hermitian conjugates. Further, row 7 and 8 can be combined into a single row by writing the value as \(-J(1 - 2s)\). The table can therefore if we also ignore the first column be compressed to

Value To Index From Index Add Hermitian conjugate
\(U_S\) {0, x, y, s} {0, x, y, s}
\(-t\) {0, x+1, y, s} {0, x, y, s} Yes
\(-t\) {0, x, y+1, s} {0, x, y, s} Yes
\(U_{Imp}\) {1, s} {1, s}
\(-J(1 - 2s)\) {1, 0} {1, 0}
\(\delta\) {0, 25, 25, s} {1, s} Yes

Once on this form, it is simple to implement the Model as follows

const int SIZE_X = 51;
const int SIZE_Y = 51;
double U_S = 1;
double U_Imp = 1;
double t = 1;
double J = 1;
double delta = 1;
Model model;
//Setup substrate.
for(int x = 0; x < SIZE_X; x++){
for(int y = 0; y < SIZE_Y; y++){
for(int s = 0; s < 2; s++){
model << HoppingAmplitude(U_S, {0, x, y, s}, {0, x, y, s});
if(x+1 < SIZE_X)
model << HoppingAmplitude(-t, {0, x+1, y, s}, {0, x, y, s}) + HC;
if(y+1 < SIZE_Y)
model << HoppingAmplitude(-t, {0, x, y+1, s}, {0, x, y, s}) + HC;
for(int s = 0; s < 2; s++){
//Setup impurity.
model << HoppingAmplitude( U_Imp, {1, s}, {1, s});
model << HoppingAmplitude(-J*(1-2*s), {1, s}, {1, s});
//Add coupling between the substrate and impurity.
model << HoppingAmplitude(delta, {0, SIZE_X/2, SIZE_Y/2, s}, {1, s}) + HC;
//Construct model.

For pedagogical reasons we here went through the process of converting the analytical Hamiltonian to a numerical Model in quite some detail. However, with a small amount of experience the second table can be read of immediately from the Hamiltonian, cutting down the work significantly. A key advice is to utilize the Hermitian conjugate to their maximum like we did above. Note in particular that we used this to only have \(x+1\) and \(y+1\) in one position for the indices, respectively (and no \(x-1\) or \(y-1\)). Doing so reduces the number of lines of code, improves readability, simplifies maintainance, and consequently reduces the chance of introducing errors.

Advanced: Using IndexFilters to construct a Model

While the already introduced concepts significantly simplifies the modeling of complex geometries, TBTK provides further ways to simplify the modeling stage. In particular, we note that in the example above, conditional statements had to be used in the first for-loop to ensure that HoppingAmplitudes were not added across the boundary of the system. For more complex structures it is useful to be able to separate the specification of such exceptions to the rule from the specification of the rule itself. This can be accomplished through the use of an IndexFilter, which can be used by the Model to accept or reject HoppingAmplitudes based on their to- and from-indices.

In this example we will show how to setup a simple two-dimensional sheet like the substrate above and in addition add a hole in the center of the substrate. We will here assume that the size of the substrate and the radius of the hole has been specified using three global variables SIZE_X, SIZE_Y, and RADIUS. A good implementation would certainly remove the need for global variables, but we use this method here to highlight the core concepts since removing the global variables requires a bit of extra code that obscures the key point. Moreover, for many applications such parameters would be so fundamental to the calculation that it actually may be justified to use global variables to store them.

The first step toward using an IndexFilter is to create one. The syntax for doing so is as follows.

class Filter : public AbstractIndexFilter{
Filter* clone() const{
return new Filter();
bool isIncluded(const Index &index){
double radius = sqrt(
pow(index[0] - SIZE_X/2, 2)
+ pow(index[1] - SIZE_Y/2, 2)
if(radius < RADIUS)
return false;
else if(index[0] < 0 || index[0] >= SIZE_X)
return false;
else if(index[1] < 0 || index[1] >= SIZE_Y)
return false;
return true;

The experienced C++ programmer recognizes this as an implementation of an abstract base class and is encouraged to write more complicated filters. In this case we note that it is important that the function clone() returns a proper copy of the Filter, which in this case is trivial since there are no member variables. However, for the developer not familiar with such concepts, it is enough to view this as a template where the main action takes place in the function isIncluded(). This function is responsible for returning whether a given input index is included in the Model or not. As seen we first calculate the radius for the Index away from the center of the sheet. Next we check whether the calculated radius is inside the specified RADIUS for the hole, or if it falls outside the boundaries of the sample. If either of these are true we return false, while otherwise we return true.

Note: When writing a Filter it is important to take into account not only the Index structure of Indices that are going to be rejected, but all Indices that are passed to the Model at the point of setup. If there for example would have been some Indices in the Model that only had one subindex, the Filter above would have been invalid since it may be passed such an Index and try to access index[1]. When writing filters it is therefore important to ensure they are compatible with the Model as a whole and are able to respond true or false for every possible Index in the Model. Since in this case we model a two-dimensional sheet where all Indices have more than two subindices this Filter is alright.

Once the Filter is specified, we are ready to use it to setup a Model

double U_s = 1;
double t = 1;
Model model;
Filter filter;
for(int x = 0; x < SIZE_X; x++){
for(int y = 0; y < SIZE_Y; y++){
for(int s = 0; s < 2; s++){
model << HoppingAmplitude(U_s, { x, y, s}, {x, y, s});
model << HoppingAmplitude( -t, {x+1, y, s}, {x, y, s}) + HC;
model << HoppingAmplitude( -t, { x, y+1, s}, {x, y, s}) + HC;

Advanced: Callback functions

Sometimes it is useful to be able to delay specification of a HoppingAmplitudes value to a later time than the creation of the Model. This is for example the case if the same Model is going to be solved multiple times for different parameter values, or if some of the parameters in the Hamiltonian are going to be determined self-consistently. For this reason it is possible to pass a function as value argument to the HoppingAmplitude rather than a fixed value. If we have indices with the structure {x, y, s}, where the last index is a spin, and there exists some global parameter \(J\) that determines the current strength of the Zeeman term, a typical callbackFunction looks like

complex<double> callbackJ(const Index &to, const Index &from){
//Get spin index.
int s = from[2];
//Return the value of the HoppingAmplitude
return -J*(1 - 2*s);

Just like when writing an IndexFilter, a certain amount of Model specific information needs to go into the specification of the callbacks. Here we have for example chosen to determine the spin by looking at the 'from'-Index, which should not differ from the spin-index of the to-Index. However, unlike when writing IndexFilters, the HoppingAmplitude callbacks only need to make sure that they work for the restricted Indices for which the callback is fed into the Model together with, since these are the only Indices that will ever be passed to the callback.

Once the callback is defined, it is possible to use it when setting up a model as follows

model << HoppingAmplitude(callbackJ, {x, y, s}, {x, y, s});

Block structure

One of the most powerful methods for solving quantum mechanical problems is block diagonalization of the Hamiltonian. What this means is that the Hamiltonian is broken apart into smaller blocks that are decoupled from each other and thereby effectively replaces one big problem with many smaller problems. The smaller problems are much simpler to solve and usually far outweighs the fact that several problems have to be solved instead. The most obvious example of this is a Hamiltonian which when written in reciprocal space separates into independent problems for each momentum \(\mathbf{k}\).

To facilitate the exploitation of block diagonal Hamiltonians TBTK has restricted support for automatically detecting when it is handed a Hamiltonian with a block diagonal structure. However, for this to work the developer has to organize the Index structure in such a way that TBTK is able to automatically take advantage of the block diagonal structure. Namely, TBTK will automatically detect existing block structures as long as the Index has the subindices that identifies the independent blocks standing to the left of the intra block indices. For a system with say three momentum subindices, one orbital subindex, and one spin subindex, where the Hamiltonian is block diagonal in the momentum space subindices, an appropriate Index structure is {kx, ky, kz, orbital, spin}. If for some reason the Index structure instead is given as {kx, ky, orbital, kz, spin}, TBTK will only be able to automatically detect kx and ky as block-indices, and will treat all of the three remaining subindices orbital, kz, and spin as internal indices for the blocks (at least as long as the Hamiltonian does not happen to also be block diagonal in the orbital subindex).