RT-TDDFT
We refer the interested reader to the following paper detailing the Ehrenfest dynamics with real-time TDDFT implementation within FHI-aims:
In Ehrenfest dynamics, non-adiabatically coupled motion between the electronic and ionic subsystems is simulated. Mobile atoms and thus mobile basis functions lead to a modified differential equation for the electrons. A short summary, including the relevant equations, is provided within the manual.
In this tutorial we cover the basics of running RT-TDDFT with Ehrenfest dynamics in FHI-aims through a simple example: bond dissociation of \(H_2\).
Workflow
We are going to carry out:
- Prepare geometry.in file: Set up \(H_2\) structure.
- Ground-State DFT Calculation: Compute a converged ground state for \(H_2\) molecule.
- RT-TDDFT Calculation with Ehrenfest Dynamics Perform RT-TDDFT calculation, using the length gauge to induce dissociation of the molecule.
- Visualize Results Using the outputted trajectory file, visualize the bond breaking.
Materials
Related starting files
| file or dir name | description |
|---|---|
control.in |
FHI-aims control file, used as a template for all calculations |
geometry.in |
FHI-aims geometry file of \(H_2\). |
submit.sh |
Example script for submission of calculations to HPC (only an example, all calculations can be run locally) |
Prepare geometry.in file
It's important to start a calculation from a reasonable structure. The H_2 molecule is simple enough to start from an educated guess, knowing the bond length should be 0.74 Angstrom. If we define the bond in the z direction our file would be:
atom 0.00 0.00 0.00 H
atom 0.00 0.00 0.74 H
Prepare control.in file
The control.in for running with Ehrenfest dynamics will look very similar to the other RT-TDDFT examples provided, with a few additional runtime commands. Here we cover the basic choices and a few of the output options. Additional details are provided in the manual.
# 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 and Ehrenfest
RT_TDDFT_input_units atomic
RT_TDDFT_propagation 300 0.2 0.4
RT_TDDFT_propagator exponential_midpoint
RT_TDDFT_td_field 0 100 3 0.1 0 50 10 0 0 0.2
RT_TDDFT_ehrenfest default 0.2 0.4
RT_TDDFT_td_field_gauge length
RT_TDDFT_output_energies T T
RT_TDDFT_ehrenfest_output_trajectory T T
RT_TDDFT_write_file_prefix H2
RT_TDDFT_output_level 2
#intermediate species defaults
################################################################################
#
# FHI-aims code project
# Volker Blum, 2017
#
# Suggested "intermediate" defaults for H atom (to be pasted into control.in file)
#
################################################################################
species H
# global species definitions
nucleus 1
mass 1.00794
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 1e-4
#
radial_base 24 7.0
radial_multiplier 2
angular_grids specified
division 0.1930 50
division 0.3175 110
division 0.4293 194
division 0.5066 302
division 0.5626 434
# division 0.5922 590
# division 0.6227 974
# division 0.6868 1202
# outer_grid 770
outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 1 s 1.
# ion occupancy
ion_occ 1 s 0.5
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Basis constructed for dimers: 0.5 A, 0.7 A, 1.0 A, 1.5 A, 2.5 A
#
################################################################################
# "First tier" - improvements: -1014.90 meV to -62.69 meV
hydro 2 s 2.1
hydro 2 p 3.5
# "Second tier" - improvements: -12.89 meV to -1.83 meV
hydro 1 s 0.85
# hydro 2 p 3.7
# hydro 2 s 1.2
for_aux hydro 3 d 7
# "Third tier" - improvements: -0.25 meV to -0.12 meV
# hydro 4 f 11.2
# hydro 3 p 4.8
# hydro 4 d 9
# hydro 3 s 3.2
The first two sections of the file are identical to the isolated water case in tutorial-2. Only a few commands are necessary to run Ehrenfest dynamics with RT-TDDFT.
RT_TDDFT_ehrenfest default tells FHI-aims to run the RT-TDDFT simulation with Ehrenfest dynamics, and at the moment default must be chossen to specify the scheme. The final two inputs (0.2 and 0.4) specify the time step for geometry update and force computation and the printing of the trajectory. Currently, running with identical time steps for Ehrenfest and RT-TDDFT is a good idea to prevent and errors.
The only other Ehrenfest input we are including here is RT_TDDFT_ehrenfest_output_trajectory T T. This tells FHI-aims to print the time-dependent coordinates, velocities and forces to a file.
Additionally, we have made some changes within the RT-TDDFT section to simulate the bond-breaking.
We are now using the exponential_midpoint propagator to propagate the wavefunctions, simply to highlight there are multiple propagator options availiable within FHI-aims. Additional options are avaliable in the manual.
For RT_TDDFT_propagation we have changed the total simulation time to 300 a.u. the time step for RT-TDDFT to 0.2 a.u. and the timestep for printing the energy to 0.4 a.u. i.e. every second step.
Here we have different time steps for progagation and printing to highlight that data need not be printed for every step. Depending on the application, printing less data may be the optimal route.
The most complicated difference here is the change to the external potenital. We are now using RT_TDDFT_td_field 0 100 3 0.1 0 50 10 0 0 0.2 to apply the field. The first two inputs specify starting application of the electric field at t=0 a.u. and ending it at t=100 a.u.. The third input the 3 tells aims to use localized field oscillation via sine-modulated Gaussian. The 0.1 0 50 10 tell aims to apply the field with frequency of 0.1 a.u. centered at 50 a.u. time with a width of 10 a.u.. The final three inputs 0 0 0.2 tell aims to apply the pulse with strength 0.2 a.u. in the z cartesian direction. More details on RT_TDDFT_td_field are included in the manual.
Finally, we are using the intermediate specie default for hydrogen to create the basis.
Run the simulation
Now we are ready to run the simulation! To run on an HPC, change the path of the aims.x executable along with other system specific parameters and run,
sbatch submit.sh
N processes by typing
srun -n N aims.x > H2.out 2>&1
Once the calculation is finished, the directory will contain an H2.out file with Have a nice day printed near the end upon successful completion.
Check and visualize the result
You should have the following files:
Output files
| file or dir name | description |
|---|---|
H2.rt-tddft.ext-field.dat |
External field in time |
H2.rt-tddft.energy.dat |
Change in energy throughout RT-TDDFT simulation |
H2.rt-tddft-ehrenfest.trajectory.xyz |
Time-dependent coordinates, velocities and forces |
H2.out |
FHI-aims output stream |
We can use the printed H2.rt-tddft-ehrenfest.trajectory.xyz file to visualize the bond breaking. View the first 4 lines and the last for lines of the file. You'll see something like this:
# first four lines
[chrsshpr@dogwood-login2 Tutorial-4]$ ls
2
# Time (fs) = 0.00000000 | Format: x y z vx vy vz Fx Fy Fz | Units: [R] = A, [V] = A/ps, [F] = eV/A
H 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 -0.8162441356E-33 -0.3516573040E-30 -0.4091614821E+00
H 0.0000000000E+00 0.0000000000E+00 0.7400000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.1576402252E-30 -0.3516573040E-30 0.4091614821E+00
...
# last four lines
# Time (fs) = 7.26632830 | Format: x y z vx vy vz Fx Fy Fz | Units: [R] = A, [V] = A/ps, [F] = eV/A
H 0.2280863015E-15 0.7781705430E-15 -0.8017080191E+00 0.1557963492E-12 0.1874865692E-12 -0.1662735068E+03 -0.3948680256E-12 0.1262998012E-11 -0.5554473806E+01
H 0.2602893088E-15 0.1189165653E-14 0.1542625810E+01 0.1009047665E-13 0.6486078452E-12 0.1613178071E+03 -0.2739859416E-12 -0.1449727193E-12 -0.4832057210E+01
If one has access to a visualization program, such as VMD, we can make a movie using the trajectory file to see the dissociation. An example is shown below:
We are finished the RT-TDDFT tutorial!