Real-Time Time-Dependent Density Functional Theory Implementation of Electronic Circular Dichroism Applied to Nanoscale Metal-Organic Clusters

Electronic circular dichroism (ECD) is a powerful spectroscopical method for investigating chiral properties at the molecular level. ECD calculations with the commonly used linear-response time-dependent density functional theory (LR-TDDFT) framework can be prohibitively costly for large systems. To alleviate this problem, we present here an ECD implementation for the projector augmented-wave method in the real-time-propagation TDDFT (RT-TDDFT) framework in the open-source GPAW code. Our implementation supports both local atomic basis set and real-space finite-difference representations of wave functions. We benchmark our implementation against an existing LR-TDDFT implementation in GPAW for small chiral molecules. We then demonstrate the efficiency of our local atomic basis set implementation for a large hybrid nanocluster.


I. INTRODUCTION
Chirality is an essential property in several branches of science and technology. Chiral molecules play a fundamental role in biological activities; 1 for example, DNA double-helices are right-handed and amino acids are left-handed. Chirality is also critical in pharmaceuticals. For example, R-enantiomer thalidomide is effective against morning sickness for pregnant women, but the S-species produce fetal deformations. 2,3 In addition, chiral molecules and chiral nanomaterials have many potential applications in catalysis, sensors, spintronics, optoelectronics, and nanoelectronics. [4][5][6][7][8][9][10][11][12][13] Therefore, the determination of the handedness of chiral systems is of paramount importance. Chiral molecules absorb left and right circular polarizations of light differently. This difference is probed in electronic circular dichroism (ECD) spectroscopy. ECD is defined as ∆ = L − R , where L and R are the molar extinction coefficients for left and right circularly polarized light, respectively. Since ECD is highly sensitive to small details in the atomic structure of molecules and unique for each conformation, it is a powerful technique for characterizing chiral systems and for distinguishing enantiomers. 14 The ECD spectroscopy accompanied by computational modeling can help to get insightful knowledge of the atomic structures of chiral biomolecules and nanoclusters. [15][16][17] For example, in our earlier work, we have identified by comparing simulated and measured ECD that the Ag + -mediated guanine duplex has the lefthanded helix configuration, 16 while the Ag + -mediated cytosine has the right-handed helix configuration. 17 Most computational ECD approaches are based on the linear-response formalism. 18 Time-dependent densityfunctional theory (TDDFT) 19 has become the linear-response method of choice due to its favorable balance of accuracy and computational cost, compared to quantum chemical approaches such as coupled cluster and configuration interaction methods. 18,20,21 In linear-response TDDFT (LR-TDDFT), the Casida equation 22,23 is solved in the basis of Kohn-Sham (KS) particle-hole transitions in the frequency domain. 24,25 A full ECD spectrum requires the calculation of a large number of transitions from occupied to unoccupied states.This becomes computationally prohibitive for large systems with a high density of states, resulting in an unfavorable O(N 5 ) scaling, where N is the system size.
An alternative to LR-TDDFT is real-time-propagation time-dependent density-functional theory (RT-TDDFT). In RT-TDDFT, the system is subjected to an initial perturbation and the KS wave functions are propagated in the time domain by numerically integrating the timedependent KS equations. The real-time approach captures the same information as LR-TDDFT, for small initial perturbations, and incorporates nonlinear spectral information for larger initial perturbations. 23,[26][27][28] RT-TDDFT scales as O(N 2 ), but suffers from a large prefactor. LR-TDDFT is usually faster for small systems such as small organic molecules. For large molecules, clusters, and nanoparticles, RT-TDDFT becomes more cost-effective than LR-TDDFT. 23,29 RT-TDDFT computation of ECD has been implemented for a variety of basis sets: real-space grids, 27,30 Gaussian-type atomic orbitals, 29 and a mix of Gaussiantype and plane-wave basis sets. 31 In this work, we present a RT-TDDFT ECD implementation in the open-source GPAW package. 32,33 Our implementation uses the projector augmented wave (PAW) method 34 and supports both localized basis sets (LCAO mode) 35 and real-space grids (grid mode). 32,33 We verify our implementation by comparing our results to those calculated with the existing LR-TDDFT implementation in GPAW. We also benchmark the LCAO mode against accurate real-space grid calculations.
In GPAW, the time-dependent density and potential are expressed on a uniform grid, and the matrix elements of the potential are evaluated on this grid. 32 The smoothness of these quantities allows for a coarse grid spacing. The LCAO-PAW pseudo wavefunctions can form a local and efficient representation suitable for systems with hundreds of atoms. 35 Previous work has shown that LCAO RT-TDDFT in GPAW is capable of simulating the optical spectrum of a silver cluster of more than 500 atoms (Ag 561 ). 36,37 In this work, we demonstrate the efficiency of our LCAO RT-TDDFT ECD implementation for a large ligand-protected Ag 78 cluster 38 consisting of over 1000 atoms and over 4000 electrons.
The rest of this paper is organized as follows. In the Methods section we illustrate our implemented methodologies by showing how the ECD is calculated from timedependent magnetic dipole moment and how the timedependent magnetic dipole moment is calculated with RT-TDDFT in GPAW. The important information considering all simulations is described in the Computational Methods section. In the Results section we demonstrate the capability of our implementations to predict ECD spectrum for four test cases. Finally, we give a summary in the Conclusion section.

