Dynamically adjustable HoppingAmplitudes using callback functions in TBTK

Most recent TBTK release at the time of writing: v1.0.3

In this blog post we will take a look at how to specify models without necessarily knowing all the model parameters at the point of specification. This can be very useful when the Hamiltonian is time dependent or depends on some parameter that is to be solved for self-consistently.

We will here consider a third situation in which adjustable HoppingAmplitudes can be utilized by setting up a model with a variable potential. We will use this to repeatedly calculate the probability density and eigenvalues for the lowest energy states for seven different types of one-dimensional potentials. Namely, the infinite square well, square well, double square well, harmonic oscillator, double well, and step potential, as well as a potential with two regions with different potential separated by a barrier.


The parameters that are needed are the number of sites to use for the discretized Schrödinger equation, the number of states for which we are going to calculate the probability density, and the nearest neighbor hopping amplitude t. We define these as follows.

//Number of sites and states to extract the
//probability density for.
const int SIZE_X = 500;
const int NUM_STATES = 7;
double t = 1;

We further create an enum that will be used to specify the type of the potential.

//Enum type for specifying the type of the
enum PotentialType{

//The type of the potential.
PotentialType potentialType;

Callback function

Normally a HoppingAmplitude is specified as

HoppingAmplitude(value, toIndex, formIndex);

However, it is also possible to pass a function instead of a value as first argument. This function will later be called by the HoppingAmplitude to determine the value, passing the toIndex and fromIndex as arguments. That is, the HoppingAmplitude will “call back” later to determine the value.

For a function to classify as a valid callback function, it needs to take two indices as arguments and return a complex. In addition, it also needs to be able to handle every pair of indices that it will be called with. That is, every pair of indices which it is packed together with into a HoppingAmplitude when the model is specified.

In our case we will deal with a one-dimensional system and since the callback is used for the potential term, it will only be packed together with indices that are equal with each other. When implementing the callback for the potential, we therefore only need to consider one of the two indices and can extract the position from its first (and only) subindex.

//Callback that dynamically returns the
//current potential on a given site.
complex<double> potential(
	const Index &toIndex,
	const Index &fromIndex
	int x = fromIndex[0];

	case InfiniteSquareWell:
		return infiniteSquareWell(x);
	case SquareWell:
		return squareWell(x);
	case DoubleSquareWell:
		return doubleSquareWell(x);
	case HarmonicOscillator:
		return harmonicOscillator(x);
	case DoubleWell:
		return doubleWell(x);
	case Step:
		return step(x);
	case Barrier:
		return barrier(x);
		Streams::out << "Error. This should never happen.";

After having extracted the position from the fromIndex, this function calls one of seven potential functions depending on the current potentialType to determine the value of the potential at the given point. The implementation of the potential functions is not explicitly shown here, but can be found in the complete example linked from the bottom of this post.


The model that we will consider is the discretized one-dimensional Schrödinger equation which we specify as follows.1

//Create the Model.
Model model;
for(unsigned int x = 0; x < SIZE_X; x++){
	//Kinetic terms.
	model << HoppingAmplitude(
	if(x + 1 < SIZE_X){
		model << HoppingAmplitude(
			{x + 1},
		) + HC;

	//Potential term.
	model << HoppingAmplitude(

The kinetic terms result from discretizing the Laplace operator and letting t = \hbar^2/2mdx^2, while the potential term has been added by passing the callback function potential() as the first parameter to the HoppingAmplitude.

Set up the Solver and PropertyExtractor

We next set up the Solver and PropertyExtractor as follows.

//Setup the Solver and PropertyExtractor.
Solver::Diagonalizer solver;

The main loop

In the main loop, we will carry out the same calculation for each potential. We therefore begin by creating a list of potential types to perform the calculations for, and a list of filenames to which to save the results.

//List of potentials to run the calculation
vector<PotentialType> potentialTypes = {

//List of filenames to save the results to.
vector<string> filenames = {

The main loop itself is as follows.

//Run the calculation and plot the result for
//each potential.
	unsigned int n = 0;
	n < potentialTypes.size();
	//Set the current potential.
	potentialType = potentialTypes[n];

	//Run the solver.

	//Calculate the probability density
	//for the first NUM_STATES.
	Array<double> probabilityDensities(
		unsigned int state = 0;
		state < NUM_STATES;
			unsigned int x = 0;
			x < (unsigned int)SIZE_X;
			complex<double> amplitude
				= propertyExtractor.getAmplitude(
			probabilityDensities[{state, x}]
				= pow(abs(amplitude), 2);

	//Calculate the eigenvalues.
	Property::EigenValues eigenValues
		= propertyExtractor.getEigenValues();

	//Plot the results.

Here we loop over each type of potential and set the global parameter potentialType to the current potential in line 9. The solver is then executed, which results in the HoppingAmplitudes being requested internally, which triggers the callback functions to be called and return the values for the current potential. Next, the probability density is calculated2, followed by the calculation of the eigenvalues.

Finally, the results are plotted and saved to file. The plotting routine involves significant post processing of the data and the full code for achieving this can be found in the code linked below. In particular, the probability densities are first rescaled such that NUM_STATES probability densities can be stacked evenly on top of each other without overlapping or extending beyond the plotting window. Then the probability densities are shifted by the respective eigenvalue of the corresponding state to simultaneously convey information about the probability density and the eigenvalues in relation to the potential.


Below the results are presented. The thick red line indicates the potential, while the black lines are the probability densities. Each probability density has been offset by the corresponding states eigenvalue to put its energy in relation to the potential. The eigenvalues are also plotted with thin red lines. In agreement with the basic principles of quantum mechanics, we see how the states are oscillating in the regions where the potential is smaller than the energy of the state, while they decay exponentially into the regions where the potential is larger.

Infinite square well

Square well

Double square well

Harmonic oscillator

Double well



Full code

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

  1. For more details see for example Using TBTK to calculate the density of states (DOS) of a 1D, 2D, and 3D square lattice and Direct access to wave function amplitudes and eigenvalues in TBTK
  2. See Direct access to wave function amplitudes and eigenvalues in TBTK for a more detailed description of the calculation of the probability density.

Leave a Reply

Your email address will not be published. Required fields are marked *