Andrew WildmanElectron dynamics can be used to simulate a variety of interesting chemical processes. Here we will use a couple of examples to get spectra from some chemical systems. While spectra are often easier to extract using response methods, this serves as an example to get you acquainted with electronic dynamic simulations in ChronusQ and a starting point for simulating more interesting phenomena.
Sodium D line splitting
First we'll consider an atomic system that can showcase relativistic effects. The P manifold of the neutral sodium atom splits into ^1P_{\frac{3}{2}} and ^1P_\frac{1}{2} states, separated by a few meV. The absorption spectrum of a system can be expressed as:1
{\boldsymbol \sigma}(\omega) \propto \omega \text{Im}\left[ \text{Tr} \left( {\boldsymbol \alpha}(\omega) \right) \right]The dynamic polarizability, \boldsymbol{\alpha}(\omega) can be extracted as the dipole impulse response to an electric dipole perturbation. In general, to get a spectrum, one must perform three separate electronic dynamics simulations, Fourier transforming the dipole signal after an impulse perturbation to get the impulse response function.1 For this system, since it is a single atom, we will assume that the dynamic polarizability is isotropic, so we will only need a single simulation.
Input
The input we'll use is as follows:
#
# Na doublet electron dynamics
#
[Molecule]
charge = 0
mult = 2
geom:
Na 0.0 0.0 0.0
[QM]
reference = X2CBLYP
job = RT
[SCF]
MaxIter = 400
[BASIS]
basis = 6-311+g(d,p)
[INTS]
Alg = InCore
[RT]
TMAX = 12000.
DELTAT = 0.05
FIELD:
StepField(0.,0.0001) Electric 0. 0.005 0.
[MISC]
nsmp = 4
mem = 2GBBecause the splitting we are looking for is very small, a long simulation is required to resolve the difference between the states. Here we choose 12,000 au. (~290 fs) Our "impulse perturbation" is a static field that is on for only a single time step, and we are perturbing in the y Cartesian direction, although any direction should give the same results. Since we're doing a long simulation on a small system, we'll specify that we want the ERIs stored directly in memory (INTS.ALG = INCORE) so we don't spend time every step forming the ERIs. We'll also speed up the calculation with multiple cores, but these computational resources should still be small enough to be done on a desktop computer. Finally, we're giving the SCF some extra time to converge the guess, since this is an atomic system and can only be given a poor guess. (core Hamiltonian guess)
Output
For the output of a real time electron dynamics job, after the SCF results, the real time dynamics options are printed:
*** Parsing RT options ***
Defaulting to MMUT integration algorithm
Defaulting to Magnus 2 restart for MMUT
================================================================================
Real-Time Propagation Settings:
* Simulation Parameters:
Simulation Time: 12000.0000000 ħ / Eh
Number of Steps: 240000
Step Size: 0.0500000 ħ / Eh
0.0012094 fs
* Integration Parameters:
Electronic Integration: Modified Midpoint Unitary Transformation (MMUT)
Restarting MMUT every 50 steps with a(n) Explicit 2nd Order Magnus step
* Perturbation:
Field 1: Electric Dipole Field
* Amplitude (AU) { 0.0000000, 0.0050000, 0.0000000 }
* Step Field TON = 0.0000000 TOFF = 0.0001000
* Misc Parameters:
Matrix Exponential Method: Eigen DecompositionHere we can confirm all of the options we set above. Next comes the time-dependent data:
================================================================================
Time (AU) Energy (Eh) Dipole (X) Dipole (Y) Dipole (Z)
--------------------------------------------------------------------------------
*** Restarting MMUT ***
0.0000 -162.4724136286 0.00000044 -0.00000195 -0.00000110
*** Restarting MMUT ***
0.0500 -162.4724134985 0.00000043 -0.00002602 -0.00000109
*** Restarting MMUT ***
0.1000 -162.4724134985 0.00000043 -0.00006518 -0.00000109
0.1500 -162.4724134985 0.00000043 -0.00010242 -0.00000109
0.2000 -162.4724134985 0.00000043 -0.00014467 -0.00000109
0.2500 -162.4724134985 0.00000043 -0.00017484 -0.00000109
0.3000 -162.4724134986 0.00000043 -0.00020060 -0.00000109
0.3500 -162.4724134983 0.00000043 -0.00022986 -0.00000109
0.4000 -162.4724134988 0.00000043 -0.00024640 -0.00000109
0.4500 -162.4724134984 0.00000044 -0.00026020 -0.00000109
0.5000 -162.4724134984 0.00000044 -0.00027781 -0.00000109
0.5500 -162.4724134987 0.00000044 -0.00028612 -0.00000109
0.6000 -162.4724134985 0.00000044 -0.00029473 -0.00000108
0.6500 -162.4724134984 0.00000044 -0.00031008 -0.00000108
0.7000 -162.4724134988 0.00000044 -0.00031880 -0.00000108
0.7500 -162.4724134984 0.00000044 -0.00033168 -0.00000108
0.8000 -162.4724134985 0.00000044 -0.00035183 -0.00000108
0.8500 -162.4724134986 0.00000044 -0.00036656 -0.00000108
0.9000 -162.4724134987 0.00000044 -0.00038622 -0.00000108
0.9500 -162.4724134984 0.00000044 -0.00041144 -0.00000108
1.0000 -162.4724134986 0.00000044 -0.00042893 -0.00000108We can see that the time-dependent energy and electric dipole moment of the system are printed. While this is useful for ensuring that the simulation is proceeding as you expected, this data is saved to full precision (and with many more properties) to the binary file.
The first thing to notice is that the energy jumps from step 0 to step 1. This is expected - we have a time dependent external field, so molecular energy should not be conserved. After that, some energy fluctuations are observed, but they are rather small. By checking the end of the simulation:
11999.5000 -162.4724135131 0.00001433 -0.00110638 -0.00014043
11999.5500 -162.4724135131 0.00002048 -0.00111576 -0.00013357
11999.6000 -162.4724135131 0.00002533 -0.00112343 -0.00012923
11999.6500 -162.4724135131 0.00002703 -0.00112735 -0.00012679
11999.7000 -162.4724135131 0.00002838 -0.00113122 -0.00012445
11999.7500 -162.4724135131 0.00002945 -0.00113498 -0.00012469
11999.8000 -162.4724135131 0.00002923 -0.00113705 -0.00012638
11999.8500 -162.4724135131 0.00003019 -0.00114101 -0.00012710
11999.9000 -162.4724135132 0.00003193 -0.00114640 -0.00012913
*** Restarting MMUT ***
11999.9500 -162.4724135131 0.00003380 -0.00115155 -0.00013125
*** Restarting MMUT ***
*** Saving data to binary file ***
12000.0000 -162.4724135131 0.00003782 -0.00115962 -0.00013090we see that the total energy drift over the 12,000 au is 1.47e^{-8}\ E_h, which certainly suffices to ensure that the response is physically relevant.
Now we can use the data on the binary file to ensure we get the full precision on the dipole signal, extract it, use the Padé approximant to the Fourier transform,2 and plot the spectral region relevant to the sodium d line, and we get the spectrum below:
Water spectrum
In this example, we will use a molecular system that does require 3 different simulations to calculate the spectrum: gas phase water.
Input
The first input we'll use has a perturbation along the x Cartesian axis.
#
# Water B3LYP spectrum
#
[Molecule]
charge = 0
mult = 1
geom:
O 0. -0.07579184359 0.
H 0.866811829 0.6014357793 0.
H -0.866811829 0.6014357793 0.
[QM]
Reference = RB3LYP
Job = RT
[BASIS]
Basis = cc-PVDZ
[INTS]
Alg = InCore
[RT]
TMax = 620.
DeltaT = 0.05
Field:
StepField(0.,0.0001) Electric 0.0001 0. 0.
[MISC]
NSMP = 4
Mem = 2GBWe use many of the techniques above to speed up the calculation. Here, we need a significantly shorter time to resolve the spectrum. We'll use 620 au. (~15 fs)
The other two simulations just change the axis along which the system is perturbed. For y:
Field:
StepField(0.,0.0001) Electric 0. 0.0001 0.and for z:
Field:
StepField(0.,0.0001) Electric 0. 0. 0.0001Output
The same output sections exist for this simulation as above. We see that the energy doesn't change over the course of the simulation, (other than the initial step) so we are ready to extract the spectrum. Extracting and damping the signal with a damping factor of 0.001 au, then transforming to frequency space with the Padé approximant to the Fourier transform,2 and finally tracing over the diagonal of the dynamic polarizability, we get the following spectrum:
References
-
Goings, J. J., Lestrange, P. J., & Li, X. (2018). Real‐time time‐dependent electronic structure theory. Wiley Interdisciplinary Reviews: Computational Molecular Science, 8(1), e1341.
↩ ↩ 2 -
Bruner, A., LaMaster, D., & Lopata, K. (2016). Accelerated broadband spectra using transition dipole decomposition and Padé approximants. Journal of Chemical Theory and Computation, 12(8), 3741-3750.
↩ ↩ 2