II. METHODS
In RT-TDDFT, the KS wave functions are propagated in time in response to a time-dependent potential starting from an initial state, here chosen to be the ground state. The time-dependent KS equation is defined as where H KS (t) is the time-dependent KS Hamiltonian and ψ n (r, t) is a time-dependent KS wave function. A common practice in time-propagation schemes is to use the weak δ-kick approach 26 to calculate the linear-response functions. After perturbing H KS (t) by the δ-kick at t = 0, Eq. (1) is propagated using the semi-implicit Crank-Nicolson method, whose numerical reliability has been demonstrated previously. 36 In this work, we implement and thoroughly benchmark the calculation of time-dependent magnetic moment within the time-propagation framework for obtaining the ECD spectrum. In the following, we derive the relevant equations within the PAW method. 34 This derivation partly follows the one shown by Varsano et al. 30 A. ECD from the Induced Time-Dependent Magnetic Moment A commonly used experimental quantity to measure ECD is the difference in molar extinction coefficients, which is given by Here ω is the the energy of the incident light, c the speed of light, the reduced Planck constant, N A Avogadro's constant, and R(ω) cgs the rotatory strength in cgs units. The quantity that characterizes ∆ (ω) and therefore the ECD spectrum is the rotatory strength. The relationship between rotatory strength in cgs units and rotatory strength in atomic units (denoted as R(ω)) is where e is the elementary charge, m e the mass of an electron, and α the fine structure constant. We will work in atomic units and perform the required unit conversions afterwards.
The rotatory strength is defined through the optical rotatory response tensor by where index k enumerates Cartesian coordinates (k ∈ {x, y, z}). Next we will derive the relationship between the optical rotatory response tensor and the timedependent magnetic dipole moment, which is the quantity calculated in TDDFT.
For a system in an external electric field E = [E x , E y , E z ] (no external magnetic field present), the induced time-dependent magnetic dipole moment in direction j, m j (t), has the expansion where β jk is the jk component (indices j, k ∈ {x, y, z}) of the optical rotatory response tensor, and E k is the electric field component. The response tensor β jk describes the induced magnetic dipole moment in the j direction for a perturbing electric field in the k direction. The response is causal, which means that β jk (t) vanishes for negative values of t. In the weak-field limit, the time-dependent magnetic moment is dominated by the first-order term given by the linear-response functions β jk (t). By using the convolution theorem and the properties of Fourier transforms, the first-order term of Eq. (5) can be written in the frequency domain as To resolve all components of β jk (t), we perform the δkick 26 in all three Cartesian directions using a perturbing electric field of the form E (k) (t) = κkδ(t). The k superscript in parenthesis indicates the kick direction, to be distinguished from the component subscript. We keep the intensity κ weak to restrict our calculations to the linear-response regime.
In the frequency domain, the delta kick becomes a constant for all frequencies Eq. (6) then simplifies to Combining Eqs. (4) and (8), we get the equation for rotatory strength expressed with magnetic dipole moment In our methodology, m In principle, the integration interval goes from zero to infinity. In practice, a finite propagation time (T ) suffices by introducing an artificial lifetime ω → ω + i σ 2 2 t, where σ is the parameter that determines the line width of the Gaussian line shape. Introducing this into Eq. (10) gives For a desired value of σ, the propagation time T needs to be large enough so that e − σ 2 2 T 2 ≈ 0.
The magnetic moment is defined by the following operator (in atomic units) The expectation value of the time-dependent magnetic moment is obtained as where f n is the occupation number of the n:th KS state and ψ n (r, t) is the time-evolved KS wave function.
In the PAW method 34 the wave functions ψ n (r, t) have decomposition whereψ n (r, t) is a smooth pseudo wave function and ai φ a i (r) −φ a i (r) p a i |ψ n (t) a local correction inside an atomic augmentation sphere.p a i is a localized projector function and φ a i andφ a i are partial and pseudo partial waves, respectively. These quantities are specific to PAW. In the PAW formalism, the expectation value in Eq. (13) becomes where the augmentation-sphere contribution is For evaluating ∆M a ij , the required matrix elements of the form φ a i |r × ∇|φ a j are evaluated in two atom-centered parts as φ a i |(r − R a ) × ∇|φ a j + R a × φ a i |∇|φ a j , where R a is the coordinate of atom a.
In the LCAO expansion, the time-dependent pseudo wave functionψ n (r, t) is written as a linear combination of atom-centered basis functions ϕ µ (r − R a ) where c µn (t) are the time-dependent expansion coefficients. With this LCAO expansion, Eq. (15) can be written compactly as where ρ µν (t) = n f n c µn (t)c * νn (t) is the KS density matrix in the LCAO basis. The matrix elements M µν are given by the pseudo and augmentation contributions In real-space grid mode the magnetic moment is calculated using equation Eq. (15). In LCAO mode, Eq. (17) is used. The matrix elements M µν are time independent and calculated only once before the time propagation.
After a complete time propagation, the recorded m(t) is transformed to frequency domain as a post processing step according to Eq. (11) at each desired ω value. Then the rotatory strength is calculated according to Eq. (9).
We note that research of gauge origin issues is not within the scope of this work and the reader is suggested to explore a recent detailed investigation on the matter. 31

