In the previous tutorial, we calculated the total energy of silicon and explored its Brillouin zone convergence. This tutorial illustrates calculations of the electronic band structure, specifically, the variation of the Kohn-Sham eigenvalues along a special kpoint path in the Brillouin zone. It will also introduce an alternate algorithm for converging the electronic state, the self-consistent field (SCF) method.
First lets specify the common specifications for bulk silicon (same as previous tutorial):
#Save the following to common.in: lattice face-centered Cubic 10.263 ion-species GBRV/$ID_pbe.uspp elec-cutoff 20 100 ion Si 0.00 0.00 0.00 0 ion Si 0.25 0.25 0.25 0
and a total energy calculation:
#Save the following to totalE.in: include common.in kpoint-folding 8 8 8 #Use a Brillouin zone mesh electronic-SCF #Perform a Self-Consistent Field optimization dump-name totalE.$VAR dump End State BandEigs #State and band eigenvalues dump End ElecDensity #Save the self-consistent electron density dump End EigStats #Get eigenvalue statistics
which we run with
mpirun -n 4 jdftx -i totalE.in | tee totalE.out
Note the new command electronic-scf, which uses a different algorithm from the variational minimize to solve the Kohn-Sham equations. Essentially, Kohn-Sham eigenvalues and orbitals are calculated for an input potential, and a new density and potential are constructed from those orbitals. This output potential would in general be different from the input one, so the above step is repeated with successive guesses for the input potential until the input and output potentials become identical (within some threshold). In practice, the guess for the input potential at one step is based on input and output potentials from several previous steps (Pulay algorithm).
Examine the output file: instead of ElecMinimize lines, you will find lines starting with SCF, which report the change in energy and eigenvalues between "cycles". This algorithm is usually faster than the variational minimize, but unlike the default minimize algorithm, it is not guaranteed to converge.
Note that whether we used SCF or Minimize is not important for the remainder of this tutorial demonstrating band structure calculations; we will interchangeably use both algorithms for total energy calculations from now on to gain familiarity with both.
Next, we list high-symmetry points in the Brillouin zone laying out a path along which we want the band structure:
#Save the following to bandstruct.kpoints.in kpoint 0.000 0.000 0.000 Gamma kpoint 0.000 0.500 0.500 X kpoint 0.250 0.750 0.500 W kpoint 0.500 0.500 0.500 L kpoint 0.000 0.000 0.000 Gamma kpoint 0.375 0.750 0.375 K
For common crystal structures, you can find the high-symmetry points easily on the web, eg. this course website. For a more complete listing, you can consult a crystallography database eg. the Bilbao database, but to use that, you will first need to find the space group of your crystal.
From the above path specification, we can generate a sequence of points along the path and a plot script by running bandstructKpoints (in the jdftx/scripts directory):
bandstructKpoints bandstruct.kpoints.in 0.05 bandstruct
This should generate a file bandstruct.kpoints containing kpoints along the high-symmetry path and a GNUPLOT script bandstruct.plot. The second parameter, dk, of bandstructKpoints specifies the typical distance between kpoints in (dimensionless) reciprocal space; decreasing dk will produce more points along the path, which will take longer to calculate, but produce a smoother plot. Note that the last column of banstruct.kpoints.in is a label for the special Brillouin zone point which is used to label the plot.
Now we can run a band structure calculation along this path with the input file:
#Save the following to bandstruct.in include common.in include bandstruct.kpoints #Get kpoints along high-symmetry path created above fix-electron-density totalE.$VAR #Fix the electron density (not self-consistent) elec-n-bands 10 #Number of bands to solve for dump End BandEigs #Output the band eigenvalues for plotting dump-name bandstruct.$VAR #This prefix should match the final paramater given to bandstructKpoints
which we run using:
mpirun -n 4 jdftx -i bandstruct.in | tee bandstruct.out
which produces bandstruct.eigenvals. Finally running
gnuplot --persist bandstruct.plot
will generate a plot of the electronic band structure.
Note that the y-axis (energy) is in Hartrees and its absolute value is not particularly meaningful because the electrostatic potential in a 3D periodic system does not have an unambiguous zero reference point. Typically electronic band structure plots are referenced to the valence band maximum (VBM or HOMO) energy at zero for insulators, or the Fermi level for metals. Look at the eigStats output from the total energy calculation (totalE.eigStats) to identify the VBM (HOMO) energy and replace the final line of the auto-generated bandstruct.plot with the following lines:
set xzeroaxis #Add dotted line at zero energy set ylabel "E - VBM [eV]" #Add y-axis label set yrange [*:10] #Truncate bands very far from VBM plot for [i=1:nCols] "bandstruct.eigenvals" binary format=formatString u 0:((column(i)-VBM)*27.21) w l
replacing VBM in the last line with the relevant value from totalE.eigStats. Note that this modification also converts the energy from Hartrees to eV for the plot. Now rerunning gnuplot produces the plot shown below and at the top of this page. Notice that at the Gamma point, the lowest band is single while the next three higher bands are degenerate: these line up with the s and p valence orbitals on the Silicon atoms. These degeneracies change in different parts of the Brillouin zone: the XW segment has two pairs of degenerate bands, while the WL and Gamma-K segments have no degeneracies.