In this example, we will plot the electronic band structure and density of states (DOS) of Cu₂O at the DFT-PBE0/TZVP level of theory.

Prerequisites: You need to be in the directory that contains the wavefunction file Cu2O_Pn-3m_PBE0_TZVP_opt.w from the Geometry optimization of Cu2O (the wavefunction file attached here is for CRYSTAL17 and may not work in other versions). 

Table of contents

Band widths

First it makes sense to run a bandwidth calculation (BWIDTH) which yields the highest and lowest energies of all bands. This helps is to find out the relevant bands for the band structure calculation. We also include PPAN keyword in the input to get Mulliken population analysis at the same time (it is computationally very cheap).

Create the following property calculation input file Cu2O_Pn-3m_PBE0_TZVP_opt.d3 in the directory that contains the wavefunction file Cu2O_Pn-3m_PBE0_TZVP_opt.w.

NEWK
0 8
8 8 8
1 0
BWIDTH
0 0
END

Line 1, NEWK: Request CRYSTAL to calculate the Fock/Kohn-Sham eigenvectors at a number of k-points in the reciprocal space. The lines 0 8 and 8 8 8 are equal to the original SHRINK values (see the chapter on geometry optimization), the original k-mesh is enough to obtain the band widths. Line 1 0 requests CRYSTAL to calculate the Fermi energy using this k-mesh (1) and tells that no additional printing is required (0).

Line 5, BWIDTH: Option 0 0 requests analysis starting from the first valence band (0) and including 4 virtual bands (0). Note that the definition of valence bands is not very clear especially when some atoms are described using effective core potentials. 

Next, submit the job:

jsub -np 8 crystal Cu2O_Pn-3m_PBE0_TZVP_opt.d3

Input files for property calculations always have the extension .d3 and runcrys script uses file extension .prop.o for the resulting output.

The calculation finishes very quickly and produces the following output file Cu2O_Pn-3m_PBE0_TZVP_opt.prop.o:

 BAND LIMITS FROM BAND  39 TO BAND  70

 BAND  39 EMIN -9.725863E-01 EMAX -9.583699E-01
 BAND  40 EMIN -9.590734E-01 EMAX -9.454801E-01
 BAND  41 EMIN -4.597114E-01 EMAX -4.337285E-01
 BAND  42 EMIN -4.597114E-01 EMAX -4.337285E-01
 BAND  43 EMIN -4.597114E-01 EMAX -4.198848E-01
 BAND  44 EMIN -4.291165E-01 EMAX -3.688226E-01
 BAND  45 EMIN -4.113627E-01 EMAX -3.591029E-01
 BAND  46 EMIN -4.113627E-01 EMAX -3.591029E-01
 BAND  47 EMIN -3.591029E-01 EMAX -3.115155E-01
 BAND  48 EMIN -3.293013E-01 EMAX -3.115155E-01
 BAND  49 EMIN -3.293013E-01 EMAX -3.091945E-01
 BAND  50 EMIN -3.166618E-01 EMAX -3.016122E-01
 BAND  51 EMIN -3.166618E-01 EMAX -2.992606E-01
 BAND  52 EMIN -3.166618E-01 EMAX -2.984846E-01
 BAND  53 EMIN -3.018289E-01 EMAX -2.916771E-01
 BAND  54 EMIN -3.018289E-01 EMAX -2.902164E-01
 BAND  55 EMIN -2.975386E-01 EMAX -2.834116E-01
 BAND  56 EMIN -2.975386E-01 EMAX -2.834116E-01
 BAND  57 EMIN -2.901128E-01 EMAX -2.834116E-01
 BAND  58 EMIN -2.870673E-01 EMAX -2.781059E-01
 BAND  59 EMIN -2.861799E-01 EMAX -2.781059E-01
 BAND  60 EMIN -2.861799E-01 EMAX -2.756278E-01
 BAND  61 EMIN -2.796075E-01 EMAX -2.636825E-01
 BAND  62 EMIN -2.775669E-01 EMAX -2.594586E-01
 BAND  63 EMIN -2.743033E-01 EMAX -2.215842E-01
 BAND  64 EMIN -2.573174E-01 EMAX -1.888888E-01
 BAND  65 EMIN -2.268851E-01 EMAX -1.888888E-01
 BAND  66 EMIN -2.178131E-01 EMAX -1.888888E-01
 BAND  67 EMIN -1.039795E-01 EMAX -1.743018E-02
 BAND  68 EMIN -6.641208E-02 EMAX  1.245247E-02
 BAND  69 EMIN -6.641208E-02 EMAX  4.670397E-02
 BAND  70 EMIN  2.638270E-02 EMAX  2.190251E-01
 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT BWIDTH      TELAPSE        0.36 TCPU        0.36

