Density Functional Theory Formalism and Applications Mauro Stener

Density Functional Theory Formalism and Applications Mauro Stener Dipartimento di Scienze Chimiche e Farmaceutiche Università degli Studi di Trieste Via L. Giorgieri 1, 34127 TRIESTE E-mail: stener@units. it XIII-EM-TCCM-European Master in Theoretical Chemistry and Computational Modelling 13 th International Intensive Course, Perugia (I), 3 rd – 28 th September 2018

Summary 1. THEORY - DFT Review of DFT basic theory and equations 2. THEORY - TDDFT with the aim to underline what is general from Linear Response formalism. 3. COMPUTATIONAL - Some numerical aspects of the implementation in the ADF code, most important computational choices, how to assign spectral features. New TDDFT algorithm. 4. APPLICATIONS: 4. a. Plasmons in metal clusters. 4. b. Core electron excitations (NEXAFS).

ADF program: a) DFT Kohn-Sham (KS) Equations 1. LCAO formulation (STO basis set) 2. Numerical integrals 3. Density fitting Auxiliary basis functions

ADF program: available STO basis set STO (Slater Type Orbital) : more physical than GTO (Gaussian) ADF employs STO for both basis set for wavefunction and auxiliary basis set for density fitting STO (Slater Type Orbital) : more physical than GTO (Gaussian) SZ: Single Zeta (Minimal Basis) 1 function for each shell (1 s, 2 p ecc. ) DZ: Double Zeta, 2 functions for each shell DZP: DZ + polarization: higher angular momentum (es. 3 d for C) TZ Triple Zeta, TZP Triple Zeta + polarization ET: Even Tempered (very large basis set, dangerous for numerical stability) ZORA: relativistic basis set

ADF program: b) TDDFT Equations Diagonalization of W matrix furnishes discrete excitation energies and intensities Gaussian broadening TDDFT results, lowest n eigenvalues

TDDFT equations: i and j run over Nocc a and b run over Nvirt TDDFT results: FI eingenvectors give 1) Intensity (Oscillator Strength f) 2) Allow analysis of excited state in terms of 1 h 1 p excited configurations: i

M. Stener, A. Nardelli, R. De Francesco and G. Fronzoni, JPCC 111 (2007) 11862 Interband sp d 1 transition many excited config.

Efficient numerical algorithms (ADF) to solve: Computational effort: storage CPU i and j run over Nocc a and b run over Nvirt ØDavidson iterative diagonalization: n lowest eigenvalues ØEfficient management of matrix-vector Pout=K • Pin product K matrix is not calculated and is not stored! N 4 scaling!

Density fitting procedure: fn Auxiliary basis functions for the density an: fitting coefficients such that:

Numerical integration (Gaussian weights and points) Analytical! A favourable N 3 scaling is gained, run over three index: i a n

TDDFT spectra and numerical stability issues ØDavidson iterative diagonalization: n roots (lowest eigenvalues) ØStable up to n = 350, for n = 400 numerical instability shows up ØNumerical stability: improved by tight SCF convergence criterion: |FP-PF|<10 -8

Au 64+: TEST (BASIS SET) • negligible difference DZ TDDFT / DZ geo opt • computational economy TZP TDDFT / TZP geo opt DZ TDDFT / TZP geo opt with DZ optimized geometry Computational details: TZP TDDFT / DZ geo opt • 100 roots • TDDFT LB 94 XC potential • geo opt VWN XC potential • All Electron basis set

Au 64+: TEST (FROZEN CORE) • negligible differences • Computational saving around 53% from AE to FC 4 f Core orbital frozen up to Au 4 f Computational details: • 100 roots • LB 94 XC potential • TZP basis set

Au 64+: TEST (EXCHANGE CORRELATION POTENTIAL) • SAOP is usually considered more accurate than LB 94 • negligible differences • LB 94 is chosen, since SAOP does not support FC! LB 94 XC potential Computational details: • 100 roots • AE TZP basis set

