Absorption spectra of a water molecule using RT-TDDFT

drawing

In this tutorial we are going to simulate the absorption spectra of an isolated water molecule using the length gauge to introduce an electric field.

Workflow

We are going to carry out:

Geometry Optimization: Perform a geometry optimization to obtain the equilibrium structure of the water. This is covered in the Basics tutorial and for simplicity we start from the pbe optimized structure.
Ground-State DFT Calculation: Compute a converged ground state for the isolated water molecule
RT-TDDFT Calculation with Length Gauge Perform 3 RT-TDDFT calculations, using the length gauge to introduce a delta kick in the x, y and z direction seperatly.
Calculate Dipole Strength Equation: Using the outputed dipole moments from the three simulations, calculate the absorption strength function for water.

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 water.
`submit.sh`	Example script for submission of calculations to HPC (only an example, all calculations can be run locally)

Geometry Optimization

As with all calculations it's important to start from a reasonable structure, and isolated water is no different. Rather than perform a geometry optimization here, instead we use the pbe optimized structure provided in the Basics of Running AIMS tutorial. Please refer to this tutorial for details on the geometry optimization. Using said structure we can create the geometry.in file:

geometry.in

atom      -0.00000000     -0.07320103      0.00000000 O
atom       0.76731112     -0.67039949      0.00000000 H
atom      -0.76731112     -0.67039949      0.00000000 H

Using a text editor such as vim one can view the contents of this file. More details on the syntax of the geometry.in can be found in the manual.

Setting up the simulation

First, using the text editor, open the control.in file to see what parameters are used for both the DFT and RT-TDDFT calculations.

control.in

#  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-4
    sc_accuracy_eev    1E-2
    sc_accuracy_etot   1E-6
    sc_iter_limit      100
# RT-TDDFT
    RT_TDDFT_input_units atomic
    RT_TDDFT_propagation 250 0.1 0.1
    RT_TDDFT_propagator crank_nicolson
    RT_TDDFT_td_field 0 0.1 1 0 0 0 0 0 0 0.001
    RT_TDDFT_td_field_gauge length
    RT_TDDFT_output_energies T T
    RT_TDDFT_output_dipole T T
    RT_TDDFT_write_file_prefix z
    RT_TDDFT_output_level 2

The first section of the control file (denoted by physical model) is used to describe the physical system we are looking at. We are using the PBE XC functional, no spin, a non-relatistic calculation and an uncharged water molecule.

The second section of the control file contains the parameters for convergence of the ground state. The occupation_type gaussian 0.01 specifies a gaussian broadening scheme to find the Fermi level and the KS state occupations using the paramters recommened in the manual to prevent any fractional occupations.

mixer pulay specifies the default pulay mixing algorithm for the electron density to achieve SCF convergence. The final four lines of the SCF section, all starting with sc_ denote convergence parameters for the SCF calculation; charge density, eigenvalue sum, total energy and number of iterations respectively. Additonal information can be found in the manual.

The third sections contains the RT-TDDFT parameters. Here we are using atomic units, and RT_TDDFT_propagation states we will propagate the system for 250 atomic units of time using a 0.1 a.u. time step for both the simulation and the printing of observables. RT_TDDFT_propagator specifies the numerical integration scheme, here we are using Crank-Nicolson. The timestep and propagator should be chossen with care, and more information can be found in the manual.

RT_TDDFT_td_field RT_TDDFT_td_field_gauge specify the gauge and the external field applied to the system.

For RT_TDDFT_td_field the usage is as follows:

RT_TDDFT_td_field t_start t_end type freq cycle center width Ex Ey Ez

Full documentation on the input can be found in the manual. Here we are specifying a constant external field starting at \(t=0\) and ending at \(t=0.1\), applied in the z-direction at a stength of 0.001 a.u. This ammounts to a delta kick in the z-direction, exciting all frequencies. For the gauge we are using the length gauge, which is the default for isolated systems. Details on gauge selection can be found in the manual.

Additonally, we are outputting the energy and dipole into files with the prefix z to denote the response to the delta kick in the z-direction.

Finally, we are using the Tier2_aug2 specifies defaults, which are appeneded to the end of the file, for the basis set. One can also use atom-centered orbital (NAO) basis sets and numerical definitions provided in species_defaults.

Species Defaults: Choice of the numerical settings, including the basis set

FHI-aims uses numerically tabulated NAO basis sets to represent the quantum mechanical orbitals and wave functions. These basis functions are specific to each chemical element ("species") used in the simulation.

FHI-aims also employs a variety of other numerical parameters that are specific to each atom and/or chemical element: Angular-momentum expansion order of the electrostatic (Hartree) potential, atom-centered real-space integration grids, etc.

Predefined files that contain these numerical definitions for each chemical element are located in species_defaults directory of the FHI-aims code distribution. At least four different levels of precision are provided for each chemical element: light, intermediate, tight, and really_tight. A more accurate choice of defaults means a higher computational cost of the simulation.

Running the Simulation(s)

Now we have a control.in template for running RT-TDDFT simulation with a delta kick. However, for a complete absorption spectra we need to perturbe the system in all 3 Cartesian directions, not just z, as done in the control.in file above. To apply the delta kick in a different direction, we can adjust theRT_TDDFT_td_field command in the input. This is done in the three kick subdirectories, and each contains a delta kick in the cartessian direction listed in the directory name. Using a text editor one can see the only difference between the control.in files within each directory is the RT_TDDFT_td_field line, where the direction of the delta kick is changed. Prefixes for the outputted files are also changed for simplicity.

Now we are ready to run the simulations! Each job can be run locally or on a HPC. For the latter submit.sh provides a template for running all three jobs together on an HPC.

In the case of the latter one can change the path the aims.x executable along with other system specific parameters and run,

sbatch submit.sh

This will submit all simulations.

For the case of the former, one can run a parallel simulation with N processes by typing

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

This command must be excuted within each kick directory to run the delta kick in the 3 Cartesian directions. While this can be run locally, the calculation will be very slow. In this case one should switch to light species defaults to test the code rather then running with the much larger Tier2_aug2 basis.

Once the calculations are finished, each subdirectory should contain an aims.out file with Have a nice day printed at the end if the calculation finished successfully. Additionally each subdirectory will contain (with different prefixes depending on the directory):

Output files

file or dir name	description
`x.rt-tddft.dipole.dat`	Electronic dipole in time
`x.rt-tddft.ext-field.dat`	External field in time
`x.rt-tddft.energy.dat`	Change in energy throughout RT-TDDFT simulation

Note on running

The binary name aims.x should be replaced with whatever is the name of the FHI-aims binary file compiled by you (including the corresponding path, i.e., the location of the directory in which that file is located).

The mpirun command facilitates the parallel execution of the code and can have different names that depend on the particular computer system and/or queueing system used. The actual name and usage are usually documented by the computing center in question.

Processing the results

Once all simulations have finished, we can process the indiviudal dipoles to get a complete absorption spectra by running:

python eval_tddft.py -pol x_kick/x.rt-tddft.dipole.dat y_kick/y.rt-tddft.dipole.dat z_kick/z.rt-tddft.dipole.dat -field x_kick/x.rt-tddft.ext-field.dat -fa -d 0 poly -w

The -pol argument tells the script to calculate the absorption spectra from the three dipole files, and -field specifies the strength of the applied field for each delta kick. We only need to pass the code one ext-field file, as long as we applied a delta kick of identical magnitude in each direction (which we did). The -fa argument tells the script to take the FT of the field numerically, the -d 0 poly arugument specifies a polynomial damping function starting at \(t=0\) to reduce wiggles in the specta and -w tells the script to write results to file. The theory behind

This produces the following files in the x_kick directory:

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

Plotting the spectra

Now the absorption strength is written to, we can plot the spectrum. For a quick check one can view the spectra in gnuplot using the following command.

cd x_kick
gnuplot
plot 'x.rt-tddft.dipole.abs_strength.dat' u 1:2 w l; set xrange [0:40]

For atomic and molecular systems, determining the precise spectral location of electronic transitions as well as their oscillator strengths is a major concern. The script ifit abs.py offers graphical support for this task, as described above. One can run the script and view the same energy range using the following command:

python ifit_abs.py -f x_kick/x.rt-tddft.dipole.abs_strength.dat -t 24 -e 0 40

The command runs the script with the absorption stength data as input. The t tells the script it was a 24fs simulation and -e 0 40 tells the script to plot energies in the range 0 to 40 eV. We can now see the absorption spectrum for an energy range of 0 to 40 eV! The spectra should look like the following:

drawing

To determine the exact energy of individual peaks, click on one of the highest peaks within the spectra. For comparison, the osscilator strengths from a LR-TDDFT simulation of water are provided within the solutions folder.