Part 3: Periodic Systems
This tutorial introduces the basic concepts of DFT calculations for periodic system geometries with FHI-aims. In a first step, we will carry out a structure optimization of a Si bulk crystal in the diamond structure. We will then consider spin-orbit coupling (SOC) effects in the energy band structure of GaAs.
Set up the geometry.in file for Si in the diamond structure
The syntax for the geometry.in file remains the same as previous two parts, with only a few additions:
lattice_vector 0.00000 2.71500 2.71500
lattice_vector 2.71500 0.00000 2.71500
lattice_vector 2.71500 2.71500 0.00000
atom_frac 0.00000 0.00000 0.00000 Si
atom_frac 0.25000 0.25000 0.25000 Si
The first three lines starting with lattice_vector define the three unit cell vectors \(a_i\) (\(i\)=1,2,3). For our example, we use the primitive unit cell of Si in the diamond structure. The following lines define the position and the species of the atoms for this system, however, this time in so-called fractional coordinates (atom_frac) \(f_1\), \(f_2\), \(f_3\). The Cartesian coordinates \(x\), \(y\), \(z\) of each atom follow by multiplying fractional coordinates and lattice vectors: \((x,y,z)= f_1 * a_1 + f_2 * a_2 + f_3 * a_3\).
Fractional coordinates are especially convenient when working with high-symmetry crystals. Nevertheless, one could also simply use Cartesian coordinates (in units of Å) by using the keyword atom at the beginning of each line.
unit cell choice
You may know that silicon crystallizes in the diamond structure, which is a face-centered cubic unit cell with eight atoms in a cube. It would appear that the unit cell shown above, however, does not define a cube (the unit cell vectors are not at 90 degree angles) and it only contains two atoms. What happened? In fact, the eight atom unit cell often shown in textbooks is not the smallest possible choice of unit cell that one can use to define the Si crystal in the diamond structure. The eight-atom cell is known as the "conventional unit cell", but the above choice of lattice vectors and coordinates leads to the exact same lattice of the crystal - only using two atoms, not eight. This is the meaning of the "primitive unit cell" term (used above), which is computationally more efficient.
atom and atom_frac
One can use the atom or the atom_frac keywords simultaneously in a geometry.in file. However, some specific functionality, sometimes external to FHI-aims (e.g., creating input files using the ASE framework), may not support mixing them - keep this in mind.
Visualizing Periodic Structures
A visual verification of the structure is usually key to avoid simple mistakes. This is even more important for periodic systems as nearest neighbors and atom distances cannot be identified easily. There are many tools to visualize crystal structures, e.g jmol, VESTA, and Avogadro to just name a few free software packages.
Here we will use GIMS, in particular its Structure Builder. One can upload the geometry file created from the above example and inspect the structure. Create a supercell to ensure that nearest neighbors and distances across the unit cell boundary look correct. Furthermore, GIMS provides additional information about the structure, such as lattice parameters, Bravais lattice type, space group number, and occupied Wyckoff positions. Finally, one can also create a snapshot of the structure and download it.
Set up the control.in file for Periodic Structures Relaxation
The control.in for structure optimization of a periodic system looks similar to the non-periodic case. First, we specify the runtime choices, which should look like this:
xc pbe
relativistic atomic_zora scalar
k_grid 8 8 8
relax_geometry bfgs 5e-3
relax_unit_cell full
k_grid, which is followed by three integer numbers. These three numbers specify the number of k-points along each reciprocal lattice vector. Their order corresponds to the order of lattice vectors specified in the geometry.in file. FHI-aims will create a three-dimensional, even-spaced, \(\Gamma\)-centered \(k\)-grid using the number of points specified for each reciprocal space direction.
Alternatively, the keyword k_grid_density followed by a single floating point number can be used. FHI-aims will create again an even-spaced, \(\Gamma\)-centered k-grid, with the same density along all reciprocal lattice vectors.
There is no absolute recommendation that one can make regarding whether to use k_grid or k_grid_density.
- In general, the keyword
k_grid_densityis as a practical option if one employs a sufficiently large k-point density. To obtain well converged calculations (with respect to the k-grid), it is often preferable to fulfill the criteria \(n_i * a_i > 50 Å\), for the lattice vector length \(a_i\) and the number of k-points \(n_i\) along the corresponding k-space direction \(i\). Usingk_grid_densityprevents the accidental use of over-converged settings, e.g., for larger supercell or slab calculations, where usually a much smaller number of k-points along the longer lattice vectors is sufficient. - To guarantee consistent k-grid numbers between different systems, e.g., when comparing total energies between different supercell sizes, it is preferable to use the
k_gridkeyword, which allows one to exactly control the number of k-points along each k-space direction. For example, if one doubles the unit cell size in each direction, you can divide the number of k-points in all directions by a factor of 2.
The necessary k-space grid can vary based on the target physical observable; for example, optical properties or densities of states generally require denser k_grid settings than just the total energy.
The 8x8x8 k-grid used in our current example is chosen based on our previous experience with Si bulk crystal. However, as just mentioned above, for any practical applications, one must ensure the convergence of physical properties of interest with respect to the k_grid setting. For example, a test using k_grid 12 12 12 or k_grid 16 16 16 would show how certain properties of interest, such as the total energy or the density of states (see below) reach numerical convergence as the resolution of the k-space grid is increased.
As before, the line starting with relax_geometry in the above control.in example requests a structure relaxation. However, unlike in non-periodic calculations, in case of periodic systems we do have six more degrees of freedom besides atomic positions: the lengths a, b, c of the lattice vectors \(a_1\), \(a_2\), \(a_3\), as well as the angles \(\alpha\), \(\beta\), and \(\gamma\) between the lattice vectors. These lattice degrees of freedom should be taken into account for structure relaxation.
Lattice vectors optimization in FHI-aims is specified by the keyword relax_unit_cell. There are several options for this keyword, but here we just want to freely optimize the lattice vectors, thus, setting relax_unit_cell full. If relax_unit_cell full is not set, then the relax_geometry keyword will keep the lattice fixed constant and only the atomic positions are optimized.
Before running the relaxation, do not forget to attach the species defaults for Si (light settings). Once the input files are ready, execute the calculation.
Optimization of lattice vectors: Analytical stress
Evaluating the force and stress components requires more time and memory compared to, e.g., total energy calculations. These components are therefore only evaluated once for each geometry, after the self-consistent field cycle is considered converged for the total energy.
Details about the implementation of the stress tensor in FHI-aims can be found in: Knuth, Franz, et al. "All-electron formalism for total energy strain derivatives and stress tensor components for numeric atom-centered orbitals." Comput. Phys. Commun. 190 (2015). DOI: https://doi.org/10.1016/j.cpc.2015.01.003
Quick look at the Result
As for the molecular test cases, the Output Analyzer of GIMS is the preferred method to visualize the structure optimization output data from FHI-aims.
Alternatively, at the command line, a useful perl script get_relaxation_info.pl in the [FHI-aims directory]/utilities folder can be used to print some important energies and residual forces along the relaxation trajectory from the main output file. In a terminal window, get_relaxation_info.pl is used as follows (please adapt [FHI-aims directory]):
[FHI-aims directory]/utilities/get_relaxation_info.pl aims.out
# Step Total energy [eV] E-E(1) [meV] Free energy [eV] F-F(1) [meV] max. force [eV/AA]
1 -15802.64745204 0.000000 -15802.64745204 0.000000 0.111748
2 -15802.65124091 -3.788873 -15802.65124091 -3.788873 0.076222
3 -15802.65421196 -6.759923 -15802.65421196 -6.759923 0.003273 converged.
geometry.in.next_step, we can see that the unit cell volume changed from 40.03 Å\(^3\) to 41.07 Å\(^3\).
FHI-aims does not assume any symmetry by default for relaxation. In general, this is the safer choice since it lets FHI-aims explore the full potential energy surface of the system during the relaxation. However, if one wishes to analyze structural and electronic properties of a system of known or constrained symmetry, it is helpful to restrict the degrees of freedom to only the parameters allowed by that symmetry. For example, the final result from our first relaxation only reveals very tiny deviations from the expected high-symmetry geometry:
#
# This is the geometry file that corresponds to the current relaxation step.
# If you do not want this file to be written, set the "write_restart_geometry" flag to .false.
# aims_uuid : F25E040D-A04F-46CE-AB4B-F9C455337798
#
lattice_vector -0.00000003 2.73831638 2.73831638
lattice_vector 2.73831641 0.00000001 2.73831640
lattice_vector 2.73831641 2.73831640 0.00000001
atom_frac 0.00000000 0.00000000 -0.00000000 Si
atom_frac 0.25000000 0.25000000 0.25000000 Si
#
# What follows is the current estimated Hessian matrix constructed by the BFGS algorithm.
# This is NOT the true Hessian matrix of the system.
# If you do not want this information here, switch it off using the "hessian_to_restart_geometry" keyword.
#
trust_radius 0.2000000030
hessian_file
Si is the actual global minimum.
In FHI-aims, there are different ways to incorporate symmetry in your relaxation. In particular, the symmetry related keywords rlsy_symmetry and use_symmetric_forces might be of interest for highly-symmetric systems or in order to constrain the symmetry of a solid in a particular way. To learn more about symmetry for relaxation of solids, please follow this tutorial.
Followup relaxation of solids with hybrid density functionals
Similar to relaxation of molecules, we perform a followup relaxation using the HSE06 XC functional (a hybrid functional), however, only using light species defaults for this example (and for the sake of this tutorial). For better converged results, we recommend intermediate settings (the results of which you can find in the directory Tutorial/3-Periodic-Systems/solutions/Si/HSE06_followup_relaxation_intermediate). We use the geometry obtained from the light relaxation above, which can be found in geometry.in.next_step. The control.in file should look like this:
xc hse06 0.11
hse_unit bohr
relativistic atomic_zora scalar
k_grid 8 8 8
relax_geometry bfgs 5e-3
relax_unit_cell full
[attach light species defaults for Si]
Follow the steps from the first part of the tutorial about restarting relaxations. It is identical for solids.
If you run the calculation only on few cores, you have some time to get a coffee (with four cores it takes about 15 minutes). Hybrid functionals are at the next rung of the Jacob's ladder. The HSE06 functional is physically more accurate than semilocal functionals for semi-conductors and insulators. However, it requires larger computational resources (time and memory).
After the calculation is finished, run the script get_relaxation_info.pl again. The output will look similar to the following example:
# Step Total energy [eV] E-E(1) [meV] Free energy [eV] F-F(1) [meV] max. force [eV/AA]
1 -15804.82402894 0.000000 -15804.82402894 0.000000 0.065454
2 -15804.82752017 -3.491229 -15804.82752017 -3.491229 0.000646 converged.
Similar to PBE functional, symmetry might be used to enhance the relaxation of highly-symmetric systems in case of HSE06 functional. For more information regarding symmetry in FHI-aims please follow this tutorial.
Post-Processing: Band structure and density of states
Looking into the output file, we already find some helpful information about the electronic structure such HOMO, LUMO, and estimated overall HOMO-LUMO gap energies. However, to get all the fine details regarding the electronic structure of a system, we have to request further output in our FHI-aims calculation: the band structure and the density of states.
FHI-aims requires explicit definition of the band path segment (a band path segment is a line segment in the k-space, which usually connects high symmetry points). This segment is determined by its start and end point in fractional k-space coordinates (first six floating numbers in control.in) and the number of points for the segment. Optionally, one can specify the name of the start and end point at the end of the corresponding line. The general syntax to request a band structure output looks like this:
output band k_xstart k_ystart k_zstart k_xend k_yend k_zend n_points symbol_start symbol_end
To request the output for the total DOS, we need to specify an energy interval (remember that FHI-aims is an all-electron code; outputting the whole KS-eigenvalue spectrum would create large and, in most cases, not very helpful files). Our former calculations help to identify the energy region of interest. Here, we are interested in the DOS around the gap. Looking into the output file from the last calculation (HSE06 relaxation; look for the last appearance of Writing Kohn-Sham eigenvalues. which corresponds to the converged structure):
Writing Kohn-Sham eigenvalues.
K-point: 1 at 0.000000 0.000000 0.000000 (in units of recip. lattice)
State Occupation Eigenvalue [Ha] Eigenvalue [eV]
1 2.00000 -66.589703 -1811.99801
2 2.00000 -66.589700 -1811.99794
3 2.00000 -5.378504 -146.35653
4 2.00000 -5.378177 -146.34763
5 2.00000 -3.684800 -100.26852
6 2.00000 -3.684800 -100.26852
7 2.00000 -3.684800 -100.26852
8 2.00000 -3.684187 -100.25184
9 2.00000 -3.684187 -100.25184
10 2.00000 -3.684187 -100.25184
11 2.00000 -0.705218 -19.18995
12 2.00000 -0.219401 -5.97022
13 2.00000 -0.219401 -5.97022
14 2.00000 -0.219401 -5.97021
15 0.00000 -0.096810 -2.63433
16 0.00000 -0.096810 -2.63433
17 0.00000 -0.096810 -2.63432
18 0.00000 -0.064767 -1.76240
19 0.00000 0.109755 2.98658
20 0.00000 0.109755 2.98658
Here, we find that -20 eV to 10 eV is a reasonable region around the band gap, so we can request:
output dos -20.0 10.0 15001 0.1
output dos e_min e_max n_points broadening
dos with a broadening parameter of 0.1 eV. If the goal is to resolve the individual peaks in the DOS, you can also experiment with smaller values such as 0.05 eV. The energy interval is determined by the first two float numbers, and the last integer specifies the number of points for this interval.
Different methods for DOS
Here, the k-space integration for calculating the DOS is performed using the Gaussian broadening method, as indicated by dos. For smaller unit cells, we typically recommend the tetrahedron method using dos_tetrahedron, which can resolve the individual peaks of the DOS more accurately for higher k-point densities. Note that sparse k-grids can lead to artifacts on the DOS when the tetrahedron method is used, however. Further details on different DOS methods are covered below.
Setting up all these pieces manually is tedious and error prone. Therefore, we recommend tools that automate this step. The band structure workflow in GIMS provides such a tool which uses the Setyawan/Curtarolo conventions for band paths. To reproduce the control.in example given below, you need to use a k-space density of 30 k-points ⋅ Å for the band path. GIMS also lets you to define the DOS output.
For Si-diamond, your final control.in file should be similar to the following version:
xc hse06 0.11
hse_unit bohr
hybrid_xc_coeff 0.25
relativistic atomic_zora scalar
k_grid 12 12 12
output dos -20 20 15001 0.1
output band 0.00000 0.00000 0.00000 0.50000 0.00000 0.50000 35 G X
output band 0.50000 0.00000 0.50000 0.50000 0.25000 0.75000 17 X W
output band 0.50000 0.25000 0.75000 0.37500 0.37500 0.75000 12 W K
output band 0.37500 0.37500 0.75000 0.00000 0.00000 0.00000 37 K G
output band 0.00000 0.00000 0.00000 0.50000 0.50000 0.50000 30 G L
output band 0.50000 0.50000 0.50000 0.62500 0.25000 0.62500 21 L U
output band 0.62500 0.25000 0.62500 0.50000 0.25000 0.75000 12 U W
output band 0.50000 0.25000 0.75000 0.50000 0.50000 0.50000 24 W L
output band 0.50000 0.50000 0.50000 0.37500 0.37500 0.75000 21 L K
output band 0.62500 0.25000 0.62500 0.50000 0.00000 0.50000 12 U X
[attach light species defaults for Si]
Noteworthy features of this control.in file include:
- The k-point grid numbers are increased to obtain better converged results for the band structure and density of states. This normally increases the computation time but here the calculation will still run reasonably fast, as we do not request the computation of forces and stresses.
To start the calculation, we use the relaxed structure of previous step which can be found in geometry.in.next_step. Copy this file into the geometry.in of you current simulation and execute the run. After the calculation has finished, you should find the following file in the simulation directory:
KS_DOS_total_raw.dat
KS_DOS_total.dat
aims.out
band1001.out
band1002.out
band1003.out
band1004.out
band1005.out
band1006.out
band1007.out
band1008.out
band1009.out
band1010.out
control.in
geometry.in
The new files carry the following information:
KS_DOS_total_raw.dat: The raw total DOS data.KS_DOS_total.dat: The total DOS shifted by the chemical energy of the electrons (Fermi level).band****.out: The band structure data for each line segment in the same order as specified in thecontrol.infile. The four digits followingbandhave the following meaning:- The first digit specifies the spin channel (here, we have only one spin channel, thus, all file names start with
1). - The next three digits specifies the band segment number (as ordered in the
control.in). For this calculation we requested ten segments, hence the files are numbered from 1 to 10.
The next step in a post-processing calculation is to plot the data. We facilitate this process by providing a script in the FHI-aims utility suite clims, which can be called with clims-aimsplot command. Further information regarding how to install clims is accessible in the preparation section of this tutorial. The script can be called from the directory including all of your output files simply by:
clims-aimsplot --band --total_dos
clims-aimsplot --help
clims-aimsplot will not save an image file by default, and one must either use the controls in the pop-up, or include the --save_only option, which will not produce a pop-up, but instead save a file aimsplot.png. The final command used to produce the plot below, including a narrowing of the plot energy window from -5 eV to 5 eV is
clims-aimsplot --band --total_dos --emin -5 --emax 5 --save_only
It is also possible to plot the band structure and DOS using the GIMS output analyzer. Just select all files from your calculation directory and GIMS sorts out the needed data.
The band gap can be extracted from the output files of the band structure calculation. In particular, for indirect gap materials such as present case, where the conduction band minimum and valence band maximum does not occur at the same k-point, it is more accurate to estimate the band gap by examining the band structure. The reason for this is that these two k-points may not be present in the original k-point grid (specified by k_grid) used for the simulation. We find a band gap of 1.206 eV for the present calculation.
We note that the DOS here is not very well resolved, and contains a number of peaks. To obtain a better resolution, smoother DOS, we can utilise the keyword dos_kgrid_factors. This performs a Fourier interpolation of the k-grid used in the SCF cycle to yield a denser k-grid without self-consistency, upon which the DOS is calculated, in what is sometimes called a perturbative DOS calculation. We can, for example, include in our control.in file
dos_kgrid_factors 3 3 3
This is, however, a much more computationally expensive calculation, and we do not recommend running it locally on your PC. The result may be found in the solutions folder, 3-Periodic-Systems/solutions/Si/HSE06_dos_kgrid_factors. Further details on computing the DOS, including the tetrahedron method are presented in the following section.
Spin-orbit coupling
Relativistic effects, such as spin-orbit coupling (SOC), become more important as the elements forming the material become heavier. GaAs is a well-known example that shows the impact of spin-orbit coupling on the electronic properties; therefore we use it as a reference example for introducing spin-orbit coupling here.
The first part of this section will also serve as a recap of our standard workflow for structure relaxations. Perform the following steps for GaAs:
- Find the structure and initial guess of the structure for GaAs.
- Perform a PBE relaxation with light species defaults.
- Perform a HSE06 followup relaxation with light (or, better, intermediate) species defaults.
You can find the solution for the above three steps in the folder Tutorial/3-Periodic-Systems/solutions/GaAs. At the end of these steps, we have an optimized structure of GaAs, which can be used to study SOC effects. Note that the HSE06 followup relaxation with light species defaults will take approximately 30 minutes on 4 cores.
At the time of writing, FHI-aims offers a second-variational SOC correction to the self-consistent scalar-relativistic eigenstates. More details regarding this approach are provided in the following paper:
William P. Huhn and Volker Blum, "One-hundred-three compound band-structure benchmark of post-self-consistent spin-orbit coupling treatments in density functional theory", Phys. Rev. Materials 1, 033803 (2017).
To activate the SOC correction, the include_spin_orbit keyword should be added to the control.in file. Since Si and GaAs have the same Bravais lattice, we can also use the same band path through the Brillouin zone as in the previous section. In addition, we will also request a species-projected DOS to better understand the composition of the band edges.
xc hse06 0.11
hse_unit bohr
hybrid_xc_coeff 0.25
relativistic atomic_zora scalar
k_grid 8 8 8
output dos -20 10 15001 0.1
output species_proj_dos -20 10 15001 0.1
output band 0.00000 0.00000 0.00000 0.50000 0.00000 0.50000 33 G X
output band 0.50000 0.00000 0.50000 0.50000 0.25000 0.75000 17 X W
output band 0.50000 0.25000 0.75000 0.37500 0.37500 0.75000 12 W K
output band 0.37500 0.37500 0.75000 0.00000 0.00000 0.00000 35 K G
output band 0.00000 0.00000 0.00000 0.50000 0.50000 0.50000 29 G L
output band 0.50000 0.50000 0.50000 0.62500 0.25000 0.62500 20 L U
output band 0.62500 0.25000 0.62500 0.50000 0.25000 0.75000 12 U W
output band 0.50000 0.25000 0.75000 0.50000 0.50000 0.50000 23 W L
output band 0.50000 0.50000 0.50000 0.37500 0.37500 0.75000 20 L K
output band 0.62500 0.25000 0.62500 0.50000 0.00000 0.50000 12 U X
include_spin_orbit
[attach light species defaults for Ga and As]
Start the calculation using the above control.in (don't forget to attach the "light" or, better, "intermediate" species defaults) and the geometry that was optimized using the HSE06 functional.
After the calculation has finished, you should see the following files in your folder:
As_l_proj_dos.dat
As_l_proj_dos.dat.no_soc
As_l_proj_dos_raw.dat
As_l_proj_dos_raw.dat.no_soc
Ga_l_proj_dos.dat
Ga_l_proj_dos.dat.no_soc
Ga_l_proj_dos_raw.dat
Ga_l_proj_dos_raw.dat.no_soc
KS_DOS_total_raw.dat
KS_DOS_total_raw.dat.no_soc
KS_DOS_total.dat
KS_DOS_total.dat.no_soc
aims.out
band1001.out
band1001.out.no_soc
band1002.out
band1002.out.no_soc
band1003.out
band1003.out.no_soc
band1004.out
band1004.out.no_soc
band1005.out
band1005.out.no_soc
band1006.out
band1006.out.no_soc
band1007.out
band1007.out.no_soc
band1008.out
band1008.out.no_soc
band1009.out
band1009.out.no_soc
band1010.out
band1010.out.no_soc
control.in
geometry.in
- All files (for both band and DOS outputs) ending with
*.datcontain the data with the SOC correction. - All files (for both band and DOS outputs) ending with
*.dat.no_soccontain the data without SOC correction, only for reference purposes. - For DOS outputs, the files including
*_l_proj_doscontain the species-projected DOS.
To visualize the band structure and DOS files, you can either execute clims-aimsplot in the results directory or upload it to GIMS' output analyzer. The visualized data should look like following:
where the command to produce the top figure was
clims-aimsplot --band --emin -5 --emax 5
clims-aimsplot --band --species_dos --total_dos --emin -5 --emax 5
In the top figure, the black lines correspond to the band structure including SOC effects and the red lines without the SOC correction (SR=scalar relativistic). To spot the effect of SOC, look at the valence band maximum (VBM) at the \(\Gamma\)-point. Before the spin-orbit correction, the states at the VB were degenerate. The SOC correction splits the bands at the VBM according to their symmetry.
In the bottom figure, we visualize the band structure and the total (black) and the species projected DOS of As (purple) and Ga (brown). The valence bands are predominantly contributed by the As atoms, while the conduction bands have contributions from both As and Ga atoms. Note that the conduction band minimum (CBM) is not yet well captured in the DOS (the onset of the DOS is, in fact, at the CBM, but a more highly resolved plot would show this). In fact, the k_grid 8 8 8 used in our example, while computationally manageable, is not yet sufficient to completely resolve this important part of the DOS. In GaAs the CBM is very steep and consequently a much denser k-grid would be needed to capture its correct value. In the following section, we briefly outline the methods available to better converge the DOS.
DOS Convergence
In FHI-aims, there are three approaches one can use to obtain a better resolved DOS:
-
Use a larger
k_gridin the SCF cycle. -
Use the tetrahedron method for the DOS (keyword
output dos_tetrahedron). -
Fourier interpolation of the k-grid to yield a denser k-grid without self-consistency for the DOS (keyword
dos_kgrid_factors). As mentioned in Si section, this is sometimes referred to as a perturbative DOS calculation.
Option 1 is the obvious one, however it can rapidly become more computationally expensive, and in the following we will explore options 2 and 3 in more detail, with GaAs as the example. In particular, we will present the following 4 examples:
- Perturbative Gaussian DOS with a 8x8x8 SCF k-grid
- Perturbative Gaussian DOS with a 12x12x12 SCF k-grid
- Tetrahedron DOS with a 12x12x12 SCF k-grid
- Perturbative tetrahedron DOS with a 12x12x12 SCF k-grid
These calculations would take a while to run on a personal computer, therefore for convenience, the results following submission to a HPC can be found in the directory 3-Periodic-Systems/solutions/GaAs_DOS. Note that these calculations (with intermediate, rather than light species defaults, which we would recommend for publication) also appear in the /testcases/GaAs_HSE06+SOC/ folder in your FHI-aims installation.
1. Perturbative Gaussian DOS with 8x8x8 SCF k-grid
In order to switch on the perturbative approach to compute the DOS on a denser post-SCF k-grid we use the keyword
dos_kgrid_factors 4 4 4
output dos keyword. This produces an interpolated k-grid that is denser by a factor of 4 in each direction.
This adds computational expense, however is much cheaper than calculating a 32x32x32 SCF
k_grid. The DOS from this calculation is shown below, compared to that without dos_kgrid_factors
2. Perturbative Gaussian DOS with 12x12x12 SCF k-grid
Next, let us observe the effect of using a larger SCF k-grid (k_grid 12 12 12), as well as changing the broadening parameter for the Gaussian DOS from 0.1 eV, to 0.05 eV. The choice of broadening value affects the DOS; too large, and you risk missing features. For example around the conduction band minimum, the gap to the valence band maximum is narrowed by choosing a larger broadening value. This is especially true for materials that have a sharp CBM, such as GaAs, as displayed in the plot below:
3. Tetrahedron DOS with 12x12x12 SCF k-grid
The tetrahedron method for the DOS can be requested in FHI-aims via the keyword
in control.in
output dos_tetrahedron -20.0 10.0 15001
dos_kgrid_factors
We can observe that the sharpness of the CBM is well reproduced, however there are a number of spikes in the plot, which are not present in the Gaussian method, but are well known to occur with the tetrahedron method.
4. Perturbative Tetrahedron DOS with 12x12x12 SCF k-grid
We can also combine the tetrahedron DOS approach with the perturbative DOS calculation. In this case we use
dos_kgrid_factors 2 2 2
output dos_tetrahedron and k_grid 12 12 12.
The DOS from this calculation is shown below, compared with the tetrahedron result without the perturbative DOS
Here we can observe that many of the spikes no longer appear in the plot.
To summarise, we recommend the use of the dos_kgrid_factors keyword, if affordable, with either the tetrahedron or Gaussian based DOS method to get a converged DOS, especially for materials that have a difficult to converge DOS like GaAs. Choosing a large enough dos_kgrid_factors, as well as a small enough broadening with the Gaussian method to not obscure features may be somewhat a trial and error process. Therefore we also recommend the use of the keyword
elsi_restart write scf_converged
elsi_restart read
Projected DOS calculations
Please note that as of version 240717 the atom and species projected DOS for the Gaussian method DOES include contributions from dos_kgrid_factors.
However, as of version 240920 the atom and species projected DOS for the tetrahedron method does NOT include contributions from dos_kgrid_factors, i.e. it is calculated on the SCF k-grid. This is intended to be included in
the near future.
Mulliken-projected band structure
To gain more insights into the electronic structure, it is also possible to output the Mulliken-projected band structure. The Mulliken projection is the projection of electronic states onto contributions from individual atoms, using the basis functions that are attached to each atom.
To activate the Mulliken-projected band structure in output of FHI-aims, one only needs to change the band keyword to band_mulliken, as demonstrated in the following snippet:
xc hse06 0.11
hse_unit bohr
hybrid_xc_coeff 0.25
relativistic atomic_zora scalar
k_grid 8 8 8
output dos -20 10 15001 0.1
output species_proj_dos -20 10 15001 0.1
output band_mulliken 0.00000 0.00000 0.00000 0.50000 0.00000 0.50000 33 G X
output band_mulliken 0.50000 0.00000 0.50000 0.50000 0.25000 0.75000 17 X W
output band_mulliken 0.50000 0.25000 0.75000 0.37500 0.37500 0.75000 12 W K
output band_mulliken 0.37500 0.37500 0.75000 0.00000 0.00000 0.00000 35 K G
output band_mulliken 0.00000 0.00000 0.00000 0.50000 0.50000 0.50000 29 G L
output band_mulliken 0.50000 0.50000 0.50000 0.62500 0.25000 0.62500 20 L U
output band_mulliken 0.62500 0.25000 0.62500 0.50000 0.25000 0.75000 12 U W
output band_mulliken 0.50000 0.25000 0.75000 0.50000 0.50000 0.50000 23 W L
output band_mulliken 0.50000 0.50000 0.50000 0.37500 0.37500 0.75000 20 L K
output band_mulliken 0.62500 0.25000 0.62500 0.50000 0.00000 0.50000 12 U X
include_spin_orbit
[attach light species defaults for Ga and As]
Run the calculation with the above settings (do not forget to attach the light species defaults) for the same converged geometry as before.
At the time of writing this tutorial, GIMS does not support visualizing the Mulliken-projected band structure. However, one can simply use clims-aimsplot for plotting. Your final plot should be similar to the following:
which was obtained using the command
clims-aimsplot --bandmlk --species_dos --total_dos --emin -5 --emax 5 --scale_mulliken_marker 2 --legend_offset 0.9 0.4
The darker the color of the dots, the more do they correspond to that species. Using the Mulliken-projected band structure, the contributions of each species to the band structure of different parts of the Brillouin zone becomes clearer.
Solutions
You find all the solution to all the above exercises by clicking on the button below.