Relativistic effects in Au compounds 6 s shell Large relativistic contraction of the Au 6 s shell P. Pyykko and J. P. Desclaux, Acc. Chem. Res. 12 (1979) 276 Strong relativistic narrowing of the 5 d – 6 s gap J. P. Desclaux and P. Pyykko, Chem. Phys. Lett. 39 (1976) 300

Relativistic effects: transformation • in ADF: ZORA (Zero Order Regular Approximation) • ZORA: well behaved over the nuclei • Two components: Spin-Orbit (SO) coupling included • If SO is neglected: Scalar Relativistic (SR)

Gold bimetallic clusters: M@Au 12 • Icosaedral bimetallic gold clusters: Au cage with encapsulated heteroatom WAu 12 Mo. Au 12 • First theoretically predicted, then synthesized and characterized by spectroscopy • Analysis of the spin orbit coupling on optical spectra

WAu 12: scalar relativistic electronic structure P. Pyykko and N. Runeberg, Angew. Chem. In. Ed. 41 , 2174(2002). M. Stener, A. Nardelli, and G. Fronzoni J. Chem. Phys. 128, 134307 (2008)

Spin-orbit interaction in molecules • Similar to atoms: lower degeneracy • States classified according to Double Groups I Ih 2 • Example: Ih h Ag E 1 g(1/2) T 1 g E 1 g(1/2) + Gg(3/2) T 2 g Ig(5/2) Gg E 2 g(7/2) + Ig(5/2) Hg Gg(3/2) + Ig(5/2) Au E 1 u(1/2) T 1 u E 1 u(1/2) + Gu(3/2) T 2 u Iu(5/2) Gu E 2 u(7/2) + Iu(5/2) Hu Gu(3/2) + Iu(5/2)

WAu 12: spin-orbit electronic structure Exp: photodetachment of WAu 12 -

TDDFT electronic excitations (valence) WAu 12 SR ZORA TZ 2 P LB 94 Excitation energy (e. V)

WAu 12: spin-orbit effect over spectrum SR SO low energy • red shift • intensity reduction SO SO above 5 e. V: strong configuration remixing

WAu 12: Scalar Relativistic vs Spin-Orbit TDDFT 1. Discrete lines are very sensitive to SO 2. Broadened spectrum is much less sensitive to SO 3. When “density” of lines is high SO can be reasonably neglected 4. SO inclusion is computationally very demanding

A New TDDFT Algorithm for Large Systems and Plasmons in Metal Clusters 1. Conventional (Casida like) TDDFT for large systems (bottleneck: diagonalization of large matrix with Davidson algorithm) 2. New TDDFT algorithm: inspired by previous work on photoionization to overcome problem 1) 3. Some applications to metal clusters 4. Just published: O. Baseggio, G. Fronzoni and M. Stener, J. Chem. Phys. , 143 (2015) 024106

The new algorithm to solve TDDFT equations Before to describe the new algorithm, it is worth mentioning other recent TDDFT strategies, which are very promising for applications to large systems: • the explicit time-propagation technique over real space grids, implemented in the OCTOPUS program [K. Yabana and G. F. Bertsch, Phys. Rev. B 54, 4484 (1996), M. A. L. Marques, A. Castro, G. F. Bertsch and A. Rubio, Comput. Phys. Comm. , 151 (2003) 60. ] • the superoperator formulation of the TDDFT, which allows the calculation of the dynamical polarizability by means of a very efficient Lanczos method, implemented with plane waves basis set [B. Walker, A. M. Saitta, R. Gebauer and S. Baroni, Phys. Rev. Lett. 96 (2006) 113001. ] • the very recent simplified TDA scheme of S. Grimme [J. Chem. Phys. 138 (2013) 244104] • the Linear-scaling TDDFT of Zuehlsdorff [T. J. Zuehlsdorff, N. D. M. Hine, J. S. Spencer, N. M. Harrison, D. J. Riley, and P. D. Haynes, J. Chem. Phys. 139 (2013) 064104] • TDDFT time-propagation with Transition Contribution Map (TCM) of Hakkinen [Sami Malola, Lauri Lehtovaara, Jussi Enkovaara and Hannu Häkkinen, ACS Nano 7 (2013) 10263]