From the output we can now see the band numbers and energies. Explanation:

• The bands 39–66 are occupied and the bands 67–70 are empty.
• The energies are in atomic units (Hartree). 
• The Hartree to eV conversion factor is 27.2114 (CODATA2018).
• These bands span energies from –26.5 eV (-9.725863E-01 a.u.) to 6.0 eV (2.190251E-01 a.u.).
• The top of the valence bands is at –5.1 eV (BAND 66 EMAX = -1.888888E-01 a.u.).

Electronic band structure

Electronic band structure is typically visualized as a band structure plot, where x-axis shows the wavevector k and y-axis shows the energy of each band at each k.

Reciprocal space path

To obtain a band structure plot, we need to define a path in the reciprocal space (inside the first Brillouin zone).

This is the reciprocal space path information SeeK-path yields if you upload a CIF of Cu₂O (section Reciprocal space and Brillouin-zone information):

Suggested path
Γ—X—M—Γ—R—X|R—M

High-symmetry points (scaled units)
Label	  k1               k2               k3
Γ     0.0000000000     0.0000000000     0.0000000000
M     0.5000000000     0.5000000000     0.0000000000
R     0.5000000000     0.5000000000     0.5000000000
X     0.0000000000     0.5000000000     0.0000000000
X1    0.5000000000     0.0000000000     0.0000000000

Creating the input file

The band structure calculation will yield a number of files, so it is better to create a new subdirectory called prop, where we copy the wavefunction file  Cu2O_Pn-3m_PBE0_TZVP_opt.w:

mkdir prop
cp Cu2O_Pn-3m_PBE0_TZVP_opt.w prop/Cu2O_Pn-3m_PBE0_TZVP_band.w
cd prop

Create the following property calculation input file Cu2O_Pn-3m_PBE0_TZVP_band.d3 in the directory:

BAND
SEEKPATH cP G-X-M-G-R-X|R-M
6 2 600 57 74 1 0
0 0 0   0 1 0    G -> X 
0 1 0   1 1 0    X -> M 
1 1 0   0 0 0    M -> G
0 0 0   1 1 1    G -> R 
1 1 1   0 1 0    R -> X 
1 1 1   1 1 0    R -> M 
ENDPROP

Line 1: Keyword BAND requests band structure calculation.

Line 2: Title for the job. It is highly recommended to include here (1) the source of the reciprocal space path and (2) the actual path. In this case, the path is from SeeK-path, the path is for primitive cubic Bravais lattice (cP), and the actual path is Γ-X-M-Γ-R-X|R-M. The wavevector k = 0 is called the Γ point and in the input file it is denoted as G to avoid Greek symbols in the input file. G-X means path from Γ to X and X|R means that this segment is not plotted (this is to save space, because before this segment we plotted R-X, which is similar to X-R). You can visualize the reciprocal space paths in SeeK-path or with the KVEC tool of the Bilbao Crystallographic Server.

Line 3: Options for the band structure calculation:

• 6 → number of path segments in reciprocal space (max. 20)
• 2 → Shrinking factor. The coordinates of the k-vectors are given in relative to this shrinking factor. For example, when the shrinking factor is 2, the coordinate 0 1/2 0 is expressed as 0 1 0 in the actual band path input.
• 600 → total number of k points along the path. Larger number of points means smoother band structure, but also requires more computational resources. Typically, something like 100 points per segment already gives a very smooth band structure. Here we have 6 segments and use 600 points.
• 57 → The first band to include in the calculation. The EMIN of this band is 2.75 eV below the top of the valence bands.
• 74 → The last band to include in the calculation. 10 valence bands (57-66) and 8 empty (conduction) bands (67-74) are included in the calculation. In this case, we will want to illustrate the band gap and few bands around it. 
• 1 → Produce files for plotting
• 0 → Do not print eigenvalues (energies)

Lines 4-9: Coordinates of the k-vectors, given relative to the shrinking factor

• 0 1 0 = 0 1/2 0 (X) 
• 1 1 0 = 1/2 1/2 0 (M)
• 1 1 1 = 1/2 1/2 1/2 (R)

Submitting the job

After this, the input is ready and the job can be submitted:

