Skip to content

Absorption Spectra of Crystalline Silicon

drawing

This tutorial shares many similarities with the water molecule tutorial, albeit we are now investigating crysalline silicon. Our calculations thus have to be adjusted to deal with a system in periodic boundary conditions.

Workflow

We are going to carry out:

  1. Geometry Optimization: Perform a geometry optimization to obtain the equilibrium structure. How this is done is covered in the Basics tutorial and we start from the PBE optimized structure.
  2. Ground-State DFT Calculation: Compute a converged ground state for the silicon.
  3. RT-TDDFT Calculation with Velocity Gauge Perform 3 RT-TDDFT calculations, using the velocity gauge to introduce a delta kick in the x, y and z direction seperatly.
  4. Calculate Absorption Spectra: Using the outputed current densities from the three simulations, calculate the imaginary dielectric for silicon.

Materials

Input files can be found here

Related starting files

file or dir name description
eval_tddft.py Script recieves RT-TDDFT data from FHI-aims, calculates values for an explicitly given electrical field, and finally computes the time-dependent polarisability and the absorption strength function. This can be found in the FHIaims directory (in the FHIaims/utitilies/rt-tddft subdirectory). Alternatively, a copy is provided with the tutorial.
ifit_abs.py Script offers graphical support for this task: by clicking on a peak in the spectrum, an algorithm based on Pad´e approximants bruner-2016 outputs the precise peak position. This can be found in the FHIaims directory (in the FHIaims/utitilies/rt-tddft subdirectory). Alternatively, a copy is provided with the tutorial.
control.in FHI-aims control file, used as a template for all calculations
geometry.in FHI-aims geometry file of silicon.
submit.sh Example script for submission of calculations to HPC (only an example, all calculations can be run locally)

Key differences when considering periodic systems

  1. The input, RT_TDDFT_td_field_gauge length, does not work for periodic systems, thus we must move to a the velocity gauge. We can do this by employing the keyword RT_TDDFT_td_field_gauge velocity.
  2. We now must include keywords in the control.in file and geometry.in file specifying the k-grid and periodicity of the system. This is covered within the Basics tutorial.
  3. Instead of computing the dipole strength function, we are now computing the imaginary dielectric function to determine the absorption spectra.

Prepare the geometry.in file

The syntax for the geometry.in file remains the same as the previous section, with only a few additions:

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

The first three lines starting with lattice_vector define the three unit cell vectors \(a_i\) (\(i\)=1,2,3). The following lines define the position and the species of the atoms for this system. Again we are using the PBE optimized geometry and lattice from the Basics tutorial.

Fractional coordinates

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.

The position of atoms can also be defined in so-called fractional coordinates (atom_frac) \(f_1\), \(f_2\), \(f_3\). Fractional coordinates are especially convenient when working with high-symmetry crystals.

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

Practical Recommendation

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.

Setting up and Running the Simulation(s)

With the optimized geometry, we can immediately set up our DFT/RT-TDDFT simulations(s). The control.in file is provided below (omitting light species defaults used from the basis set): 

#  Physical model
  xc                 pbe
  spin               none
  relativistic       none
  charge             0.
#  SCF convergence
  occupation_type    gaussian 0.01
  mixer              pulay
    n_max_pulay        10
    charge_mix_param   0.5
  sc_accuracy_rho    1E-5
  sc_accuracy_eev    1E-3
  sc_accuracy_etot   1E-6
  sc_iter_limit      100
#  For periodic boundary conditions
  symmetry_reduced_k_grid .false.
  k_grid 4 4 4
# RT-TDDFT
        RT_TDDFT_input_units atomic
        RT_TDDFT_propagation 500 0.1 0.1
        RT_TDDFT_propagator crank_nicolson
        RT_TDDFT_td_field 0 0 1 0 0 0 0 0 0 0.001
        RT_TDDFT_td_field_gauge velocity
        RT_TDDFT_output_energies T T
        RT_TDDFT_output_current T T        
        RT_TDDFT_write_file_prefix z
        RT_TDDFT_output_level 2