C. Computational Methods
For our calculations in this work, we used the PBE exchange-correlation functional, 39 unless otherwise mentioned. The molecules were placed into a cubic unit cell with the vacuum size of 8 Å. The real-space grid spacing was chosen as h = 0.2 Å. We tested coarser settings with h = 0.3 Å in RT/LCAO mode for the Ag + -mediated guanine duplex case to demonstrate that such a coarser grid is sufficient for calculating the ECD spectrum within LCAO mode, where the grid is used to represent only real-space density and potential. 35,36 Per atom, the electronic configuration of valence electrons is H(1s 1 ) O(2s 2 2p 4 ), C(2s 2 2p 2 ), N(2s 2 2p 3 ), S(3s 2 3p 4 ), P(3s 2 3p 3 ), F(2s 2 2p 5 ) and Ag(4d 10 5s 1 ). The remaining electrons were treated as a frozen core. The default PAW dataset package 0.9.20000 was used for all the atoms.
In the LCAO mode, the default GPAW double-zeta polarized (dzp) basis sets 35 were used for all other elements, unless otherwise mentioned. For Ag, the optimized double-zeta basis set (so-called "p-valence" basis set) was used for Ag atoms. In this basis set, the default p-type polarization function is replaced with a bound unoccupied p-type orbital and its split-valence complement. The inclusion of 5p orbitals in the valence improves the chemistry and photochemistry as showed in a previous work. 36 To test the effects of basis set in ECD simulations, more complete basis sets were constructed by adding diffuse augmentation functions through truncated numerical Gaussian-type orbitals (NGTOs) to the default dzp basis sets. 40 We denote these basis sets dzp+NGTOs. Our approach follows a recent study of introducing augmentation functions that demonstrated good results for Bethe-Salpeter equation (BSE) and LR-TDDFT calculations for molecules with numeric atom-centered orbitals. 41 Gaussian basis function exponential parameters, the ζ-parameters, were taken from aug-cc-pvdz basis sets tabulated in Basis Set Exchange. 42,43 The parameters are tabulated in Supplementary Table S1.
For comparison, we also calculated the ECD by existing LR-TDDFT methods in GPAW. The LR-TDDFT approach of GPAW chooses a cut-off for the Kohn-Sham single-particle excitations and diagonalizes the Casida matrix, hence there exists a cut-off parameter in these calculations. We choose high cut-off (> 20 eV) to compare with our RT-TDDFT results in this work. The effect of the cut-off to the convergence of LR-TDDFT is discussed in Supplementary Note S1.
An artificial life-time for the electron dynamics was introduced via Gaussian line shape with σ = 0.2 eV in all figures unless otherwise mentioned. In this work, the rotatory strength is presented in units 10 −40 erg · esu · cm · Gauss −1 eV −1 = 10 −40 cgs eV −1 . The reported computational run times are obtained with Intel Xeon Gold 6230 processors with Mellanox HDR InfiniBand interconnect as installed in Puhti supercomputer at CSC -Finnish IT Center for Science. To support open data-driven chemistry and materials science, 44 we will upload all calculations of this work to the Novel Materials Discovery (NOMAD) laboratory and open-access Zenodo repository.

