Part 1: Optical Conductivity of Aluminum
In this example, we will calculate the optical conductivity \(\sigma(\omega)\) of a 4-atom aluminum cell (fcc structure) at an extremely high electronic temperature of \(k_BT=1\) eV and we will then compare the results with those in the literature. We will use two different k-grids for out Al cell, a \(6\times6\times6\) k-grid and an \(18\times18\times18\) k-grid. We will compute the \(18\times18\times18\) k-grid in two ways - either with, or without Fourier interpolation. Note that in this example, the nuclear motion is neglected.
1. 6x6x6 k-grid
First, we will compute the optical conductivity of Al on a \(6\times6\times6\) k-grid. Let us create a directory to perform this calculation in and go into it:
Like with any other calculation in FHI-aims, we require the input filesgeometry.in and control.in.
Geometry
An example geometry.in is provided below; it can be visualized with GIMS. Note that the cell is in equilibrium, i.e. no displacements are included.
lattice_vector 4.0389296899999998 0.0000000000000000 0.0000000000000000
lattice_vector 0.0000000000000000 4.0389296899999998 0.0000000000000000
lattice_vector 0.0000000000000000 0.0000000000000000 4.0389296899999998
atom 0.0000000000000000 0.0000000000000000 0.0000000000000000 Al
atom 0.0000000000000000 2.0194648449999999 2.0194648449999999 Al
atom 2.0194648449999999 0.0000000000000000 2.0194648449999999 Al
atom 2.0194648449999999 2.0194648449999999 0.0000000000000000 Al
Control tags
occupation_type fermi 1
compute_kubo_greenwood 0.06 1 -20.57 9.44 0.0 5.0 2000 a a
kg_width 0.01 0.10 10
dos_kgrid_factors 1 1 1
First, we specify that a Fermi distribution shall be used to determine the Fermi level and the occupation numbers at the desired electronic temperature of \(k_BT=1\) eV:
Here1 means that a Fermi function broadening of \(k_BT=1\) eV is applied during the SCF cycle. Please be aware that conductivity is directly related to the carrier density via the occupation numbers. Accordingly,
it is essential to use the same carrier densities when comparing results, e.g., to literature.
The tag to invoke the KG calculation is:
- Here
0.02specifies the broadening width of the broadening function used to represent the delta function in the KG formula.
Warning
This parameter is a legacy of the old KG code and is deprecated in the current KG implementation. Accordingly, it is present for backward compatibility and will be removed in the future.
Currently, the broadening widths of the delta function should be specified by the tag kg_width introduced below.
-
1indicates that a Fermi temperature of \(k_BT=1\) eV is applied in the evaluation of the KG formula. -
-20.57 9.44(in eV) means the energy window (related to the internal zero) within which the momentum-matrix elements are considered. -
0.0 5.0 2000means that the optical conductivity \(\sigma(\omega)\) is calculated between \(\hbar\omega_{min} = 0.0\) eV and \(\hbar\omega_{max} = 5.0\) eV with 2000 steps. -
a ameans that the output shall report the averaged conductivity over all three Cartesian coordinates \(\sigma_{ii}(\omega)\).
The kg_width tag:
0.01 0.10 10 indicates 10 evenly-spaced broadening widths from \(0.01eV\) to \(0.10eV\). The default broadening function is the Gaussian function.
Tag dos_kgrid_factors invokes a Fourier interpolation to produce a denser post-SCF k-grid, as also used for instance for density-of-state calculations. Here we simply set:
In order to obtain the complete control.in file, we also need to copy the "light" species default for Al:
[FHI-aims-directory] is the location of your FHI-aims directory. Webinar users should find this location specified in the README file in your home directory.
A Valid control.in for this example
xc pbe
relativistic atomic_zora scalar
k_grid 6 6 6
# type[letter] temperature[eV]
occupation_type fermi 1
# kubo_broadening[eV] , Fermi-Temperature[eV] , E_min[eV] , E_max[eV] , w_min[eV] , w_max[eV] , n_w_points[1] , spatial directions[letter]
compute_kubo_greenwood 0.02 1 -20.57 9.44 0.0 5.0 2000 a a
kg_width 0.01 0.10 10
dos_kgrid_factors 1 1 1
################################################################################
#
# FHI-aims code project
# VB, Fritz-Haber Institut, 2009
#
# Suggested "light" defaults for Al atom (to be pasted into control.in file)
# Be sure to double-check any results obtained with these settings for post-processing,
# e.g., with the "tight" defaults and larger basis sets.
#
# 2020/09/08 Added f function to "light" after reinspection of Delta test outcomes.
# This was done for all of Al-Cl and is a tricky decision since it makes
# "light" calculations measurably more expensive for these elements.
# Nevertheless, outcomes for P, S, Cl (and to some extent, Si) appear
# to justify this choice.
#
################################################################################
species Al
# global species definitions
nucleus 13
mass 26.9815386
#
l_hartree 4
#
cut_pot 3.5 1.5 1.0
basis_dep_cutoff 1e-4
#
radial_base 41 5.0
radial_multiplier 1
angular_grids specified
division 0.6594 110
division 0.8170 194
division 0.9059 302
# division 1.0363 434
# division 1.1443 590
# division 1.2621 770
# division 2.8177 974
# outer_grid 974
outer_grid 302
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 3 s 2.
valence 3 p 1.
# ion occupancy
ion_occ 3 s 1.
ion_occ 2 p 6.
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Constructed for dimers: 2.0 A, 2.5 A, 3.0 A, 3.75 A, 4.5 A
#
################################################################################
# "First tier" - improvements: -199.47 meV to -10.63 meV
ionic 3 d auto
ionic 3 p auto
hydro 4 f 4.7
ionic 3 s auto
# "Second tier" - improvements: -5.35 meV to -1.57 meV
# hydro 5 g 7
# hydro 3 d 6
# hydro 2 s 11.6
# hydro 2 p 0.9
# "Third tier" - improvements: -0.63 meV to -0.20 meV
# hydro 5 f 7.6
# hydro 4 p 7.2
# hydro 4 s 3.7
# hydro 4 d 7.6
# "Fourth tier" - improvements: -0.17 meV to -0.08 meV
# hydro 4 d 13.6
# hydro 5 g 11.2
# hydro 4 d 0.9
# hydro 1 s 0.4
# hydro 4 p 0.1
# hydro 5 f 9.8
# Further basis functions that fell out of the optimization - noise level...
# hydro 4 p 5
Run the calculation
If you are running the calculation on a compute cluster, for example in the webinar, please use a command like the following:
Example for the submission script on the webinar cluster
#!/bin/bash -l
# Standard output and error:
#SBATCH -o ./slurm.stdout
#SBATCH -e ./slurm.stderr
# Initial working directory:
#SBATCH -D ./
# Job Name:
#SBATCH -J FHI-aims
# Number of nodes and MPI tasks per node:
#SBATCH --nodes=1
# HPC7a
#SBATCH --tasks-per-node=192
#
# Wall clock limit:
#SBATCH --time=00:05:00
module purge
module load fhi-aims vibes
date
mpirun aims.x > aims.out
If running this calculation on a local machine, a parallel calculation with N processes can be launched by:
your_aims.x should be replaced by your FHI-aims binary file.
Checking the result
After the FHI-aims calculation has finished, there should be 20 output files named KG_full_L11_out_Gaussian_******_a_a.out and KG_full_Mobility_out_Gaussian_******_a_a.out. Here L11 denotes the respective Onsager coefficient,
i.e., the conductivity, and ****** denotes the Gaussian broadening value. An example output for a 0.06 eV Gaussian broadening looks like this:
# Gaussian broadening used, width 6.0000E-02 eV
# Number of electrons above Fermi-energy: 0.102689661033E+01 0.155857548643E+23 (ch.carr./cm^-3)
# Number of holes below Fermi-energy: 0.114726698070E+01 0.174126798601E+23 (ch.carr./cm^-3)
# Frequency Inter-band Full Intra-band Full Both-band Full Inter-band Hole Intra-band Hole Both-band Hole Inter-band Elec Intra-band Elec Both-band Elec
# eV (Ohm cm)^-1 (Ohm cm)^-1 (Ohm cm)^-1 (Ohm cm)^-1 (Ohm cm)^-1 (Ohm cm)^-1 (Ohm cm)^-1 (Ohm cm)^-1 (Ohm cm)^-1
2.500000E-003 3.141151E+002 2.605452E+005 2.608593E+005 4.975609E+001 1.086660E+005 1.087158E+005 2.643590E+002 1.518792E+005 1.521436E+005
5.000000E-003 1.694669E+002 2.598676E+005 2.600370E+005 2.936314E+001 1.083834E+005 1.084127E+005 1.401038E+002 1.514842E+005 1.516243E+005
7.500000E-003 1.219543E+002 2.587421E+005 2.588641E+005 2.310366E+001 1.079140E+005 1.079371E+005 9.885066E+001 1.508281E+005 1.509270E+005
1.000000E-002 9.878901E+001 2.571747E+005 2.572735E+005 2.044541E+001 1.072602E+005 1.072807E+005 7.834360E+001 1.499144E+005 1.499928E+005
1.250000E-002 8.542119E+001 2.551733E+005 2.552587E+005 1.929025E+001 1.064255E+005 1.064448E+005 6.613095E+001 1.487478E+005 1.488139E+005
1.500000E-002 7.700853E+001 2.527484E+005 2.528254E+005 1.894665E+001 1.054142E+005 1.054331E+005 5.806188E+001 1.473342E+005 1.473923E+005
1.750000E-002 7.148317E+001 2.499122E+005 2.499837E+005 1.912591E+001 1.042313E+005 1.042504E+005 5.235726E+001 1.456809E+005 1.457333E+005
2.000000E-002 6.781859E+001 2.466792E+005 2.467470E+005 1.969121E+001 1.028829E+005 1.029026E+005 4.812737E+001 1.437963E+005 1.438445E+005
2.250000E-002 6.545200E+001 2.430657E+005 2.431312E+005 2.057416E+001 1.013758E+005 1.013964E+005 4.487784E+001 1.416899E+005 1.417348E+005
2.500000E-002 6.405337E+001 2.390897E+005 2.391538E+005 2.174138E+001 9.971752E+004 9.973926E+004 4.231199E+001 1.393722E+005 1.394145E+005
2.750000E-002 6.342054E+001 2.347708E+005 2.348342E+005 2.317949E+001 9.791622E+004 9.793940E+004 4.024106E+001 1.368546E+005 1.368948E+005
3.000000E-002 6.342684E+001 2.301300E+005 2.301934E+005 2.488762E+001 9.598069E+004 9.600557E+004 3.853922E+001 1.341493E+005 1.341879E+005
3.250000E-002 6.399297E+001 2.251897E+005 2.252537E+005 2.687350E+001 9.392021E+004 9.394708E+004 3.711948E+001 1.312695E+005 1.313066E+005
3.500000E-002 6.507103E+001 2.199731E+005 2.200382E+005 2.915124E+001 9.174455E+004 9.177370E+004 3.591978E+001 1.282286E+005 1.282645E+005
...... ...... ...... ...... ...... ...... ...... ...... ...... ......
- The 1st line indicates the broadening parameter corresponding to this output file, in this case 0.06 eV.
- The 2nd and 3rd lines indicate the carrier density of electrons and holes respectively. The first number is the total number of free electrons/holes in the simulation cell, and the second number is the corresponding electron/hole density per cm\(^3\).
- The 4th and 5th lines indicate the output quantities and their unit, which will be introduced in details below.
Below the first 5 headlines, there are 10 Columns of data.
- The 1st column is the frequency \(\omega\). Here the first frequency ( \(\omega = 0\) ) is not printed, to avoid a division by zero error as can be seen from the Kubo-Greenwood formula in the Theory section.
Columns 2-10 are different components of the conductivity in units (\(\Omega \cdot {\rm cm}\))\(^{-1}\). In the Kubo-Greenwood formula, the summation is over all k-points \(\bf k\) and all energy eigenstate \(m\) and \(n\), i.e.
Therefore, one can separate different components from this summation and obtain:
- The 2nd column: inter-band ( \(m \neq n\) ) full conductivity (both electron + hole contributions), i.e. \(\sum_{{\bf k},(m\neq n)n}\)
- The 3rd column is intra-band ( \(m = n\) ) full conductivity (both electron + hole contributions), i.e. \(\sum_{{\bf k},(m = n)n}\)
- The 4th column is inter+intra-band full conductivity, i.e. \(\sum_{{\bf k}mn}\)
- The 5-7th columns are the corresponding inter, intra and inter+intra-band hole (only sum over the valence bands (VB) ) conductivity \(\sigma_h\) , i.e. \(\sum_{{\bf k}mn \in {\rm VB}}\).
- The 8-10th columns are the corresponding inter, intra and inter+intra-band electron (only sum over the conduction bands (CB) ) conductivity \(\sigma_e\) , i.e. \(\sum_{{\bf k}mn \in {\rm CB}}\).
The mobility files have exactly the same structure as the conductivity output files, but with corresponding mobilities in units \({ {\rm cm}^2}/{\rm Vs}\).
Warning
Currently only the inter-band contributions are thoroughly tested and validated, so please only use conductivity/mobility data in the 2,5 and 8th column in production calculations !
In this example, we are interested in the full conductivity of Al, thus we will plot the second column of the above file KG_full_L11_out_Gaussian_0.0600_a_a.out.
Visualize the result
The optical conductivity \(\sigma(\omega)\) can be plotted using any plotting software. Below, we provide an example using python with the matplotlib and numpy libraries.
Example python visualization code for \(6\times6\times6\) k-grid
import numpy as np
import matplotlib.pyplot as plt
filename = f"KG_full_L11_out_Gaussian_0.0600_a_a.out"
cond_al = np.loadtxt(filename)
# convert conductivity from 1/(Omega*cm) to 1/(Omega*m)
unit_trans = 100
fig,ax = plt.subplots()
ax.plot(cond_al[5:,0],cond_al[5:,1]*unit_trans/1e6
,color='tab:red'
, label=r"k-grid $6^3$ full SCF")
plt.xlim(0,5)
plt.ylim(0,7)
plt.legend()
plt.xlabel(r'$\hbar\omega$ (eV)',fontsize=12)
plt.ylabel(r'Conductivity $(10^{-6} \ {\Omega \cdot m})^{-1}$',fontsize=12)
plt.savefig("Conductivity_k666.pdf")
If we copy this to a file plot_conductivity_k666.py and execute with
Conductivity_k666.pdf which looks like:
Webinar Users
If you are performing this tutorial in the webinar, in order to visualise Conductivity_k666.pdf you first need to copy the file to your local machine from the cluster.
This can be done by running the following command on your local machine:
In comparison, the conductivity spectra of FCC Al reported in literature is:
In general, the shape of the spectrum fits quite well; slight differences in peak heights and positions can originate from the usage of different basis sets, xc functional, and broadening parameters.
2. 18x18x18 k-grid
Fourier interpolation
The KG calculation above was performed on a sparse k-grid using \(6 \times 6 \times 6\) points. For converged calculations, the Brillouin zone integration needs much denser k-grids to converge. Due to the locality of the real-space matrix elements, the evaluation of the KG formula can be performed on a dense k-grid via Fourier interpolation, even if the SCF calculation is performed on a relatively sparse k-grid. For instance, this technique is also used for band-structure and density-of-states calculations.
To showcase this, we now perform a KG calculation on a dense \(18 \times 18 \times 18\) k-grid via Fourier interpolation.
Let us create a directory k181818_fourier_interp and move to it:
geometry.in and control.in files from the previous calculation
and modify the dos_kgrid_factors tag in control.in from
to
Here, the SCF calculation still uses a sparse \(6 \times 6 \times 6\) k-grid, but the tag dos_kgrid_factors 3 3 3 requests an interpolated k-grid that is three times denser in each direction (\(3 \times 6 = 18\)) for the evaluation of the KG formula. We can then submit and run this calculation as we did in the previous part with e.g.
SCF k-grid
Using Fourier interpolation saves both memory and computational time compared to a naive k_grid 18 18 18 setting, which we will compare the result to in the following. Let us first create a directory k181818_scf and go to it:
k_grid tag in control.in from
to
and submit and run the calculation as previous with e.g.
Visualize the result
Once finished the output files are in the same structure as described previously. The reliablity of the Fourier interpolation can be validated by comparing the results with
and full SCF We can move to the directory above thek181818_fourier_interp and k181818_scf folders, and use the following python script plot_conductivity_k181818.py to visualise the conductivity spectra of the two approaches, like we did in the \(6\times6\times6\) k-grid case.
Example python visualization code for \(18\times18\times18\) k-grid
import numpy as np
import matplotlib.pyplot as plt
filename2 = f"k181818_fourier_interp/KG_full_L11_out_Gaussian_0.0600_a_a.out"
filename3 = f"k181818_scf/KG_full_L11_out_Gaussian_0.0600_a_a.out"
cond_al2 = np.loadtxt(filename2)
cond_al3 = np.loadtxt(filename3)
# convert conductivity from 1/(Omega*cm) to 1/(Omega*m)
unit_trans = 100
fig,ax = plt.subplots()
ax.plot(cond_al3[5::10,0],cond_al3[5::10,1]*unit_trans/1e6
,'o',color='blue', markerfacecolor='white'
, label=r"k-grid $18^3$ full SCF")
ax.plot(cond_al2[5:,0],cond_al2[5:,1]*unit_trans/1e6
,color='red'
, label=r"k-grid $18^3$ via Fourier interpolation")
plt.xlim(0,5)
plt.ylim(0,2.2)
plt.legend()
plt.xlabel(r'$\hbar\omega$ (eV)',fontsize=12)
plt.ylabel(r'Conductivity $(10^{-6} \ {\Omega \cdot m})^{-1}$',fontsize=12)
plt.savefig("Conductivity_k181818.pdf")
After executing
you should have a fileConductivity_k181818.pdf which looks like the following
Webinar users will once more have to copy this file over to their local machine before visualising, as decribed above.
The figure demonstrates the reliability of the Fourier interpolation. In comparison, the literature result is:
In general, the shape of the spectrum fits quite well; slight differences in peak heights and positions can originate from the usage of different basis sets, xc functional, and broadening parameters.
Warning
In production calculations, the SCF k-grid must be converged before performing a Fourier interpolation. The converged SCF k-grid is system-dependent and needs to be tested. For instance, a setting of k_grid 1 1 1 would be too sparse to produce a correct Fourier interpolation.
An important rule of thumb is: for larger cells, a sparser SCF k_grid is sufficient.
In production calculations, the species defaults should also be converged, for example computing with intermediate or tight defaults.
The light basis used here is purely for illustrative purposes, to enable calculations to be done in a relatively quick timeframe.
Solutions
You find all the solution to all the above exercises by clicking on the button below.
Appendix: Parallelization Scheme
In a KG calculation, the total number of k-points calculated on each node nkp_per_node is:
On each node, the k-points are calculated in loops. The number of k-points calculated each time is the total number of SCF k-grid (nkp_scf = 6^3) or the number of tasks on each node (n_tasks) when n_tasks < nkp_scf.
Thus, use more HPC nodes can speed up the KG calcualtion.
- For example, use 3 HPC nodes and 36 tasks on each node, the KG calculation with Fourier interpolation will finish in \(6^3\times 3^3 / 3 /36 = 54\) loops. Now if one increase the number of HPC nodes to 6, the same KG calculation will finsih in \(27\) loops.
Warning
This is different from the parallelization scheme in the density of state (DOS) calculation with Fourier interpolation, where FHI-aims only calculate nkp_scf k-points each time, no matter how many nodes one use. This is to ensure that one will not encounter any memory problem when doing Fourier interpolation, since we know that SCF can be done successfully with nkp_scf k-points.
But in KG implementation, we change to the above mentioned parallelization scheme to make our calculation more efficient. Aware that, with this scheme although the efficiency is increased one may encounter memory problem for large system.