Exploit linearity of the problem: defines the kernel K defines the susceptibility KS The Response Equation becomes: To solve : represent the response equation in the auxiliary density fitting basis functions, dimension of fitting set (density-based formulation), at variance with density-matrix Casida formulation (1 h 1 p space) M. Stener, G. Fronzoni and P. Decleva, J. Chem. Phys. , 122 234301(1 -11) (2005).

KS is known for the non-interacting (KS) system: Infinite summations are involved: not practical for numerical implementation as it is.

Extract the spectrum from polarizability: The first order density (1) calculated from TDDFT equations, with = R+i I basis set: auxiliary density fitting functions! The “only” problem is to calculate efficiently KS!

Change the double sum: The energy grid consists of P+1 knots: and are defined P intervals: Now it is possible to change the double sum of previous equation (13) as follows:

Build the matrix: TDDFT equation: With this scheme, the M matrix is calculated at each photon energy , simply as a linear combination of matrices Gk with coefficients sk( ). The calculation of KS is fast (just sum of matrices) and small (over the fitting functions)! Approximations: imaginary broadening and energy discretization.

Technical details: Pair fitting (very efficient!) New integrals, easily implemented in ADF (numerical integration over grid)

Technical details: Once matrix M S and vector d are available, the TDDFT equation is solved (linear system) and the vector b is obtained, then polarizability and absorption spectrum are calculated.

Accurate polarizability and analysis It is useful to employ the Modified Sternheimer Approach (MSA) for the response: These equations are convenient to be represented in the basis of virtual KS orbitals:

Accurate polarizability and analysis So actually the polarizability (and the spectrum) are calculated by: With: In practice: - from the imaginary part of the absorption spectrum is obtained - from imaginary part of Pia the analysis in terms of one-electron excited configurations and Transition Contribution Map (TCM) is obtained [Sami Malola, Lauri Lehtovaara, Jussi Enkovaara and Hannu Häkkinen, ACS Nano 7 10263 (2013). ]

Validation of the method, comparison with ADF H 2 molecule, electric field polarization along the bond (Z) DZ basis set, LB 94 XC, the discrete lines of ADF convoluted by Lorentzian functions with Half Width Half Maximum (HWHM) = i= 0. 3 e. V, original fitting: 39 functions! ADF: Complex : Absorption spectrum Real dynamical polarizability ADF

Validation of the method, Na 2 molecule, E polarization along (Z) and perpendicular (X) to the bond. 2. 35 e. V 4 g 4 u DZP basis set, LB 94 XC, Lorentzian (HWHM) = i= 0. 15 e. V Fitting: original set (14 s 7 p 6 d 2 f 1 g) for Na (88 functions), reduced to (7 s 5 p) only 22 functions! Z 3. 38 e. V 4 g 2 u (4 g)2(4 u)0(2 u)0 X

Validation of the method, Au 2 molecule, E polarization along (Z) and perpendicular (X) to the bond. Scalar Relativistic ZORA, DZ basis set, LB 94 XC, Lorentzian (HWHM) = i= 0. 30 e. V Fitting: original set (25 s 15 p 12 d 8 f 7 g ) for Au (249 functions), reduced to (5 s 4 p 4 d 1 f ) only 44 functions! Max. error: 0. 2 e. V

Comments: 1. Technical: ADF employs STO, FC DZ is adequate, density fitting is efficient, TDDFT diagonalization is optimized. 2. Computational: ADF code is fully parallelized via standard MPI. 3. Careful choices: numerical parameters (SCF convergence and number of roots to extract in TDDFT) and proper relativistic scheme for gold or other heavy elements. 4. Outlook: use the new TDDFT algorithm to enlarge the energy range without numerical stability problems