Tutorial__GS__T07_GS_TDDFT - fmauger1/QMol-grid GitHub Wiki
Tutorial 6: Time-dependent density-functional theory
This tutorial walks through the process of setting up and simulating time-dependent density-functional theory (TDDFT) dynamics in an open-shell one-dimensional molecular model with 3 atomic centers and 5 active electrons.
Table of Contents
Introduction
For a given set of initial Kohn-Sham orbitals, $\lbrace \phi_k \rbrace_k$ , the TDDFT dynamics is described by the nonlinear system of partial differential equations, in atomic units
$$ i\partial_t \phi_k (x;t)=\hat{\mathcal{H}} _\mathrm{DFT} [ \lbrace \phi_k \rbrace_k ;t ] (x;t)~\phi_k (x;t).~~~~~~(1) $$
In the adiabatic (a.k.a. local in time) approximation $\hat{\mathcal{H}} _\mathrm{DFT}$ is the same DFT Hamiltonian operator as for ground-state calculations of tutorial 5, except for possible added external driving fields (responsible for the explicit time dependency). The QMol-grid package relies on the canonical Hamiltonian structure of the TDDFT [Mauger 2024] to integrate the dynamics of equation (1).
In this tutorial, we illustrate how to use the QMol-grid package to integrate the TDDFT dynamics of an open-shell one-dimentional molecular ion model with 3 atomic centers and 5 active electrons.
TDDFT simulation
Initial condition
In the QMol-grid package, TDDFT simulations are decoupled from setting up the initial condition, which must be done independently. For our example, we start by calculating the ground state of the same molecular system as in tutorial 5.
% Molecular model
V_1 = QMol_Va_softCoulomb('atom','X_1','charge',2,'position',-3);
V_2 = QMol_Va_softCoulomb('atom','X_2','charge',2,'position', 0);
V_3 = QMol_Va_softCoulomb('atom','X_3','charge',2,'position', 3);
% DFT model
Vext = QMol_DFT_Vext('atom',{V_1,V_2,V_3});
Vh = QMol_DFT_Vh_conv;
Vxc = {QMol_DFT_Vx_LDA_soft,QMol_DFT_Vc_LDA_soft};
DFT = QMol_DFT_spinPol( ...
'xspan', -50:.1:50, ...
'occupation', {[1 1 1],[1 1 1]}, ...
'externalPotential', Vext, ...
'HartreePotential', Vh, ...
'exchangeCorrelationPotential', Vxc, ...
'selfInteractionCorrection', 'ADSIC' );
% DFT ground state
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(DFT);
yielding (removing the license and funding disclaimers)
(...)
=== External components ==================================================
* convertUnit
V-01.02.000 (09/30/2022) F. Mauger
* fourierTool
V-01.02.002 (07/18/2023) F. Mauger
=== Build density-functional-theory (DFT) model ==========================
* Discretization OK
* Kohn-Sham orbitals OK
* External potential OK
* Hartree potential OK
* Exchange-correlation potential OK
=== Self-consistent-field (SCF) methods ==================================
* Eigen-state solver for DFT Hamiltonians MATLAB eigs function
Tolerance = 1e-12
Max. iter. = 300
Basis dim. = 100
V-01.21.000 (06/17/2024) F. Mauger
* Self-consistent-field (SCF) Anderson's mixing
DFT-SCF solver using an Anderson's mixing scheme [Anderson 1965], as
described in [Johnson 1988].
Tolerance = 1e-10
Max. iter. = 100
Mix. coeff. = 0.5
Mix. mode = density
Conv. test = density
V-01.21.000 (06/17/2024) F. Mauger
=== References ===========================================================
[Anderson 1965] D.G. Anderson, "Iterative procedures for nonlinear
integral equations," Journal of the ACM 12, 547 (1965).
[Johnson 1988] D. D. Johnson, "Modified Broyden's method for
accelerating convergence in self-consistent calculations," Physical
Review B 38, 12807 (1988).
[Mauger 2024b] F. Mauger and C. Chandre, "QMol-grid: A MATLAB package
for quantum-mechanical simulations in atomic and molecular systems,"
arXiv:2406.17938 (2024).
(...)
=== SCF iterations =======================================================
Iter. Orbital energies (-eV) Error
----- ----------------------------------------------------- --------
0 24.01 | 20.68 | 16.28 8.349e-02
24.01 | 20.68 | 16.28 8.349e-02
1 24.73 | 20.87 | 16.51 4.601e-03
24.73 | 20.87 | 16.51 4.601e-03
2 24.71 | 20.86 | 16.50 1.740e-03
24.71 | 20.86 | 16.50 1.740e-03
3 24.71 | 20.86 | 16.50 4.762e-04
24.71 | 20.86 | 16.50 4.762e-04
(...)
17 24.71 | 20.86 | 16.50 1.538e-10
24.71 | 20.86 | 16.50 1.538e-10
18 24.71 | 20.86 | 16.50 5.892e-11
24.71 | 20.86 | 16.50 5.892e-11
19 24.71 | 20.86 | 16.50 2.336e-11
24.71 | 20.86 | 16.50 2.336e-11
------------------------------------------------------------------------
Kohn-Sham orbitals converged to tolerance
=== Orbital energies =====================================================
Orbital Occ. (elec.) Energy (-eV) Error(a.u.)
------- ------------ ------------ -----------
1 1.00 24.713 2.624e-12
2 1.00 20.859 2.116e-12
3 1.00 16.502 2.735e-12
------- ------------ ------------ -----------
1 1.00 24.713 2.642e-12
2 1.00 20.859 2.110e-12
3 1.00 16.502 2.703e-12
------------------------------------------------------------------------
=== DFT-component energies ===============================================
Component Energy (a.u.) Energy (eV) [ spin up/down (eV) ]
----------- ------------- ------------- ---------------------
Kinetic 1.058 28.801 [ 14.40/ 14.40 ]
External -16.466 -448.059 [ -224.03/ -224.03 ]
Hartree 5.583 151.923
Exch.-corr. -0.564 -15.358
----------- ------------- -------------
Total -10.389 -282.692
----------------------------------------------
##########################################################################
A quick glance at the error estimates confirms that the ground-state calculation converged. Also, we see that the up- and down-spin molecular orbitals (the two groups 1-3 of orbitals) have the same energies, as expected for our neutral molecular model, which has a closed-shell electronic structure.
Next, we manually induce an excitation in the molecular cation by successively (i) replacing one of the Kohn-sham orbitals by a superposition of molecular-orbital states (excitation part) and (ii) removing an electron, going from 3 to 2, from the down-spin Kohn-Sham orbitals (ionization part).
% Induce excitation
DFT.orbital.set('orbitalDown',[DFT.KSO.KSOdw(:,1) (DFT.KSO.KSOdw(:,2)+DFT.KSO.KSOdw(:,3))/sqrt(2)]);
% Induce ionization
DFT.set('occupation',{[1 1 1],[1 1]});
We now have a non-stationary set of Kohn-Sham orbitals, leading to field-free dynamics under the flow of equation (1).
TDDFT simulation
With the DFT molecular model, including the initial condition, in hand we now move to integrating the subsequent field-free TDDFT dynamics. For this, we select a fourth-order Forest Ruth [Forest 1990] symplectic split-operator scheme [Mauger 2024].
TDDFT = QMol_TDDFT_SSO_4FR( ...
'time', 0:10:100, ...
'timeStep', 2e-2, ...
'saveDensity', true, ...
'saveDensityTime', 1);
In our example, the TDDFT object is created with:
- The first pair of arguments specifies that the integration should start at time t=0 and end at t=100 a.u. The step of 10 a.u., is unrelated to the propagation time step (next) and instead specifies the time intervals to use in progress display (see also the next tutorial for more details about input parameters).
- The second pair of arguments specifies the (fixed) time step for the propagation.
- The third pair of arguments indicates that the one-body density should be saved periodically, with the period specified by the fourth pair of arguments, i.e., every 1 a.u. in our case.
Then, we launch the TDDFT integration with
TDDFT.propagate(DFT);
yielding (again removing the license and funding disclaimers)
(...)
=== External components ==================================================
* convertUnit
V-01.02.000 (09/30/2022) F. Mauger
* fourierTool
V-01.02.002 (07/18/2023) F. Mauger
=== Build density-functional-theory (DFT) model ==========================
* Discretization OK
* Kohn-Sham orbitals OK
* External potential OK
* Hartree potential OK
* Exchange-correlation potential OK
=== Time propagation =====================================================
* The TDDFT model dynamics matches the canonical Hamiltonian formalism
discussed in [Mauger 2024].
* Symplectic split-operator propagation scheme [Mauger 2024],
with a VTV split motif.
V-01.21.000 (06/17/2024) F. Mauger
* Field free propagation
* Forest-Ruth propagation scheme [Forest 1990] 4th order
V-01.21.000 (06/17/2024) F. Mauger
* Parameters
Time interval = 0.00 to 100.00 a.u. (0.00 to 2.42 fs)
Time step = 0.020 a.u. (0.484 as)
* Restart option is disabled
=== Output results =======================================================
* Kohn-Sham orbitals and one-body density:
> Save the one-body density
=== References ===========================================================
[Forest 1990] E. Forest and R.D. Ruth, "Fourth-order symplectic
integration," Physica D: Nonlinear Phenomena 43, 105 (1990).
[Mauger 2024] F. Mauger, C. Chandre, M.B. Gaarde, K. Lopata, and K.J.
Schafer, "Hamiltonian formulation and symplectic split-operator
schemes for time-dependent density-functional-theory equations of
electron dynamics in molecules," Communications in Nonlinear Science
and Numerical Simulation 129, 107685 (2024).
[Mauger 2024b] F. Mauger and C. Chandre, "QMol-grid: A MATLAB package
for quantum-mechanical simulations in atomic and molecular systems,"
arXiv:2406.17938 (2024).
(...)
=== Time propagation =====================================================
Iteration time (a.u.) Iteration Iteration progress (%)
start finish number 0 20 40 60 80 100
--------- --------- --------- -------------------------------
0.000 10.000 1 |||||||||||||||||||||||||||||||
10.000 20.000 501 |||||||||||||||||||||||||||||||
20.000 30.000 1001 |||||||||||||||||||||||||||||||
30.000 40.000 1501 |||||||||||||||||||||||||||||||
40.000 50.000 2001 |||||||||||||||||||||||||||||||
50.000 60.000 2501 |||||||||||||||||||||||||||||||
60.000 70.000 3001 |||||||||||||||||||||||||||||||
70.000 80.000 3501 |||||||||||||||||||||||||||||||
80.000 90.000 4001 |||||||||||||||||||||||||||||||
90.000 100.000 4501 |||||||||||||||||||||||||||||||
--------- --------- --------- -------------------------------
Time propagation finished without error. All results are saved in the
TDDFT object and specified output MATLAB file(s).
##########################################################################
The "Time propagation" block summarizes the propagator options for the numerical integration of equation (1). The following "Output results" block lists the observables that are calculated and stored alongside the propagation. Here it confirms that only the one-body density is retained.
At the end of the simulation, the DFT object has been updated to contain the Kohn-Sham orbitals at t=100 a.u. The time-dependent one-body density is stored in the TDDFT object itself (see next).
Plotting the result
To conclude this tutorial, we briefly illustrate how to recover calculated observables out of the TDDFT object. Each set of observable is stored in a separate structure property in the TDDFT object, which contains (i) the exact time vector at which the quantity has been saved and (ii) the observable itself. In our case, the structure of interest is TDDFT.outDensity with the up- and down-spin densities respectively stored in the fields totalUp and totalDown. The densities are matrices with columns corresponding to the successive saved times.
To plot the spin density, defined as the difference between the up- and down-spin one-body densities, we use
figure
imagesc(TDDFT.outDensity.time,DFT.xspan,TDDFT.outDensity.totalUp-TDDFT.outDensity.totalDown)
set(gca,'box','on','FontSize',12,'LineWidth',2,'YDir','normal')
xlim(TDDFT.outDensity.time([1 end]))
ylim([-10 10])
xlabel('time (a.u.)')
ylabel('position (a.u.)')
title('spin density')
colorbar vert
producing (note that the ground-state calculation start from a random seed and thus the resulting Kohn-Sham molecular orbitals are defined with an arbitrary sign that can change from calculation to calculation, this means that the resulting figure may be flipped vertically)
We further explain how to interface input parameters and output results from calculations in the next tutorial.
References
[Forest 1990] E. Forest and R.D. Ruth, "Fourth-order symplectic integration," Physica D: Nonlinear Phenomena 43, 105 (1990).
[Mauger 2024] F. Mauger, C. Chandre, M.B. Gaarde, K. Lopata, and K.J. Schafer, "Hamiltonian formulation and symplectic split-operator schemes for time-dependent density-functional-theory equations of electron dynamics in molecules," Communications in Nonlinear Science and Numerical Simulation 129, 107685 (2024).
Notes
The results displayed in this tutorial were generated using version 01.21 of the QMol-grid package.