See more details about the Model in the API

The Model class

The Model class is a container for the Hamiltonian, chemical potential, temperature, statistics, etc. It is best understood through an example, therefore consider a 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 can be written

\(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\). Further, assume that we will study the system at temperature T = 300K, chemical potential mu = 0, and for particles with Fermi-Dirac statistics.

Letting U_0, U_1, and V be the numerical symbols for \(U_0, U_1\), and \(V\), we can implement the model as follows.

//Create the Model.
Model model;
 
//Set the temperature, chemical potential, and statistics.
model.setTemperature(300);
model.setChemicalPotential(0);
model.setStatistics(Statistics::FermiDirac);
 
//Specify the Hamiltonian.
model << HoppingAmplitude(U_0, {0}, {0});
model << HoppingAmplitude(U_1, {1}, {1});
model << HoppingAmplitude(V,   {0}, {1}) + HC;
 
//Construct a linearized basis for the Model.
model.construct();

Here the operator<< is used to pass the HoppingAmplitudes to the Model. Also note the use of "+ HC" in the line

model << HoppingAmplitude(V, {0}, {1}) + HC;

which is equivalent to

model << HoppingAmplitude(V, {0}, {1});

model << HoppingAmplitude(conj(V), {1}, {0});

The final call to model.construct() creates a linearized Hilbert space basis for the Indices passed into the Model.

A more complex example

To see how a more complex Model can be set up, consider a magnetic impurity on top of a two-dimensional substrate of size 51x51. We describe the substrate with the Hamiltonian

\(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},\)

where \(\mathbf{i}\) is a two-dimensional index, \(\sigma\) is a spin index, and \(\langle\mathbf{i}\mathbf{j}\rangle\) denotes summation over nearest neighbors. The impurity is described by

\(H_{Imp} = (U_{Imp} - J)d_{\uparrow}^{\dagger}d_{\uparrow} + (U_{Imp} + J)d_{\downarrow}^{\dagger}d_{\downarrow},\)

and is connected to the site (25, 25) in the substrate through

\(H_{Int} = \delta\sum_{\sigma}c_{(25,25)\sigma}^{\dagger}d_{\sigma} + H.c.\)

The total Hamiltonian is

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

We first note that an appropriate Index structure is {0, x, y, s} for the substrate and {1, s} for the impurity, where s means spin. Reading of the Hamiltonian terms above, we can tabulate the hopping amplitudes.

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}

The left column contains the analytic symbol, while the three other columns contain the corresponding information that is needed to create a HoppingAmplitude. Subindices that are not numeric values (or arrows) should be understood to be summed/looped over in the analytic/numeric expressions.

The line pairs (2 and 3), (4 and 5), and (9 and 10) are the Hermitian conjugates of each other. Further, line 7 and 8 can be combined into a single line if we write the amplitude as -J(1 - 2s). Taking this into account, the table can be compressed into

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

It is now straight forward to implement the Model.

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 the 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 the impurity.
    model << HoppingAmplitude(     U_Imp, {1, s}, {1, s});
    model << HoppingAmplitude(-J*(1-2*s), {1, s}, {1, s});
 
    //Add the coupling between the substrate and the impurity.
    model << HoppingAmplitude(
        delta,
        {0, SIZE_X/2, SIZE_Y/2, s},
        {1, s}
    ) + HC;
}
 
//Construct a linearized basis for the Model.
model.construct();

The conversion from the analytic Hamiltonian to the numeric Model was here done in quite some detail. However, with a bit of experience it is possible to read of the HoppingAmplitude parameters immediately from the Hamiltonian without having to write down any tables.

We finish this example with the advice to always utilize the Hermitian conjugate to its maximum, as we did above. Note in particular that we used this to only have \(x+1\) and \(y+1\) in the to-Index (and there are no \(x-1\) or \(y-1\)). In short, only forward hopping terms are explicit. This improves readability and reduces the chance of introducing errors.

Advanced: Using IndexFilters to construct a Model

In the example above, if-statements were added inside the first loop to prevent HoppingAmplitudes from being added across the boundary. The code can be made more readable by using an IndexFilter, which allows for the handling of such exceptions to be separated from the main Model specification.

We here demonstrate how an IndexFilter can be used to create a square shaped geometry with a hole inside. To do so, we begin by defining the IndexFilter.

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

The Filter needs to know the size of the the system and the holes radius, therefore the constructor takes these as arguments and stores them in its private variables. The clone()-function is required and should return a pointer to a copy of the Filter. Finally, the isIncluded()-function is responsible for returning true or false depending on whether the given Index should be included in the Model or not. In this example, it first checks whether the Index is inside the excluded radius and returns false if this is the case. Otherwise, it returns true or false depending on whether or not the Index is inside the range [0, sizeX)x[0, sizeY).

With the Filter defined, a Model with the size 51x51 and a hole radius of 10 can be set up as follows.

Model model;
model.setFilter(Filter(51, 51, 10));
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;
        }
    }
}
model.construct();

Note the call to model.setFilter() and the fact that no if-statements appear inside the loop.

We end this section with a note about substems with complex Index structures. Consider a Model with two subsystems Indexed by {0, x, spin} and {1, x, y, z, spin}, respectively. The IndexFilter needs to be able to handle both of these Index structures. For example, requesting index[4] in the isIncluded()-function before knowing that the Index belongs to subsystem 1 is an error. Therefore, an appropriate implementation of the isIncluded()-function would, in this case, look something like this.

bool isIncluded(const Index &index){
    int subsystem = index[0];
    switch(subsystem){
    case 0:
        int x = index[1];
        int spin = index[2];
        //Do the relevant checks for subsystem 0 here.
        //...
 
        break;
    case 1:
        int x = index[1];
        int y = index[2];
        int z = index[3];
        int spin = index[4];
        //Do the relevant checks for subsystem 1 here.
        //...
 
        break;
    }
}

Block structure

The Model can recognize a block diagonal Hamiltonian under the condition that the block diagonal Subindices appear to the left in the Index. Consider a translationally invariant system that is block-diagonal in kx, ky, and kz, but not in orbital or spin. If the Index structure {kx, ky, kz, orbital, spin} is used, the Model will be able to recognize the block structure. This can be utilized by some Solvers such as the BlockDiagonalizer to speed up calculations. If instead {kx, ky, orbital, kz, spin} is used, only the kx and ky Subindices will be recognized as block indices.

Next: Solvers