III. RESULTS
In this section, we present four test cases for our implementation. First, we use a benchmark molecule ((R)methyloxirane) to validate that our RT-TDDFT implementation can produce the same ECD spectra as the LR-TDDFT implementation in both LCAO and real-space grid mode. Then, we use a chiral Ag 4 string and a Ag +mediated guanine duplex G 2 −Ag 2+ 2 −G 2 structure 16,45 to demonstrate that LCAO RT-TDDFT adequately reproduces the rotatory strength of the reference grid mode calculation up to 8 eV. Finally, we apply the LCAO RT-TDDFT approach to a hybrid silver cluster [Ag 78 (S-BDPP) 6 (SR) 42 ] (hereafter denoted as Ag 78 ), where BDPP = 2,4-bis-(diphenylphosphino)pentane and SR = SPhCF 3 . 38 The whole system is considerably large for TDDFT simulations. We will show that LCAO RT-TDDFT is computationally efficient and produces ECD spectra that compare well with experimental results.
A. (R)-methyloxirane (R)-methyloxirane is one of the most typical benchmarks for optical activity calculations. 30,31 Therefore, we choose this chiral molecule as our first test system. The atomic structure is taken from the NIST database. The structure was optimized with the configuration interaction singles-doubles (CISD) method and a 6-31G* Gaussian orbital basis set. 46 The ECD of (R)-methyloxirane was calculated both with our RT-TDDFT implementation and LR-TDDFT. All RT-TDDFT calculations were propagated to T = 30 fs in steps of 5 as. The dzp+NGTO basis set was used in the LCAO simulations. The RT-TDDFT spectra look identical to the LR-TDDFT ones in both LCAO and real-space grid mode (the ECD is shown separately for LCAO and real-space grid cases in Supplementary Figure S1). The maximum difference is less than 0.5 in cgs units. The dzp+NGTO accurately predicts the four first peaks in comparison to real-space grid calculation, as shown in Figure 1, but the dzp basis doesn't give accurate results in this case (Supplementary Figure S2).

B. Ag4 string
To test our method on metallic systems, we use a chiral silver string as the second example. The Ag 4 string was built artificially, with a bond length of 2.7 Å, Ag-Ag-Ag angle of 150 • , and torsion angle of 10 • as shown in the inset of Figure 2. In the RT-TDDFT simulations, the propagation was carried out up to T = 30 fs in steps of 5 as. Figure 2a and Figure 2b show that our RT-TDDFT implementation again successfully reproduces the rotatory strength of LR-TDDFT both in LCAO and grid mode. The slight disagreement of LR-TDDFT and RT-TDDFT in grid mode above 4 eV is due to a difficult LR-TDDFT convergence, which is discussed more in detail in Supplementary Note S1. Figure 2c shows that LCAO again adequately reproduces the rotatory strength from the more accurate grid mode up to 8 eV, which is higher than energies commonly used for recording experimental spectra.

C. Ag + -mediated guanine duplex
After testing on a molecule and a silver string, we apply our method to an organic-metal hybrid system. We use one configuration of the Ag + -mediated guanine duplex (G 2 −Ag 2+ 2 −G 2 , Figure 3a) from our previous work. 16 The purpose here is to benchmark the accuracy of the LCAO method in a more complex system.
Comparing the results calculated with the two modes in Figure 3b, we find that the dzp basis set reproduces almost the same rotatory strength up to 6 eV, covering the energy window of most experimentally measured ECD spectra. The dzp+NGTO basis improves the agreement up to 8 eV as shown in Supplementary Figure S3.
The benefit of the LCAO mode is its low computational cost. For this system, the LCAO mode is over 10 times faster (9 hours on 80 cores versus 40 hours on 240 cores). Furthermore, we calculated the ECD with a coarser grid h = 0.3 Å to represent real-space density and potential and a larger time step 10 as. The ECD lies on top of the one obtained from previous RT/dzp as shown in Figure 3b. The simulation with coarser parameters took only 2.5 hours using 80 cores, which is about 50 times faster than the grid mode.