jsub -np 8 crystal Cu2O_Pn-3m_PBE0_TZVP_band.d3

In this example, the job will complete in few seconds, but for larger systems with lower symmetry and more bands, the band structure calculations are much more time-consuming (even hours).

Visualizing the band structure with CRYSPLOT

After the job has finished, you can visualize Cu2O_Pn-3m_PBE0_TZVP_band.BAND (or BAND.DAT) with CRYSPLOT → Make a plot → Band structure plot. Plot settings can be tuned from the menu on the left. Here is one recipe to create a reasonable plot:

1. Font size: Large or Very large
2. Plot title: Hide title (better add it later when creating a final figure e.g. in Powerpoint)
3. Tick labels: Customize k-points labels → Give the k-vector labels for your path. In our example, the labels are Γ, X, M, Γ, R, X|R, M (note the special label X|R for the discontinuous step in the path).
4. Y-axis unit: eV
5. Shifting Y-axis: Here the default setting is fine, but it is good to realize that the default is to shift the band energies so that the Fermi level is at 0 eV. In the case of Cu₂O with a band gap, it means that the top of the valence bands is shifted to 0 eV.
6. Fermi Energy line: If there is a band gap, I recommend hiding it. If the system is a metal, show the Fermi energy line.
7. Axis range: Do not change x-axis range. But y-axis range practically always has to be fine-tuned. Note that the range may still be written in Hartrees at this point, but you need to give it in eV. In the case of Cu₂O, Use a range of –2 to 4 eV. You should set the range in such way that all bands within that energy range have been calculated. For example, if you would set the range as  –4 to 4 eV, the band structure would look empty below –2.75 eV because the EMIN of band 57 was 2.75 eV below the top of the valence bands. If you would like to show valence bands down to –4 eV, you would need to redo the band structure calculation starting from band 47 (EMIN 4.6 eV below the top of the valence bands).
8. Grid: Without grid looks better.
9. Background color: White is usually better for further editing.

The resulting plot looks like this:

For use in presentations or publications, I would (1) change the y-axis label to "Energy (eV)", (2) change the x-axis label to "Wavevector k" (or leave it our completely), and (3) make the vertical lines and plot border more visible (black).

Visualizing the band structure with gen_bands_doss.py

It is also possible to visualize the band structure with the gen_bands_doss.py script by Marco Lorenz (2014). You need to have Python 3 available to run it. On Aalto CMAT computing clusters, you can load Python 3.7 (or later) with:

module load python/3.7

After this, you can run the script for the file Cu2O_Pn-3m_PBE0_TZVP_band.f25 as follows:

gen_bands_doss.py Cu2O_Pn-3m_PBE0_TZVP_band.f25

The script works as follows:

• It generates a gnuplot script file Cu2O_Pn-3m_PBE0_TZVP_band.plot.
• It runs gnuplot Cu2O_Pn-3m_PBE0_TZVP_band.plot, which produces EPS file Cu2O_Pn-3m_PBE0_TZVP_band.eps.
• If epstopdf is available, it runs epstopdf Cu2O_Pn-3m_PBE0_TZVP_band.eps, which produces PDF file Cu2O_Pn-3m_PBE0_TZVP_band.pdf.

The great thing about gen_bands_doss.py is that the gnuplot script file is very easy to modify to produce a nice-looking plot. For example, a slightly modified gnuplot script file Cu2O_Pn-3m_PBE0_TZVP_band.plot produces the following plot when you run gnuplot Cu2O_Pn-3m_PBE0_TZVP_band.plot:

Density of states (DOS)

A Density of states plot has energy as the x- or y-axis and density of electronic states as the other axis. We will carry out the density of states calculation in the prop directory created in the previous step. First, make another copy of the wavefunction:

cp Cu2O_Pn-3m_PBE0_TZVP_band.w Cu2O_Pn-3m_PBE0_TZVP_dos.w

Next, create the following property calculation input file Cu2O_Pn-3m_PBE0_TZVP_dos.d3 in the directory:

NEWK
0 32
32 32 32
1 0
DOSS
2 800 57 74 1 12 0
-4
1 2 3 4
-2
5 6
ENDPROP

Line 1: NEWK requests CRYSTAL to calculate the Fock/Kohn-Sham eigenvectors at a number of k-points in the reciprocal space. The lines 0 32 and 32 32 32 result in a four times denser k-mesh compared to the original 8×8×8 mesh used in the geometry optimization. DOS calculations in general require dense k sampling to produce accurate DOS. Line 1 0 requests CRYSTAL to calculate the Fermi energy using this k-mesh (1) and tells that no additional printing is required (0).

