Using TBTK to calculate the density of states (DOS) of a 1D, 2D, and 3D square lattice

Kristofer Björnson

(23rd of October 2018)

Most recent TBTK release at the time of writing: v1.0.2 (Updated to work with: v2.0.0)

In this post we will walk through how to calculate and plot the density of state (DOS) of a 1D, 2D, and 3D square lattice, using TBTK. To achieve this we will first implement three functions for creating the corresponding models. We will then create a fourth function in which the DOS is calculated for each of the three models.

This example also provides a simple demonstration of how TBTK naturally support the creation of general purpose code, in this case the calculation of the DOS, that can be used independently of the details of the model.

1 Creating the Model

1.1 One dimension

In one dimension the Hamiltonian is given by

\displaystyle H_{1D}

\displaystyle=-2t\sum_{k_{x}}\cos(k_{x})c_{k_{x}}^{\dagger}c_{k_{x}},

(1)

where t is a hopping parameter that we will set to one.

To create the model we loop over the summation index and feed each of the hopping amplitudes to the Model. Note how we let the numerical indices run over the positive discrete values 0 to 9999, while the argument of the cos function falls in the interval $[-\pi,\pi]$ .

⬇

1Model createModel1D(){

2 //Parameters.

3 const int SIZE_X = 10000;

4 double t = 1;

6 //Create the Model.

7 Model model;

8 for(int kx = 0; kx < SIZE_X; kx++){

9 double KX = 2*M_PI*kx/(double)SIZE_X - M_PI;

10 model << HoppingAmplitude(

11 -2*t*cos(KX),

12 {kx},

13 {kx}

14 );

15 }

17 return model;

18}

1.2 Two dimensions

In two dimension the Hamiltonian is given by

\displaystyle H_{2D}

\displaystyle=-2t\sum_{k_{x},k_{y}}\left(\cos(k_{x})+\cos(k_{y})\right)c_{% \mathbf{k}}^{\dagger}c_{\mathbf{k}},

(2)

where $\mathbf{k}=(k_{x},k_{y})$ .

The corresponding Model is created as follows using a $500\times 500$ grid.

⬇

1Model createModel2D(){

2 //Parameters.

3 const int SIZE_X = 500;

4 const int SIZE_Y = 500;

5 double t = 1;

7 //Create the Model.

8 Model model;

9 for(int kx = 0; kx < SIZE_X; kx++){

10 for(int ky = 0; ky < SIZE_Y; ky++){

11 double KX = 2*M_PI*kx/(double)SIZE_X;

12 double KY = 2*M_PI*ky/(double)SIZE_Y;

14 model << HoppingAmplitude(

15 -2*t*(cos(KX) + cos(KY)),

16 {kx, ky},

17 {kx, ky}

18 );

19 }

20 }

22 return model;

23}

1.3 Three dimensions

In three dimension the Hamiltonian is given by

\displaystyle H_{3D}

\displaystyle=-2t\sum_{k_{x},k_{y},k_{z}}\left(\cos(k_{x})+\cos(k_{y})+\cos(k_% {z})\right)c_{\mathbf{k}}^{\dagger}c_{\mathbf{k}},

(3)

where $\mathbf{k}=(k_{x},k_{y},k_{z})$ .

The corresponding Model is created as follows using a $200\times 200\times 200$ grid.

⬇

1Model createModel3D(){

2 //Parameters.

3 const int SIZE_X = 200;

4 const int SIZE_Y = 200;

5 const int SIZE_Z = 200;

6 double t = 1;

8 //Create the Model.

9 Model model;

10 for(int kx = 0; kx < SIZE_X; kx++){

11 for(int ky = 0; ky < SIZE_Y; ky++){

12 for(int ky = 0; ky < SIZE_Y; ky++){

13 double KX = 2*M_PI*kx/(double)SIZE_X;

14 double KY = 2*M_PI*ky/(double)SIZE_Y;

15 double KZ = 2*M_PI*kz/(double)SIZE_Z;

17 model << HoppingAmplitude(

18 -2*t*(cos(KX) + cos(KY) + cos(KZ)),

19 {kx, ky, kz},

20 {kx, ky, kz}

21 );

22 }

23 }

24 }

26 return model;

27}

2 Calculating the DOS

We are now ready to implement the main function for calculating the DOS. To take advantage of that the Hamiltonian consists of independent blocks for each $\mathbf{k}$ value, we will use the Solver::BlockDiagonalizer. Knowing that the density of states will lie in the interval $[-6,6]$ for all three Models, we next setup the PropertyExtractor to work on the energy interval $[-7,7]$ and calculate the DOS. Finally we plot the DOS, using gaussian smoothin with $\sigma=0.05$ and a convolution window of 101 energy points and save the results to file.

⬇

1int main(int argc){

2 //Initialize TBTK.

3 Initialize();

5 //Filenames to save the figures as.

6 string filenames[3] = {

7 ”figures/DOS_1D.png”,

8 ”figures/DOS_2D.png”,

9 ”figures/DOS_3D.png”,

10 };

12 //Loop over 1D, 2D, and 3D.

13 for(int n = 0; n < 3; n++){

14 //Create the Model.

15 Model model;

16 switch(n){

17 case 0:

18 model = createModel1D();

19 break;

20 case 1:

21 model = createModel2D();

22 break;

23 case 2:

24 model = createModel3D();

25 break;

26 default:

27 Streams::out << ”Error: Invalid case value.\n”;

28 exit(1);

29 }

31 //Setup and run the Solver.

32 Solver::BlockDiagonalizer solver;

33 solver.setModel(model);

34 solver.run();

36 //Setup the PropertyExtractor.

37 PropertyExtractor::BlockDiagonalizer propertyExtractor(solver);

38 propertyExtractor.setEnergyWindow(-7, 7, 1000);

39 Property::DOS dos = propertyExtractor.calculateDOS();

41 //Normalize the DOS.

42 for(unsigned int c = 0; c < dos.getResolution(); c++)

43 dos(c) = dos(c)/model.getBasisSize();

45 //Smooth the DOS.

46 const double SMOOTHING_SIGMA = 0.05;

47 const unsigned int SMOOTHING_WINDOW = 101;

48 dos = Smooth::gaussian(dos, SMOOTHING_SIGMA, SMOOTHING_WINDOW);

50 //Plot and save the result.

51 Plotter plotter;

52 plotter.plot(dos);

53 plotter.save(filenames[n]);

54 }

56 return 0;

57}

3 Results

Below we show the result of the calculation for the three different cases.

3.1 One dimensional

3.2 Two dimensional

3.3 Three dimensional

4 Full code

Full code is available in src/main.cpp in the project 2018_10_23 of the Second Tech code package. See the README for instructions on how to build and run.