D. Ag78 cluster
In this test case, we illustrate the efficiency and accuracy of our RT-TDDFT/LCAO methodology on a ligandprotected Ag 78 cluster (Figure 4a) and present a comparison between the experimentally measured 38 and the calculated ECD spectra. We used the X-ray structure 38 in calculations. The 78 silver atoms (Figure 4b) are covered by ligand molecules containing C, N, O, F, H and S atoms (Figure 4a). The total number of atoms is 1074 and the number of valence electrons is 4272. The large size and the complexity of the cluster make it an ideal system to test the computational efficiency of the RT-TDDFT/LCAO approach.
Due to its O(N 5 ) scaling LCAO LR-TDDFT was not applicable to the Ag 78 cluster with our computational resources. However, the scaling of LCAO RT-TDDFT is only O(N 3 ), which made it possible to calculate ECD spectra for Ag 78 .
In addition to the PBE exchange-correlation functional, we also used the GLLB-SC exchange-correlation potential. 47 The GLLB-SC functional was chosen because former studies show that the GLLB-SC functional provides more accurate predictions of the optical absorption spectra of Ag clusters with respect to both the local density approximation (LDA) and the generalized gradient approximations (GGA). 36 The real-time propagation was taken up to T = 30 fs in 10 as steps and a grid   5. (a) Photoabsorption spectrum of Ag78. Shaded areas represent calculated spectra and lines shifted calculated spectra. The GLLB-SC spectrum was shifted by 0.14 eV and the PBE spectrum by 0.28 eV. Gaussian broadening with σ = 0.1 eV was applied. (b) The experimental (top panel) and calculated (lower panel) ECD spectra of Ag78. Default dzp basis sets were used for other than Ag atoms. Gaussian broadening with σ = 0.2 eV was applied to approximately match the spectral linewidth with the experimental data. The calculated ECD spectra were shifted according to the shifts done for calculated absorption spectra.
spacing of h = 0.3 Å was used.
For the Ag 78 cluster, we have also calculated the photoabsorption spectrum with LCAO RT-TDDFT (Figure 5a). Both GLLB-SC and PBE reproduce the first peak of the measured absorption spectrum. However, the GLLB-SC spectrum is red-shifted by 0.14 eV and PBE by 0.28 eV. Table I and Figure 5b present the comparison between the experimentally measured 38 and the calculated ECD spectra. The TDDFT spectra are shifted to higher energies by the same amount as the optical spectra in Figure  5a). Both the PBE and GLLB-SC functional capture the main features of the experimental ECD, which are the four positive peaks (a, c, e, f) and the two negative peaks (b, d).
We now briefly discuss the differences between the theoretical spectra and the experimental spectrum. The calculated absorption and ECD spectra are shifted to lower energies most likely due to the underestimation of the energy gap between occupied and unoccupied KS states in the DFT simulations and mismatches in the Ag dband location. The underestimation is less pronounced in the GLLB-SC calculations, because GLLB-SC introduces an orbital-energy dependent localization of the exchange hole and describes Ag d-orbitals more accurately. However, the improved description of the energy gap in GLLB-SC does not remove the mismatch of peaks d and e, suggesting there may be transitions that need a better description. Furthermore, the ECD spectrum was measured in a solvent. The fact that our calculations are performed for the experimental crystal X-ray structure and without conformational sampling may contribute to the differences. Using the dzp+NGTO basis set does not remove these differences as demonstrated in Supplementary Figure S4. The calculation of Ag 78 system took 24 hours with 200 cores with the PBE exchange-correlation functional and 33 hours with the GLLB-SC functional. This is remarkably fast for TDDFT ECD calculations of such a large system.

IV. CONCLUSIONS
We present a RT-TDDFT implementation for calculating ECD in GPAW package, which supports both LCAO mode and grid mode. While RT-TDDFT/LCAO is less accurate than RT-TDDFT/GRID, our tests have shown that the LCAO method nevertheless produces matching spectra in the experimentally relevant energy ranges.
The high computational efficiency of the RT-TDDFT/LCAO is enabled by the combination of localized orbitals and the PAW method. We demonstrated the efficiency of our code by computing ECD spectra of a large hybrid nanocluster with thousand atoms, a system whose ECD is challenging to compute by RT-TDDFT/GRID or conventional linear-response formalisms.
Our RT-TDDFT implementation with localized orbitals and PAW in GPAW opens the door to study the large-scale chiral systems with good accuracy and efficiency. We expect that our open-source implementation will be advantageous for studying the chiroptical property of large systems without excessive computational cost, which will help to develop many chirality related applications.