While this file looks very similar to the molecular cases, there are some key differences. In the case of periodic systems, we need to specify the grid of crystal momentum (\(k\)) vectors used to sample the reciprocal lattice of the system. In FHI-aims, this \(k\)-space integration grid is specified by the keyword 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_density is 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\). Using k_grid_density prevents 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_grid keyword, 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.

For a fully converged spectra we need a very dense k-grid of 16x16x16 with and offset from \(\Gamma\)-centering. However, in this tutorial use a small 4x4x4 k-grid without symmetry reduction. Even with a 4x4x4 \(\Gamma\)-centered k-grid and light species defaults for the basis, a single simulation (24 fs) has a very high computational cost, and larger k-grid simulations aren't possible without access to an HPC. Thus although the k-grid is to small for an accurate result, we perform the example with it to demonsrate the steps and provide the result when using a 16x16x16 k-grid to show how important using a large enough k-grid is. A longer simulation is also needed for a smooth and converged spectra. Additionally, we will only perform the z-directional kick, and the inputs and outputs for application of the external field in the x and y directions is provided to save time.

The other key change to our input is the usage of the velocity gauge rather than the length gauge. Both gauges are physically equivalent but have different technical implications. In this implementation, only the velocity gauge can be applied in case of periodic systems. Finally, since we are using the velocity gauge, our key observable is the time dependent current and thus we add T_TDDFT_output_current T T to output the current. Now we can run the simmulation. The job can be run locally or on a HPC. For the latter submit.sh provides a template for the simulation on an HPC

Identical to the molecular case, change the path of the aims.x executable along with other system specific parameters and run, 

sbatch submit.sh
To run locally, one can run a parallel simulation with N processes by typing

cd z_kick
srun -n N aims.x > aims.out 2>&1

This command must be excuted within the z-kick subdirectory. Be aware this will take updwards of 15 minutes. It is highly recommended to run this on an HPC. The simulation will take a few minutes. The x-kick and y-kick directories contain compeleted simulation results.

Once the calculation is finished, the directory will contain an aims.out file with Have a nice day printed near the end upon successful completion. Additionally the z_kick subdirectory will contain:

Output files

file or dir name description
z.rt-tddft.ext-field.dat External field in time
z.rt-tddft.energy.dat Change in energy throughout RT-TDDFT simulation
z.rt-tddft.current.dat Change in current throughout RT-TDDFT simulation

Process and plot the results

Again we will use the eval_tddft.py file to process the data. Since our system is periodic and we used the velocity gauge, we now will input our time dependent currents: 

cd ../
python eval_tddft.py -pol x_kick/x.rt-tddft.current.dat y_kick/y.rt-tddft.current.dat z_kick/z.rt-tddft.current.dat -field z_kick/z.rt-tddft.ext-field.dat -fa -d 0 poly -w
Beside using the current the input is identical to the molecular case. The same 4 files are written:

Output files

file or dir name description
x.rt-tddft.dipole.abs_strength.dat absorption strength
x.rt-tddft.dipole.dipole_ft.dat dipole FT
x.rt-tddft.dipole.polarisability.dat x,y and z polarisability
x.rt-tddft.dipole.pow_spec.dat power spectrum

Now we can view our spectra! Since we have a periodic system, we want to use the polarisability file. 

python ifit_abs.py -f z_kick/z.rt-tddft.current.polarisability.dat -t 12 -e 0 20
The command runs the script with the polarisability as input. The -t tells the script it was a 12fs simulation and -e 0 20 tells the script to plot energies in the range 0 to 20 eV. Our spectra looks like the following:

drawing

The spectra looks quite different than those shown here. This is due to our usage of a 4x4x4 k-grid. If we were to use a larger k-grid, such as 16x16x16, our result would look much better:

drawing

However this is significantly more computationally expensive, and thus we use a small grid in the tutorial.