SHRINK
0 16
8 4 4

NEWK
0 32
32 16 16

Line 5: Keyword DOSS requests density of states calculation.

Line 6: Options for the DOSS calculation:

• 2 → DOS projections in addition to total DOS (can be 0 if you don't want any projections). The projections can show the DOS on particular atomic orbitals (for example, 3d orbitals of Cu) or on all atomic orbitals of particular atoms. We will calculate two projections that show DOS arising from Cu atoms and O atoms.
• 800 → Number of uniformly spaced energy values where DOS is calculated, from the bottom of the first band to the top of the last band. For smooth DOS plots, you should have hundreds of points (but probably not more than 1000).
• 57 → the first band (the same as for the band structure above)
• 74 → the last band (the same as for the band structure above)
• 1 → Produce files for plotting
• 12 → number of Legendre polynomials used to expand DOSS. CRYSTAL calculates DOS using a so-called Fourier-Legendre technique. Values between 10 and 18 are recommended in the manual. 12 is typically a good choice. If you see some extreme oscillations and significant negative DOS values, you may need to adjust the number of polynomials, increase the k-mesh, or increase the number of values on the energy axis (the second option on this row).
• 0 → No additional printing options

Line 7: First projected DOS: include all atomic orbitals of 4 atoms (negative sign means all orbitals of atoms listed on the next row).

Line 8: Atoms 1, 2, 3, and 4 (Cu). The numbers of the atoms refer to the contents of the primitive cell, as listed in the output file Cu2O_Pn-3m_PBE0_TZVP_opt.o.

Line 9: Second projected DOS: include all atomic orbitals of 2 atoms.

Line 10: Atoms 5 and 6 (O)

Submitting the job

After this, the input is ready and the job can be submitted:

jsub -np 8 crystal Cu2O_Pn-3m_PBE0_TZVP_dos.d3

In this example, the job will complete in few seconds, but for larger systems with lower symmetry and more bands, the density of states calculations are much more time-consuming (even hours).

Visualizing the DOS plot with CRYSPLOT

After the job has finished, you can visualize Cu2O_Pn-3m_PBE0_TZVP_dos.DOSS with CRYSPLOT → Make a plot → Density of states plot. This file is also rather straightforward to plot in any program that is able to plot XY data.

Usually, it is convenient to make a combined plot of band structure and DOS, The next to chapters deal with such combined plots.

Creating a combined band structure + DOS plot with CRYSPLOT

After you have both band structure and DOS available, you can create combined plots of them. One possibility is to use CRYSPLOT → Make a plot → Unified plot of band structure and density of states. You will need the files Cu2O_Pn-3m_PBE0_TZVP_band.BAND and Cu2O_Pn-3m_PBE0_TZVP_dos.DOSS. By following the recipe above for band structure plot, you should be able to get a plot similar to this:

Labels of DOS projections are missing, but from the largest to smallest they are as follows for the valence bands.: Total DOS > Cu > O.

Creating a combined band structure + DOS plot with gen_bands_doss.py

It is also possible to create a combined band structure + DOS plot with the gen_bands_doss.py script discussed above. You can run the script for files  Cu2O_Pn-3m_PBE0_TZVP_band.f25 and  Cu2O_Pn-3m_PBE0_TZVP_dos.f25 as follows:

gen_bands_doss.py Cu2O_Pn-3m_PBE0_TZVP_band.f25 Cu2O_Pn-3m_PBE0_TZVP_dos.f25

The script works as follows:

• It generates a gnuplot script file Cu2O_Pn-3m_PBE0_TZVP_band_multi.plot.
• It runs gnuplot Cu2O_Pn-3m_PBE0_TZVP_band_multi.plot, which produces EPS file Cu2O_Pn-3m_PBE0_TZVP_band_multi.eps.
• If epstopdf is available, it runs epstopdf Cu2O_Pn-3m_PBE0_TZVP_band_multi.eps , which produces Cu2O_Pn-3m_PBE0_TZVP_band_multi.pdf. 

It is highly recommended to modify the gnuplot script file to adjust the plot. A rather heavily modified gnuplot script file Cu2O_Pn-3m_PBE0_TZVP_band_multi.plot produces the following plot when you run gnuplot Cu2O_Pn-3m_PBE0_TZVP_band_multi.plot: