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.
atom and atom_frac
One can only use either the atom or the atom_frac keyword in a geometry.in file. Mixing those keyword in the same geometry.in file is not supported.
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 : 3B07D227-7545-412C-A1A0-2C2DBAF6241C
#
lattice_vector -0.00000003 2.73831639 2.73831639
lattice_vector 2.73831643 0.00000001 2.73831641
lattice_vector 2.73831643 2.73831641 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. We use the geometry obtained from the 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 two cores it takes about 5 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).
By examining the output, we find that the system only gained 3 meV during the geometry relaxation with the HSE06 functional. Overall, the final relaxed geometry is pretty close to the PBE result.
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 overal 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 bands 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), we find that -20 eV to 10 eV is a reasonable region around the band gap, so we can request:
output dos_tetrahedron -20.0 10.0 15001
dos_tetrahedron (this is also our recommendation for choice of integration method for smaller unit cells). The energy interval is determined by the next two float numbers, and the last integer specifies the number of points for this interval.
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
exx_band_structure_version 1
relativistic atomic_zora scalar
k_grid 12 12 12
output dos_tetrahedron -20 10 15001
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.
- For hybrid functionals, it is mandatory to specify the version of the band structure calculation by
exx_band_structure_version. Version 1 uses a faster real-space implementation; note that this keyword must be accompanied by a sufficiently densek_gridchoice. (Artifacts from an underconverged lattice Fourier transform might result for significantly underconvergedk_gridsettings.) Please refer to the manual to find more details regarding this keyword.
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_tetrahedron.dat
KS_DOS_total_tetrahedron.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_tetrahedron.dat: The raw total DOS data.KS_DOS_total_tetrahedron.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
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 GIMS result looks like this:
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.2082 eV for the present calculation.
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/Part-3/solutions/GaAs. At the end of these steps, we have an optimized structure of GaAs, which can be used to study SOC effects.
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
exx_band_structure_version 1
relativistic atomic_zora scalar
k_grid 8 8 8
output dos_tetrahedron -20 10 15001
output species_proj_dos_tetrahedron -20 10 15001
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_tetrahedron.dat
As_l_proj_dos_tetrahedron.dat.no_soc
As_l_proj_dos_tetrahedron_raw.dat
As_l_proj_dos_tetrahedron_raw.dat.no_soc
Ga_l_proj_dos_tetrahedron.dat
Ga_l_proj_dos_tetrahedron.dat.no_soc
Ga_l_proj_dos_tetrahedron_raw.dat
Ga_l_proj_dos_tetrahedron_raw.dat.no_soc
KS_DOS_total_raw_tetrahedron.dat
KS_DOS_total_raw_tetrahedron.dat.no_soc
KS_DOS_total_tetrahedron.dat
KS_DOS_total_tetrahedron.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_dos_tetrahedroncontain 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:
In this figure, the blue lines correspond to the band structure without SOC effects and the orange lines include the SOC correction. 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 right panel, one can find 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.
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
exx_band_structure_version 1
relativistic atomic_zora scalar
k_grid 8 8 8
output dos_tetrahedron -20 10 15001
output species_proj_dos_tetrahedron -20 10 15001
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:
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.