How can I create such a 3D plot using matplotlib? Like the one in Figure (a) below: enter image description here

I want to create it using DFT bandstructure calculation result

  • $\begingroup$ I don't know if this would help you but DFT calculations usually outputs energy eigenvalues at each k point defined by the k point mesh. $\endgroup$ Commented May 15, 2021 at 19:39

3 Answers 3


Consider a graphene Hamiltonian, whose dispersion looks a lot like the one in your figure. Per these notes, its k-space Hamiltonian may be written as: $$ H(k)=-t \sum_{\delta} [\cos(k\cdot\delta)\sigma_x-\sin(k\cdot\delta)\sigma_y], $$ where $k=(k_x,k_y)$, $\delta$ are nearest-neighbor vectors, $\sigma_i$ are Pauli matrices and $t$ is a hopping parameter. Please see referenced notes for vectors and other details pertaining to the model.

General idea:

Per the time-independent Schrodinger equation, the energy dispersion you seek to plot comprises just the eigenvalues of the above Hamiltonian over some region of $k$-space. I will assume that you are working with an expression for the Hamiltonian (like the one above), instead of ab initio techniques. The given $2\times 2$ Hamiltonian is a $2$-level Hamiltonian that yields $2$ instantaneous eigenvectors/wavefunctions corresponding to $2$ eigenvalues. In your figure, the sets of eigenvalues pertaining to each eigenvector are given in two different colors (red and blue).

Therefore, you simply have to solve for the Hamiltonian's eigenvalues at each $k$-point in $k$-space. You may do this 'analytically' (to find expressions for the eigenvalues using algebra, that you can then plug $k$ points into), or 'numerically' (evaluating the Hamiltonian at each $k$ point first, to numerical get eigenvalues directly).

Analytic implementation:

The 'analytic' method is just using some function $H(k)$ for the Hamiltonian to get an analytic expression for a different function $E(k)$ of the energy. You can plug $k$ into $E(k)$ to plot your figure. I personally prefer using Mathematica for this, because I like its in-build eigensystem-solving commands.

