Real-Time Time-Dependent DFT with NEO
In this tutorial, we simulate an excited-state proton transfer in an o-hydroxybenzaldehyde (oHBA) molecule, demonstrating how to run real-time time-dependent DFT with NEO (RT-NEO-TDDFT).
More details about RT-NEO-TDDFT calculations can be found in the following papers:
J. Xu, et al., "First-Principles Approach for Coupled Quantum Dynamics of Electrons and Protons in Heterogeneous Systems", Phys. Rev. Lett. 131, 238002 (2023) https://doi.org/10.1103/PhysRevLett.131.238002
J. Xu, et al., "Lagrangian Formulation of Nuclear-Electronic Orbital Ehrenfest Dynamics with Real-time TDDFT for Extended Periodic Systems". https://doi.org/10.48550/arXiv.2407.18842
Workflow
We are going to cover:
- Prepare the
geometry.infile: Set up the structural file, including ghost basis centers. - Prepare the
Control.infile: Configure parameters for the RT-NEO-TDDFT calculation with HOMO to LUMO excitation. A validcontrol.inis provided at the end. - Extra information: Additional optional parameters and instructions for Ehrefest dynamics.
- Check the results: go throught the output files for energies, positions, densities.
Prepare the geometry.in file
In this tutorial, we used restricted optimized structure of oHBA in its excited state.
To model the proton transfer process with RT-NEO-TDDFT and ensure sufficient basis functions,
we included three additional ghost proton basis function centers (indicated by empty), along with the proton basis function center located at the optimized proton position in the excited state.
An example of geometry.in file is provided below:
atom -1.310008 1.258755 0.000000 C
atom 0.019289 0.780580 0.000000 C
atom 0.322586 -0.621636 0.000000 C
atom -0.761283 -1.465586 0.000000 C
atom -2.125342 -0.985468 0.000000 C
atom -2.398971 0.362640 0.000000 C
atom 1.008708 1.658363 0.000000 O
atom 1.725826 -1.062678 0.000000 C
atom 2.682115 -0.219090 0.000000 O
atom -3.414260 0.733314 0.000000 H
atom -0.596742 -2.537818 0.000000 H
atom -2.926679 -1.713458 0.000000 H
atom -1.459388 2.331114 0.000000 H
atom 1.924136 -2.138192 0.000000 H
atom 1.844187 1.135044 0.000000 H
empty 1.981048 0.913868 0.000000 H
empty 2.117910 0.692694 0.000000 H
empty 2.366459 0.659685 0.000000 H
Empty Site in FHI-aims
In FHI-aims, a ghost atom or an empty site, where only the basis functions (but not the nucleus) of a given species are placed, is specified by the keywordempty.
An empty site is a pure basis function center with no mass, no charge, or other nuclear properties.
One should avoid using empty for physical reasons if structure relaxation is requested.
In RT-NEO-TDDFT, the empty option can be helpful for modeling proton transfer events when protons move far from their original locationm, and a single localized basis function center for each quantum proton is not sufficient to describe the dynamics.
Prepare the control.in file
Electron Input
Similar to running a NEO-DFT calculation with FHI-aims, the electronic part of an RT-NEO-TDDFT simulation is configured the same way as in a conventional real-time TDDFT (RT-TDDFT) calculation. Below is a brief overview of the electron input settings.
xc pbe
relativistic none
override_illconditioning .true.
RT_TDDFT_input_units atomic
RT_TDDFT_propagation 500 0.2 2
RT_TDDFT_propagator exponential_midpoint
RT_TDDFT_output_energies T T
RT_TDDFT_write_file_prefix oHBA
The DFT section describes the physical system, including exchange-correlation functional and relativistic effects.
We use the PBE exchange-correlation functional for this non-relativistic calculation.
The override_illconditioning flag is turned on to prevent numerical warnnings, as the empty sites (ghost atoms) are packed closely together.
The second section contains RT-TDDFT parameters.
The RT_TDDFT_input_units is a mandatory flag to avoid unit confusion, specifying atomic to ensure all inputs — such as time and external fields — are in atomic units.
The RT_TDDFT_propagation command initiates the RT-TDDFT after the DFT ground state is obtained.
It specifies the total simulation time for propagation, the time step for propagation, and the time step for output.
The general syntax is:
RT_TDDFT_propagation t_tot dt_prop dt_output
The RT_TDDFT_propagator specify the numerical integration scheme.
Here we use the expoential_midpoint options due to its simplicity.
By default, this is set to crank_nicolson.
Possible Integration Schemes - Electrons
- exponential_midpoint (EM): $$ U(t_k + \Delta t, t_k) = \exp\left( -i \Delta t \mathbf{S}^{-1} \mathbf{H}\left( t_k + \dfrac{\Delta t}{2} \right) \right) $$
- crank_nicolson (CN): $$ U(t_k + \Delta t, t_k) = \left[ 1 - i\frac{ \Delta t}{2} \mathbf{S}^{-1} \mathbf{H}\left( t_k + \dfrac{\Delta t}{2} \right) \right] \left[ 1 + i\frac{ \Delta t}{2} \mathbf{S}^{-1} \mathbf{H}\left( t_k + \dfrac{\Delta t}{2} \right) \right]^{-1} $$
- etrs (Enforced Time-Reverse symmetry): $$ U(t_k + \Delta t, t_k) = \exp\left( -i \frac{ \Delta t}{2} \mathbf{S}^{-1} \mathbf{H}\left( t_k + \Delta t \right) \right) \exp\left( -i \frac{ \Delta t}{2} \mathbf{S}^{-1} \mathbf{H}\left( t_k \right) \right) $$
- cfm4 (Commutator-Free Magnus Expansion 4)
- m4 (Magnus Expansion 4)
- runge_kutta_4 (RK4)
In all above equations:
- \(U(t_k + \Delta t, t_k)\) is the time propagation operator
- \(\Delta t\) is the time step for propagation,
- \(\mathbf{S}\) is the overlap matrix,
- \(\mathbf{H}\left( t \right)\) is the Hamiltonian matrix evaluated at time \(t\)
The RT_TDDFT_output_energies T T command ensures that energy information is written both to the FHI-aims output stream and to an additional file.
The RT_TDDFT_write_file_prefix parameter sets the prefix for all output files.
For more detailed instructions on configuring parameters for electrons in RT-TDDFT, please refer to the Real-Time Time-Dependent DFT Tutorial for FHI-aims.
NEO Input - Essential
To run RT-NEO-TDDFT, we can use the same settings as in the previous sections.
Most of the real-time propagation information (e.g., timestep, external field and simulation time) for the proton subsystem is read from the electron input in control.in.
NEO_type 2
NEO_basis_type primitive
NEO_epc_type 17
NEO_scf_proton_max_step 20
NEO_proton_hartree_rcut 20
# density fitting
NEO_exchange_type noRI
NEO_df_type 0
NEO_df_l_max 6
NEO_rt_update_occ_electron T
NEO_basis
H 15
s 4.0 1.000
p 4.0 1.000
H -16
s 4.0 1.000
p 4.0 1.000
H -17
s 4.0 1.000
p 4.0 1.000
H -18
s 4.0 1.000
p 4.0 1.000
NEO_basis_end
To initiate RT-NEO-TDDFT calculations, we set NEO_type to 2.
The NEO_epc_type is kept as epc17-2.
We use the same density fitting setting as defined in the previous section 3.
To reduce computational cost for this RT-NEO-TDDFT tutorial,
we use a minimal primitive Gaussian basis set, which contains only an s orbital and a p orbital with an exponent of 4.
In the NEO_basis section, this basis set is assigned to the transferring hydrogen atom between the two oxygen atoms and all three empty sites.
In practice, the basis set can be adjusted based on the desired accuracy.
NEO Basis for empty Sites
The syntax:
H [index] is used to assign basis functions for a quantum proton with index in geometry.in
In the case of empty sites in geometry.in, a negtive number of index is used for adding a ghost center without charge distributions.
To apply a initial exciation, we set the NEO_rt_update_occ_electron tag to .true..
This updates the electron occupation numbers, initiating the RT-NEO-TDDFT simulation with a specific electronic configuration.
The new occupations should be provided in a file called neo_occ_electron.in:
32 1
33 1
i_state_1 occ_1
i_state_2 occ_2
...
i_state_N occ_N
This let us initiate the RT-NEO-TDDFT by promoting an electron from state 32 to state 33, i.e., the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbitial (LUMO).
External Field for RT-NEO-TDDFT
One can also apply an external field to RT-NEO-TDDFT using the same RT_TDDFT_td_field keyword used for electrons.
The general syntax is:
RT_TDDFT_td_field t_start t_end type freq cycle center width Ex Ey Ez
Using the NEO_rt_ext_e_type flag to select:
- none: Ignore the exteral inputs for the proton subsystem
- scalar: Use scalar potential from RT_TDDFT_td_field for protons.
Before running the simulation, do not forget to attach the species defaults for O C and H.
In this tutorial, we will use a light defaults NAO basis to reduce computational cost.
A valid control.in for RT-NEO-TDDFT
xc pbe
override_illconditioning .true.
RT_TDDFT_input_units atomic
RT_TDDFT_propagation 500 0.2 2
RT_TDDFT_propagator exponential_midpoint
RT_TDDFT_output_energies T T
RT_TDDFT_write_file_prefix oHBA
NEO_type 2
NEO_basis_type primitive
NEO_epc_type 17
NEO_exchange_type noRI
# density fitting
NEO_df_type 0
NEO_df_l_max 6
# basis
# cut_off s
NEO_proton_hartree_rcut 20
# scf
NEO_scf_proton_max_step 20
NEO_rt_update_occ_electron T
NEO_basis
# pos in angstrom
H 15
s 4.0 1.000
p 4.0 1.000
H -16
s 4.0 1.000
p 4.0 1.000
H -17
s 4.0 1.000
p 4.0 1.000
H -18
s 4.0 1.000
p 4.0 1.000
NEO_basis_end
################################################################################
#
# FHI-aims code project
# VB, Fritz-Haber Institut, 2009
#
# Suggested "light" defaults for H 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.
#
################################################################################
species H
# global species definitions
nucleus 1
mass 1.00794
#
l_hartree 4
#
cut_pot 3.5 1.5 1.0
basis_dep_cutoff 1e-4
#
radial_base 24 5.0
radial_multiplier 1
angular_grids specified
division 0.2421 50
division 0.3822 110
division 0.4799 194
division 0.5341 302
# division 0.5626 434
# division 0.5922 590
# division 0.6542 770
# division 0.6868 1202
# outer_grid 770
outer_grid 302
################################################################################
#
# 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
# 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
################################################################################
#
# FHI-aims code project
# VB, Fritz-Haber Institut, 2009
#
# Suggested "light" defaults for C 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.
#
################################################################################
species C
# global species definitions
nucleus 6
mass 12.0107
#
l_hartree 4
#
cut_pot 3.5 1.5 1.0
basis_dep_cutoff 1e-4
#
radial_base 34 5.0
radial_multiplier 1
angular_grids specified
division 0.3326 50
division 0.5710 110
division 0.7727 194
division 0.8772 302
# division 0.9334 434
# division 0.9625 590
# division 0.9924 770
# division 1.0230 974
# division 1.4589 1202
# outer_grid 974
outer_grid 302
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 2.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 1.
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Constructed for dimers: 1.0 A, 1.25 A, 1.5 A, 2.0 A, 3.0 A
#
################################################################################
# "First tier" - improvements: -1214.57 meV to -155.61 meV
hydro 2 p 1.7
hydro 3 d 6
hydro 2 s 4.9
# "Second tier" - improvements: -67.75 meV to -5.23 meV
# hydro 4 f 9.8
# hydro 3 p 5.2
# hydro 3 s 4.3
# hydro 5 g 14.4
# hydro 3 d 6.2
# "Third tier" - improvements: -2.43 meV to -0.60 meV
# hydro 2 p 5.6
# hydro 2 s 1.4
# hydro 3 d 4.9
# hydro 4 f 11.2
# "Fourth tier" - improvements: -0.39 meV to -0.18 meV
# hydro 2 p 2.1
# hydro 5 g 16.4
# hydro 4 d 13.2
# hydro 3 s 13.6
# hydro 4 f 17.6
# Further basis functions - improvements: -0.08 meV and below
# hydro 3 s 2
# hydro 3 p 6
# hydro 4 d 20
################################################################################
#
# FHI-aims code project
# VB, Fritz-Haber Institut, 2009
#
# Suggested "light" defaults for O 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.
#
################################################################################
species O
# global species definitions
nucleus 8
mass 15.9994
#
l_hartree 4
#
cut_pot 3.5 1.5 1.0
basis_dep_cutoff 1e-4
#
radial_base 36 5.0
radial_multiplier 1
angular_grids specified
division 0.2659 50
division 0.4451 110
division 0.6052 194
division 0.7543 302
# division 0.8014 434
# division 0.8507 590
# division 0.8762 770
# division 0.9023 974
# division 1.2339 1202
# outer_grid 974
outer_grid 302
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 4.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 3.
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Constructed for dimers: 1.0 A, 1.208 A, 1.5 A, 2.0 A, 3.0 A
#
################################################################################
# "First tier" - improvements: -699.05 meV to -159.38 meV
hydro 2 p 1.8
hydro 3 d 7.6
hydro 3 s 6.4
# "Second tier" - improvements: -49.91 meV to -5.39 meV
# hydro 4 f 11.6
# hydro 3 p 6.2
# hydro 3 d 5.6
# hydro 5 g 17.6
# hydro 1 s 0.75
# "Third tier" - improvements: -2.83 meV to -0.50 meV
# ionic 2 p auto
# hydro 4 f 10.8
# hydro 4 d 4.7
# hydro 2 s 6.8
# "Fourth tier" - improvements: -0.40 meV to -0.12 meV
# hydro 3 p 5
# hydro 3 s 3.3
# hydro 5 g 15.6
# hydro 4 f 17.6
# hydro 4 d 14
# Further basis functions - -0.08 meV and below
# hydro 3 s 2.1
# hydro 4 d 11.6
# hydro 3 p 16
# hydro 2 s 17.2
Extra Information
This setion provides additional details that might be helpful for running RT-NEO-TDDFT simulations. You may skip this section if it is not directly applicable to your current task, as the following part is for informational purpose and is not used in this tutorial.
NEO Input - Extra
Below are some extra parameters that control the workflow of RT-NEO-TDDFT.
Proton Timestep
By default, the NEO subsystem uses the same time step as defined for the electron systems. However, you can control the proton timestep separately using the following command:
NEO_rt_propagate_n_step N
N electron steps.
Generally, N can be set greater than one due to the larger mass of protons as compared to electrons.
However, if N is increased, you should check energy conservation to ensure stability.
Proton Propagator
The integrator for the proton subsystem can be selected with:
NEO_rt_propagator_type type
- expoential_midpoint (EM): The exponential midpoint method. (Default)
- EM_iter : Uses an iteratively predict-and-correct method until the convergence criteria set by
NEO_rt_precor_tolare met - ruge_kutta_4 (RK4): The standard explicit Runge-Kutta 4 scheme. It works best with very small timesteps but is a good choice for accurate reference simulations due to its well-understood properties.
Proton Output
Output option includes:
- NEO energy terms (on by default).
- Dipole moment (on by default).
- Expectation value of the position operators for individual proton (Kohn-sham) states .
These can be controlled with the following tags:
NEO_rt_output_energyNEO_rt_output_dipoleNEO_rt_output_proton_center
To control the frequency of output, use:
NEO_rt_output_n_step N
N electron steps.
Note that the output frequency is based on electron steps, not proton steps.
NEO Ehrefest dynamics
NEO Ehrefest dynamcis
NEO Ehrefest is a recently implemented feature. Use at your own risk!
Ehrenfest dynamics is also avaiable in NEO. To initate Ehrenfest dynamics, use the general syntax:
RT_TDDFT_ehrenfest default dt_prop dt_output
Restarting a RT-NEO-TDDFT calculation
To restart an RT-NEO-TDDFT calculation, you can use the following steps:
-
Save necessary files: Ensure that you save all relevant files from the previous run, including electron and proton wavefunctions, coordinates and settings, etc. To write an restart file at a specific time
t_write_restart, use the general syntax incontrol.in:RT_TDDFT_restart_write t_write_restart -
Use restart tag: Same as how conventional RT-TDDFT restarts are managed, include the restart option in your
control.infile. The general syntax is:FHI-aims will run a NEO-SCF cycle, but the resulting eigenvectors and other necessary data will be overwritten from what is read from the restart files. Based on this, the RT-NEO-TDDFT simulation will be re-initialized. The restart file for NEO should be renamed asRT_TDDFT_restart_read restart_fileneo.rt-tddft.save.in -
Modify parameters: Adjust timesteps or any other parameters, if necessary, for restarting the time-dependent propagation.
Check the Results
After the FHI-aims calculation is complete, your folder should contain:
Output Files
| File name | Description |
|---|---|
oHBA.out |
FHI-aims output stream |
neo.rt-tddft.dipole.dat |
NEO output file for dipole moment |
neo.rt-tddft.energy.dat |
NEO output file for related energy terms |
oHBA.rt-tddft.energy.dat |
FHI-aims output file for energy |
The FHI-aims output stream contains running information about both the NEO-DFT SCF calculations and the real-time propagation. For example, a summary of the timing for RT-NEO-TDDFT can be found near the end of the file, as shown below:
Final timings (s): CPU | wall | [amount]
* Time propagation timings:
Matrix operations: 2.801E-01 2.750E-01 50020
-> Eigensolver calculations: 0.000E+00 0.000E+00 5002
Hamiltonian + density upd.: 3.090E+00 3.089E+00 2501
Total: 3.372E+00 3.368E+00
* Remaining timings:
Hamiltonian + density upd.: 3.095E+00 3.095E+00 5002
Output: 1.243E-02 6.000E-03
Total: 6.965E+00 6.960E+00
The oHBA.rt-tddft.energy.dat file contains all the energy terms during the real-time propagation.
The NEO.rt-tddft.energy.dat file contains the energy terms related to the proton subsystem.
The total energy (the last column of oHBA.rt-tddft.energy.dat) is a good parameter to monitor to ensure everything is proceeding as expected.
By plotting the total energy as a function of time, you can visualize how the energy fluctuate over the course of the simulation. Here is an example plot:
In this case, the solid blue line shows how energy drifted during the simulation, changing by 0.12 eV throughout. We can reduce this drift by using tighter numerical settings, such as a finer basis functions and smaller timesteps for both electron and protons.
Next, we can examine the trajectory of the quantum proton using the file neo.rt-tddft.dipole.dat.
Since we only have one quantum protons, the dipole moment also corresponds to the position operator of the proton.
Using the position of the two oxygen atoms from geometry.in, we can calcuate the distance between the quantum proton and the two oxygen atoms (acceptor and donor) and then plot them as a function of time.
Here is a PYTHON script (using matplotlib and numpy) to plot the distance:
import numpy as np
import matplotlib.pyplot as plt
O1 = np.asarray([1.008708, 1.658363, 0.000000])
O2 = np.asarray([2.682115,-0.219090, 0.000000])
bohr = 0.529177
dipole = np.loadtxt("[path_to_your_folder]/neo.rt-tddft.dipole.dat")
plt.xlabel("Time (fs)")
plt.ylabel("Distance (Å)")
plt.xlim([0,12])
plt.plot(dipole[:,0],np.sqrt(np.sum((dipole[:,1:4]*bohr-O1)**2,axis=1)),color="blue",linestyle="-",label="H...O2",lw=2)
plt.plot(dipole[:,0],np.sqrt(np.sum((dipole[:,1:4]*bohr-O2)**2,axis=1)),color="red",linestyle="-",label="H...O1",lw=2)
plt.legend(frameon=False)
It shows that we have successfully modeled the quantum proton moving from the donor to the acceptor under the HOMO to LUMO excitation using the RT-NEO-TDDFT simulation.
We are done with the RT-NEO-TDDFT tutorial!
Solutions
You find all the solution to all the above exercises by clicking on the button below.