In this example, we will plot the electronic band structure and density of states (DOS) of Cu_{2}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`

.

**Cu2O_Pn-3m_PBE0_TZVP_opt.d3**

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`

:

**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_{2}O (section *Reciprocal space and Brillouin-zone information*):

**Reciprocal space path from SeeK-path**

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:

**Cu2O_Pn-3m_PBE0_TZVP_band.d3**

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:

**Font size:**Large or Very large**Plot title:**Hide title (better add it later when creating a final figure*e.g.*in Powerpoint)**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).**Y-axis unit:**eV**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_{2}O with a band gap, it means that the top of the valence bands is shifted to 0 eV.**Fermi Energy line:**If there is a band gap, I recommend hiding it. If the system is a metal, show the Fermi energy line.**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_{2}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).**Grid**: Without grid looks better.**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:

**Cu2O_Pn-3m_PBE0_TZVP_dos.d3**

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).

Double k-mesh for Fermi energy

Note that if you were originally using a double mesh for the Fermi energy (for example, for metallic systems), you probably do not need to use a double mesh for the DOS calculation. In fact, it may lead in so many *k*-points that you run out of memory. For example, if the original SHRINK input was

SHRINK

0 16

8 4 4

the four times denser mesh for the DOS could be

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, 3
*d*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`

: