Absorption spectra of a water molecule using RT-TDDFT

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:

1. 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.
2. Ground-State DFT Calculation: Compute a converged ground state for the isolated water molecule
3. 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.
4. 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, with a magntiude small enough to ensure only linear response of the system. This electric field, which is suddenly switched on and then off, can be Fourier transformed into the frequency domain and shown to contain all frequencies with equal weight. Therefore, it should excite all \(z\)-polarized electronic eigenmodes (normal modes of vibration with corresponding eigenvalues or energies) of the molecule. As discussed later, additional simulations with delta kicks along the \(x\) and \(y\) directions can be performed to excite all eigenmodes of the molecule. 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.

Following the kick, we allow the system to evolve freely in time, and output the energy and dipole into files with the prefix z to denote the response to the delta kick in the z-direction. The dipole (as a function of time) is the key observable here as it contains the linear response of the function to the electric field, detailed in the processing results section.

Finally the basis set, while not shown here, is appended to the end of the files. We are using the Tier2_aug2 specifies defaults. 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 simulations with a delta kick. However, to obtain a complete absorption spectrum , the system must be perturbed along all three Cartesian directions, not just z, as done in the example above. The full absorption spectrum requires the response along all Cartesian directions, since experimental absorption measurements average over light polarized in all directions. For non-symmetric systems, different responses may couple to x-, y-, or z-polarized fields, so including all components ensures a direct comparison with experiment.

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. In the case of the latter one can change the path the aims.x executable to your compiled version, update other system specific parameters, and in the terminal 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 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 obtain the time-dependent dipole response in all three Cartesian directions. This response is directly analogous to the fluctuation–dissipation theorem: in both cases, the system’s natural fluctuations and its response to a small external perturbation encode the same information about its intrinsic dynamics and excitation spectrum. For an everyday analogy, consider striking a bell with a single hammer hit. The resulting combination of harmonics is analogous to the dipole response of the system: each natural frequency of the bell is excited, and the bell “rings” as a combination of these modes. By Fourier transforming this ringing, we can determine the contribution of each frequency to the overall sound.

In the molecular case, the time-dependent dipole response plays the same role. Fourier transforming this response yields a frequency-dependent spectrum, and the imaginary part of this spectrum corresponds to the absorption spectrum. This reveals which electronic or vibrational eigenmodes of the molecule are most strongly excited by the perturbation.

Here, we proceed as described above by Fourier transforming the time-dependent dipole moments to obtain the complete absorption spectrum, 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.

We have recorded the dipole response as a function of time, i.e. harmonics in terms of the analogy above. Our goal is now to determine the frequencies at which these oscillations occur. However, the simulation was stopped at an arbitrary point in time (here, after 250 a.u.), which can introduce artifacts in the Fourier transform. These artifacts appear as spurious oscillations or “wiggles.” To mitigate this, we apply a damping function to the time-domain signal (dipole) prior to performing the Fourier transform, resulting in a cleaner and more reliable spectrum.

This produces the following files in the x_kick directory. The Python script writes the output files to the same directory as the first provided dipole file:

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 6 -e 0 40

The command runs the script with the absorption stength data as input. The t tells the script it was a 6fs simulation and -e 0 40 tells the script to plot energies in the range 0 to 40 eV.

If gnuplot or an X11 viewer are unavaliable on a cluster, one can use a juypter notebook and copy the following script:

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt("x_kick/x.rt-tddft.dipole.abs_strength.dat")
x = data[:,0] 
y = data[:,1]
plt.plot(x, y)
plt.xlim(0, 40)
plt.ylim(0, 0.1)
plt.xlabel("Energy (eV)")
plt.ylabel("S(w) (1/eV)")
plt.show()

We can now see the absorption spectrum for an energy range of 0 to 40 eV! The spectra should look like the following:

The spectra will look somewhat noisey and uneven and that's from only running the simulation for 6fs. For a smoother, less noiesy spectra, a longer simulation time is needed. The following is a spectra from a 24fs (1000 a.u.) simulation:

To run the longer time simulation, simply change the RT_TDDFT_propagation parameter in all the control.in files and rerun the simulations!

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.