(* Define parameter space edges *)
(* Define Hamiltonian *)
a = 2.46; t = 1; k = {kx, ky}; 
d1 = (a/2) {1, Sqrt[3]}; d2 = (a/2) {1, -Sqrt[3]}; d3 = -a {1, 
sx = {{0, 1}, {1, 0}}; sy = {{0, -I}, {I, 0}}; (* Pauli matrices *)
hx = Cos[k.d1] + Cos[k.d2] + Cos[k.d3]; hy = 
 Sin[k.d1] + Sin[k.d2] + Sin[k.d3];
H[kx_, ky_] = -t (hx*sx - hy*sy);
(* Get analytic expressions for energies *)
eigenvals[kx_, ky_] = Eigenvalues[H[kx, ky]];
(* Plot *)
Show[Plot3D[eigenvals[kx, ky][[1]], {kx, BZx1, BZx2}, {ky, BZy1, BZy2}, 
  PlotStyle -> Red], 
 Plot3D[eigenvals[kx, ky][[2]], {kx, BZx1, BZx2}, {ky, BZy1, BZy2}, 
  PlotStyle -> Blue]]

This gives you:

enter image description here

I can imagine that a Sympy (+ matplotlib) approach would be similar. Additionally, I would not recommend trying to solve complicated expressions this way (unless you have the computational power).

Numerical implementation:

On the other hand, the 'numerical' method requires first discretizing your $k$-space into a grid, and then plugging each grid point into $H(k)$ to get a matrix with numbers as its entries (as opposed to algebraic expressions). You can easily solve for the eigenvalues, and plot all resulting eigenvalue points to get a figure like the one you seek. Please bear in mind that eigenvectors are generally complex-valued, and that when using the numerical approach, different programs assign different arbitrary, complex-valued phases to each entry. This is not an issue in our case, as we care only about eigenvalues (which are real if dealing with a Hermitian system). I personally prefer using MATLAB or Python for this approach, because they are well-suited for working with numerical datasets.

The following code might look more complicated than it needs to be, but I like working with functions:

import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

# First, define your Hamiltonian. 
# The following function returns the Hamiltonian matrix for input (kx,ky).
def Hamiltonian_Graphene(kx,ky):
    a = 2.46 # Wikipedia
    t = 1
    k = np.array((kx,ky))
    d1 = (a/2)*np.array((1,np.sqrt(3)))
    d2 = (a/2)*np.array((1,-np.sqrt(3)))
    d3 = -a*np.array((1,0))
    sx = np.matrix([[0,1],[1,0]]);
    sy = np.matrix([[0,-1j],[1j,0]]);
    hx = np.cos(k@d1)+np.cos(k@d2)+np.cos(k@d3)
    hy = np.sin(k@d1)+np.sin(k@d2)+np.sin(k@d3)
    H = -t*(hx*sx - hy*sy)
    return H

# Onto plotting the dispersion.
# The parameters m,n specify the resolution to discretize the parameter space into (e.g. break it into m=10 pieces in the kx direction).
# The parameters BZx1,BZx2,BZy1,BZy2 specify the edges of your parameter space. I used BZx1=-1,BZx2=+1,BZy1=-1,BZy2=+1 in the 
def plot_dispersion(m,n,BZx1,BZx2,BZy1,BZy2):
    # Generate a mesh
    kx_range = np.linspace(BZx1, BZx2, num=m)
    ky_range = np.linspace(BZy1, BZy2, num=n)
    # Get the number of levels with a dummy call (an NxN square matrix has N levels)
    num_levels = len(Hamiltonian_Graphene(1,1))
    energies = np.zeros((m,n,num_levels)); # initialize
    # Now iterate over discretized mesh, to consider each coordinate.
    for i in range(m):
        for j in range(n):
            H = Hamiltonian_Graphene(kx_range[i],ky_range[j]); 
            evals, evecs = LA.eig(H); # Numerically get eigenvalues and eigenvectors
    X, Y = np.meshgrid(kx_range, ky_range) # Generate actual mesh for plotting.
    # Plot! There are several ways to style this.
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.axis([X.min(), X.max(), Y.min(), Y.max()])
    transparency = 0.3 #transparency
    for band in range(num_levels): 
        ax.scatter3D(X, Y, energies[:,:,band], alpha=transparency, antialiased=False, label="c") 
# Run it!

This gives you:

enter image description here

You might get an error ComplexWarning: Casting complex values to real discards the imaginary part, but that is because of rounding errors in Python (and can be avoided if you use np.eigh() instead of np.eig, if applicable). You can try to troubleshoot ways to get rid of this, and to make resulting figures prettier. I did not include the contour plot shown at the bottom of your figure, as it basically uses the same procedure.

  • $\begingroup$ Thanks for the nice answer. But I'm really sorry that, I try to create this plot using DFT bandstructure result. So, I'm trying to create this 3D plot and its 2D BZ energy level from a DFT band structure calculation result $\endgroup$
    – Alfred
    Commented May 14, 2021 at 22:43
  • 1
    $\begingroup$ @Alfred, No worries. It might be worthwhile specifying the details in your question and adding appropriate tage. Including the software and output data file format you are working with. $\endgroup$ Commented May 14, 2021 at 22:47


VASPKIT offers a tool to plot 3D band structure using VASP. The method select K-path on a surface of the Irreducilbe Brillouin zone and calculate the K-dependent band energies on those K-points. An example 3D plot of graphene bandstructure is shown below enter image description here

For more details refer this tutorial

A few more plots enter image description here

  • 1
    $\begingroup$ Very nice plots! $\endgroup$ Commented May 15, 2021 at 21:08
  • 1
    $\begingroup$ It took me a bit of time to figure out what these ♪ in the last image were supposed to mean :) $\endgroup$
    – Ruslan
    Commented May 15, 2021 at 21:16
  • $\begingroup$ Tye plots are directly from the tutorial pages $\endgroup$
    – Thomas
    Commented May 16, 2021 at 0:56
  • $\begingroup$ @Ruslan I believe they are called eighth notes. $\endgroup$ Commented May 16, 2021 at 3:19
  • $\begingroup$ I really liked last plots. Are they also created using vaspkit? In the tutorial, it just shows the one for Graphene. @Thomas $\endgroup$
    – Alfred
    Commented May 16, 2021 at 5:16

For a general case, create you input inside a loop:

  1. Define the starting k point.
  2. Generate the k-points path.
  3. Export your input.
  4. Call and execute your DFT (using the input file generated in previous step).
  5. Read the output in the corresponding band structure file. Add the band data (1D array) to a column in a matrix.
  6. Change the starting k-point.
  7. Back to step 1.

Here, the mains steps are completely dependent on the DFT software you are using.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .