Real-Time Time-Dependent Density Functional Theory Implementation of Electronic Circular Dichroism Applied to Nanoscale Metal-Organic Clusters
Esko Makkonen, Tuomas P. Rossi, Ask Hjorth Larsen, Olga Lopez-Acevedo, Patrick Rinke, Mikael Kuisma, Xi Chen
RReal-Time Time-Dependent Density Functional Theory Implementation of ElectronicCircular Dichroism Applied to Nanoscale Metal-Organic Clusters
Esko Makkonen, Tuomas P. Rossi, Ask Hjorth Larsen, OlgaLopez-Acevedo, Patrick Rinke, Mikael Kuisma, ∗ and Xi Chen † Department of Applied Physics, Aalto University, Espoo, Finland Department of Physics, Chalmers University of Technology, Gothenburg, Sweden Simune Atomistics S.L., Donostia/San Sebastián, Spain Instituto de Física, Facultad de Ciencias Exactas y Naturales, Universidad de Antioquia, Medellín, Colombia Nanoscience Center, University of Jyväskylä, Jyväskylä, Finland (Dated: July 20, 2020)Electronic circular dichroism (ECD) is a powerful spectroscopical method for investigating chi-ral properties at the molecular level. ECD calculations with the commonly used linear-responsetime-dependent density functional theory (LR-TDDFT) framework can be prohibitively costly forlarge systems. To alleviate this problem, we present here an ECD implementation for the projec-tor augmented-wave method in the real-time-propagation TDDFT (RT-TDDFT) framework in theopen-source GPAW code. Our implementation supports both local atomic basis set and real-spacefinite-difference representations of wave functions. We benchmark our implementation against anexisting LR-TDDFT implementation in GPAW for small chiral molecules. We then demonstratethe efficiency of our local atomic basis set implementation for a large hybrid nanocluster.
I. INTRODUCTION
Chirality is an essential property in several branchesof science and technology. Chiral molecules play afundamental role in biological activities; for example,DNA double-helices are right-handed and amino acidsare left-handed. Chirality is also critical in pharmaceuti-cals. For example, R -enantiomer thalidomide is effectiveagainst morning sickness for pregnant women, but the S -species produce fetal deformations. In addition, chi-ral molecules and chiral nanomaterials have many poten-tial applications in catalysis, sensors, spintronics, opto-electronics, and nanoelectronics.
Therefore, the de-termination of the handedness of chiral systems is ofparamount importance.Chiral molecules absorb left and right circular polar-izations of light differently. This difference is probed inelectronic circular dichroism (ECD) spectroscopy. ECDis defined as ∆ (cid:15) = (cid:15) L − (cid:15) R , where (cid:15) L and (cid:15) R are the molarextinction coefficients for left and right circularly polar-ized light, respectively. Since ECD is highly sensitive tosmall details in the atomic structure of molecules andunique for each conformation, it is a powerful techniquefor characterizing chiral systems and for distinguishingenantiomers. The ECD spectroscopy accompanied by computa-tional modeling can help to get insightful knowledgeof the atomic structures of chiral biomolecules andnanoclusters.
For example, in our earlier work, wehave identified by comparing simulated and measuredECD that the Ag + -mediated guanine duplex has the left-handed helix configuration, while the Ag + -mediatedcytosine has the right-handed helix configuration. Most computational ECD approaches are based on thelinear-response formalism. Time-dependent density-functional theory (TDDFT) has become the linear- response method of choice due to its favorable balance ofaccuracy and computational cost, compared to quantumchemical approaches such as coupled cluster and con-figuration interaction methods. In linear-responseTDDFT (LR-TDDFT), the Casida equation issolved in the basis of Kohn–Sham (KS) particle–holetransitions in the frequency domain.
A full ECDspectrum requires the calculation of a large number oftransitions from occupied to unoccupied states.This be-comes computationally prohibitive for large systems witha high density of states, resulting in an unfavorable O ( N ) scaling, where N is the system size.An alternative to LR-TDDFT is real-time-propagationtime-dependent density-functional theory (RT-TDDFT).In RT-TDDFT, the system is subjected to an initial per-turbation and the KS wave functions are propagated inthe time domain by numerically integrating the time-dependent KS equations. The real-time approach cap-tures the same information as LR-TDDFT, for small ini-tial perturbations, and incorporates nonlinear spectralinformation for larger initial perturbations. RT-TDDFT scales as O ( N ) , but suffers from a largeprefactor. LR-TDDFT is usually faster for small systemssuch as small organic molecules. For large molecules,clusters, and nanoparticles, RT-TDDFT becomes morecost-effective than LR-TDDFT. RT-TDDFT computation of ECD has been imple-mented for a variety of basis sets: real-space grids,
Gaussian-type atomic orbitals, and a mix of Gaussian-type and plane-wave basis sets. In this work, we presenta RT-TDDFT ECD implementation in the open-sourceGPAW package.
Our implementation uses the pro-jector augmented wave (PAW) method and supportsboth localized basis sets (LCAO mode) and real-spacegrids (grid mode). We verify our implementation bycomparing our results to those calculated with the ex-isting LR-TDDFT implementation in GPAW. We also a r X i v : . [ phy s i c s . c h e m - ph ] J u l benchmark the LCAO mode against accurate real-spacegrid calculations.In GPAW, the time-dependent density and potentialare expressed on a uniform grid, and the matrix elementsof the potential are evaluated on this grid. The smooth-ness of these quantities allows for a coarse grid spac-ing. The LCAO-PAW pseudo wavefunctions can forma local and efficient representation suitable for systemswith hundreds of atoms. Previous work has shown thatLCAO RT-TDDFT in GPAW is capable of simulatingthe optical spectrum of a silver cluster of more than 500atoms (Ag ). In this work, we demonstrate the ef-ficiency of our LCAO RT-TDDFT ECD implementationfor a large ligand-protected Ag cluster consisting ofover 1000 atoms and over 4000 electrons.The rest of this paper is organized as follows. In theMethods section we illustrate our implemented method-ologies by showing how the ECD is calculated from time-dependent magnetic dipole moment and how the time-dependent magnetic dipole moment is calculated withRT-TDDFT in GPAW. The important information con-sidering all simulations is described in the ComputationalMethods section. In the Results section we demonstratethe capability of our implementations to predict ECDspectrum for four test cases. Finally, we give a summaryin the Conclusion section. II. METHODS
In RT-TDDFT, the KS wave functions are propagatedin time in response to a time-dependent potential startingfrom an initial state, here chosen to be the ground state.The time-dependent KS equation is defined as i ∂∂t ψ n ( r , t ) = H KS ( t )Ψ n ( r , t ) , (1)where H KS ( t ) is the time-dependent KS Hamiltonian and ψ n ( r , t ) is a time-dependent KS wave function. A com-mon practice in time-propagation schemes is to use theweak δ -kick approach to calculate the linear-responsefunctions. 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 beendemonstrated previously. In this work, we implement and thoroughly bench-mark the calculation of time-dependent magnetic mo-ment within the time-propagation framework for obtain-ing the ECD spectrum. In the following, we derivethe relevant equations within the PAW method. Thisderivation partly follows the one shown by Varsano etal. A. ECD from the Induced Time-DependentMagnetic Moment
A commonly used experimental quantity to measureECD is the difference in molar extinction coefficients,which is given by ∆ (cid:15) ( ω ) = 16 πN A π (cid:126) c ωR ( ω ) cgs . (2)Here ω is the the energy of the incident light, c the speedof light, (cid:126) the reduced Planck constant, N A Avogadro’sconstant, and R ( ω ) cgs the rotatory strength in cgs units.The quantity that characterizes ∆ (cid:15) ( ω ) and therefore theECD spectrum is the rotatory strength. The relation-ship between rotatory strength in cgs units and rotatorystrength in atomic units (denoted as R ( ω ) ) is R ( ω ) cgs = e (cid:126) m e α R ( ω ) , (3)where e is the elementary charge, m e the mass of anelectron, and α the fine structure constant. We will workin atomic units and perform the required unit conversionsafterwards.The rotatory strength is defined through the opticalrotatory response tensor by R ( ω ) = ωπc Im (cid:20) (cid:88) k β kk ( ω ) (cid:21) , (4)where index k enumerates Cartesian coordinates ( k ∈{ x, y, z } ). Next we will derive the relationship be-tween the optical rotatory response tensor and the time-dependent magnetic dipole moment, which is the quan-tity calculated in TDDFT.For a system in an external electric field E =[ E x , E y , E z ] (no external magnetic field present), the in-duced time-dependent magnetic dipole moment in direc-tion j , m j ( t ) , has the expansion m j ( t ) = 1 c (cid:88) k (cid:90) ∞−∞ β jk ( t − τ ) ∂E k ( τ ) ∂τ d τ + higher-order terms (5)where β jk is the jk component (indices j, k ∈ { x, y, z } ) ofthe optical rotatory response tensor, and E k is the elec-tric field component. The response tensor β jk describesthe induced magnetic dipole moment in the j directionfor a perturbing electric field in the k direction. The re-sponse is causal, which means that β jk ( t ) vanishes fornegative values of t .In the weak-field limit, the time-dependent magneticmoment is dominated by the first-order term given bythe linear-response functions β jk ( t ) . By using the convo-lution theorem and the properties of Fourier transforms,the first-order term of Eq. (5) can be written in the fre-quency domain as m j ( ω ) = − iωc (cid:88) k β jk ( ω ) E k ( ω ) . (6)To resolve all components of β jk ( t ) , we perform the δ -kick in all three Cartesian directions using a perturbingelectric field of the form E ( k ) ( t ) = κ ˆ k δ ( t ) . The k super-script in parenthesis indicates the kick direction, to bedistinguished from the component subscript. We keepthe intensity κ weak to restrict our calculations to thelinear-response regime.In the frequency domain, the delta kick becomes a con-stant for all frequencies E ( k ) k (cid:48) ( ω ) = κδ kk (cid:48) . (7)Eq. (6) then simplifies to β jk ( ω ) = icκω m ( k ) j ( ω ) . (8)Combining Eqs. (4) and (8), we get the equation for ro-tatory strength expressed with magnetic dipole moment R ( ω ) = 1 πκ Re (cid:20) (cid:88) k m ( k ) k ( ω ) (cid:21) . (9)In our methodology, m ( k ) j ( ω ) is calculated by Fouriertransform of m ( k ) j ( t ) , which has been obtained throughtime-propagation, m ( k ) j ( ω ) = (cid:90) ∞ e iωt m ( k ) j ( t ) d t. (10)In principle, the integration interval goes from zero toinfinity. In practice, a finite propagation time ( T ) sufficesby introducing an artificial lifetime ω → ω + i σ t , where σ is the parameter that determines the line width of theGaussian line shape. Introducing this into Eq. (10) gives m ( k ) j ( ω ) = (cid:90) T e iωt e − σ t m ( k ) j ( t ) d t. (11)For a desired value of σ , the propagation time T needsto be large enough so that e − σ T ≈ . B. Computing the Magnetic Moment m ( k ) j ( t ) in thePAW formalism The magnetic moment is defined by the following op-erator (in atomic units) ˆ m = − i c ˆ r × ˆ ∇ . (12)The expectation value of the time-dependent magneticmoment is obtained as m ( t ) = (cid:88) n f n (cid:90) ψ ∗ n ( r , t ) ˆ m ψ n ( r , t ) d r , (13)where f n is the occupation number of the n :th KS stateand ψ n ( r , t ) is the time-evolved KS wave function. In the PAW method the wave functions ψ n ( r , t ) havedecomposition ψ n ( r , t ) = ˜ ψ n ( r , t )+ (cid:88) ai (cid:104) φ ai ( r ) − ˜ φ ai ( r ) (cid:105) (cid:104) ˜ p ai | ˜ ψ n ( t ) (cid:105) , (14)where ˜ ψ n ( r , t ) is a smooth pseudo wave function and (cid:80) ai (cid:104) φ ai ( r ) − ˜ φ ai ( r ) (cid:105) (cid:104) ˜ p ai | ˜ ψ n ( t ) (cid:105) a local correction insidean atomic augmentation sphere . ˜ p ai is a localized pro-jector function and φ ai and ˜ φ ai are partial and pseudopartial waves , respectively. These quantities are specificto PAW. In the PAW formalism, the expectation valuein Eq. (13) becomes m ( t ) = (cid:88) n f n (cid:104) ˜ ψ n ( t ) | ˆ m | ˜ ψ n ( t ) (cid:105) + (cid:88) naij f n (cid:104) ˜ ψ n ( t ) | ˜ p ai (cid:105) ∆ M aij (cid:104) ˜ p aj | ˜ ψ n ( t ) (cid:105) , (15)where the augmentation-sphere contribution is ∆ M aij = (cid:104) φ ai | ˆ m | φ aj (cid:105) − (cid:104) ˜ φ ai | ˆ m | ˜ φ aj (cid:105) . For evaluating ∆ M aij , the required matrix elements of the form (cid:104) φ ai | r × ∇| φ aj (cid:105) are evaluated in two atom-centered partsas (cid:104) φ ai | ( r − R a ) × ∇| φ aj (cid:105) + R a × (cid:104) φ ai |∇| φ aj (cid:105) , where R a isthe coordinate of atom a .In the LCAO expansion, the time-dependent pseudowave function ˜ ψ n ( r , t ) is written as a linear combinationof atom-centered basis functions ϕ µ ( r − R a )˜ ψ n ( r , t ) = (cid:88) µ c µn ( t ) ϕ µ ( r − R a ) , (16)where c µn ( t ) are the time-dependent expansion coeffi-cients. With this LCAO expansion, Eq. (15) can be writ-ten compactly as m ( t ) = (cid:88) µν ρ νµ ( t ) M µν , (17)where ρ µν ( t ) = (cid:80) n f n c µn ( t ) c ∗ νn ( t ) is the KS density ma-trix in the LCAO basis. The matrix elements M µν aregiven by the pseudo and augmentation contributions M µν = ˜ M µν + ∆ M µν (18) ˜ M µν = (cid:104) ϕ µ | ˆ m | ϕ ν (cid:105) (19) ∆ M µν = (cid:88) aij (cid:104) ϕ µ | ˜ p ai (cid:105) ∆ M aij (cid:104) ˜ p aj | ϕ ν (cid:105) . (20)In real-space grid mode the magnetic moment is calcu-lated using equation Eq. (15). In LCAO mode, Eq. (17)is used. The matrix elements M µν are time independentand calculated only once before the time propagation.After a complete time propagation, the recorded m ( t ) is transformed to frequency domain as a post processingstep according to Eq. (11) at each desired ω value. Thenthe rotatory strength is calculated according to Eq. (9).We note that research of gauge origin issues is notwithin the scope of this work and the reader is suggestedto explore a recent detailed investigation on the matter. C. Computational Methods
For our calculations in this work, we used the PBEexchange-correlation functional, unless otherwise men-tioned. The molecules were placed into a cubic unit cellwith the vacuum size of 8 Å. The real-space grid spac-ing was chosen as h = 0 . Å. We tested coarser settingswith h = 0 . Å in RT/LCAO mode for the Ag + -mediatedguanine duplex case to demonstrate that such a coarsergrid is sufficient for calculating the ECD spectrum withinLCAO mode, where the grid is used to represent onlyreal-space density and potential. Per atom, the electronic configuration of valenceelectrons is H(1s ) O(2s ), C(2s ), N(2s ),S(3s ), P(3s ), F(2s ) and Ag(4d ). Theremaining electrons were treated as a frozen core. Thedefault PAW dataset package 0.9.20000 was used for allthe atoms.In the LCAO mode, the default GPAW double-zetapolarized (dzp) basis sets were used for all other el-ements, unless otherwise mentioned. For Ag, the opti-mized double-zeta basis set (so-called "p-valence" basisset) was used for Ag atoms. In this basis set, the defaultp-type polarization function is replaced with a bound un-occupied p-type orbital and its split-valence complement.The inclusion of 5p orbitals in the valence improves thechemistry and photochemistry as showed in a previouswork. To test the effects of basis set in ECD simulations,more complete basis sets were constructed by adding dif-fuse augmentation functions through truncated numeri-cal Gaussian-type orbitals (NGTOs) to the default dzpbasis sets. We denote these basis sets dzp+NGTOs.Our approach follows a recent study of introducingaugmentation functions that demonstrated good resultsfor Bethe-Salpeter equation (BSE) and LR-TDDFTcalculations for molecules with numeric atom-centeredorbitals. Gaussian basis function exponential parame-ters, the ζ -parameters, were taken from aug-cc-pvdz basissets tabulated in Basis Set Exchange. The parame-ters are tabulated in Supplementary Table S1.For comparison, we also calculated the ECD by ex-isting LR-TDDFT methods in GPAW. The LR-TDDFTapproach of GPAW chooses a cut-off for the Kohn-Shamsingle-particle excitations and diagonalizes the Casidamatrix, hence there exists a cut-off parameter in thesecalculations. We choose high cut-off ( >
20 eV) to com-pare with our RT-TDDFT results in this work. The ef-fect of the cut-off to the convergence of LR-TDDFT isdiscussed in Supplementary Note S1.An artificial life-time for the electron dynamics was introduced via Gaussian line shape with σ = 0 . eV inall figures unless otherwise mentioned. In this work, therotatory strength is presented in units − erg · esu · cm · Gauss − eV − = 10 − cgs eV − .The reported computational run times are obtainedwith Intel Xeon Gold 6230 processors with MellanoxHDR InfiniBand interconnect as installed in Puhti su-percomputer at CSC – Finnish IT Center for Science.To support open data-driven chemistry and materialsscience, we will upload all calculations of this work tothe Novel Materials Discovery (NOMAD) laboratory andopen-access Zenodo repository. III. RESULTS
In this section, we present four test cases for our im-plementation. First, we use a benchmark molecule ( ( R ) -methyloxirane) to validate that our RT-TDDFT imple-mentation can produce the same ECD spectra as the LR-TDDFT implementation in both LCAO and real-spacegrid mode. Then, we use a chiral Ag string and a Ag + -mediated guanine duplex G − Ag − G structure to demonstrate that LCAO RT-TDDFT adequately re-produces the rotatory strength of the reference gridmode calculation up to 8 eV. Finally, we apply theLCAO RT-TDDFT approach to a hybrid silver clus-ter [ Ag ( S -BDPP) (SR) ] (hereafter denoted as Ag ),where BDPP = 2,4-bis-(diphenylphosphino)pentane andSR = SPhCF . The whole system is considerably largefor TDDFT simulations. We will show that LCAO RT-TDDFT is computationally efficient and produces ECDspectra that compare well with experimental results. A. ( R ) -methyloxirane ( R ) -methyloxirane is one of the most typical bench-marks for optical activity calculations. Therefore, wechoose this chiral molecule as our first test system. Theatomic structure is taken from the NIST database. The
FIG. 1. Rotatory strength of ( R ) -methyloxirane calculated byRT-TDDFT and LR-TDDFT in both LCAO and real-spacegrid modes. FIG. 2. Rotatory strength of a Ag -string (inset figures) calculated with (a) the LCAO mode and (b) the grid mode. (c) Com-parison between the two modes. structure was optimized with the configuration interac-tion singles-doubles (CISD) method and a 6-31G* Gaus-sian orbital basis set. The ECD of ( R ) -methyloxirane was calculated bothwith 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 wasused in the LCAO simulations. The RT-TDDFT spec-tra look identical to the LR-TDDFT ones in both LCAOand real-space grid mode (the ECD is shown separatelyfor LCAO and real-space grid cases in SupplementaryFigure S1). The maximum difference is less than 0.5 incgs units. The dzp+NGTO accurately predicts the fourfirst peaks in comparison to real-space grid calculation,as shown in Figure 1, but the dzp basis doesn’t give ac-curate results in this case (Supplementary Figure S2). B. Ag string To test our method on metallic systems, we use a chiralsilver string as the second example. The Ag string wasbuilt artificially, with a bond length of 2.7 Å, Ag-Ag-Agangle of 150 ◦ , and torsion angle of 10 ◦ as shown in theinset of Figure 2. In the RT-TDDFT simulations, thepropagation was carried out up to T = 30 fs in steps of5 as.Figure 2a and Figure 2b show that our RT-TDDFT im-plementation again successfully reproduces the rotatorystrength of LR-TDDFT both in LCAO and grid mode.The slight disagreement of LR-TDDFT and RT-TDDFTin grid mode above 4 eV is due to a difficult LR-TDDFTconvergence, which is discussed more in detail in Supple-mentary Note S1.Figure 2c shows that LCAO again adequately repro-duces the rotatory strength from the more accurate gridmode up to 8 eV, which is higher than energies commonlyused for recording experimental spectra. C. Ag + -mediated guanine duplex After testing on a molecule and a silver string, we applyour method to an organic-metal hybrid system. We useone configuration of the Ag + -mediated guanine duplex(G − Ag − G , Figure 3a) from our previous work. The purpose here is to benchmark the accuracy of theLCAO method in a more complex system.Comparing the results calculated with the two modesin Figure 3b, we find that the dzp basis set reproduces al-most the same rotatory strength up to 6 eV, covering theenergy window of most experimentally measured ECDspectra. The dzp+NGTO basis improves the agreementup to 8 eV as shown in Supplementary Figure S3.The benefit of the LCAO mode is its low computa-tional cost. For this system, the LCAO mode is over10 times faster (9 hours on 80 cores versus 40 hours on240 cores). Furthermore, we calculated the ECD witha coarser grid h = 0 . Å to represent real-space densityand potential and a larger time step 10 as. The ECDlies on top of the one obtained from previous RT/dzp asshown in Figure 3b. The simulation with coarser param-eters took only 2.5 hours using 80 cores, which is about50 times faster than the grid mode.
D. Ag cluster In this test case, we illustrate the efficiency and accu-racy of our RT-TDDFT/LCAO methodology on a ligand-protected Ag cluster (Figure 4a) and present a com-parison between the experimentally measured and thecalculated ECD spectra. We used the X-ray structure in calculations. The 78 silver atoms (Figure 4b) are cov-ered by ligand molecules containing C, N, O, F, H andS atoms (Figure 4a). The total number of atoms is 1074and the number of valence electrons is 4272. The largesize and the complexity of the cluster make it an idealsystem to test the computational efficiency of the RT-TDDFT/LCAO approach.Due to its O ( N ) scaling LCAO LR-TDDFT was notapplicable to the Ag cluster with our computational FIG. 3. (a) Structure and (b) rotatory strength ofG − Ag − G . LCAO calculations use dzp basis sets. Here † notes settings with grid parameter h = 0 . Å, and propaga-tion of 30 fs in steps of 10 as. resources. However, the scaling of LCAO RT-TDDFTis only O ( N ) , which made it possible to calculate ECDspectra for Ag .In addition to the PBE exchange-correlation func-tional, we also used the GLLB-SC exchange-correlationpotential. The GLLB-SC functional was chosen be-cause former studies show that the GLLB-SC functionalprovides more accurate predictions of the optical absorp-tion spectra of Ag clusters with respect to both the localdensity approximation (LDA) and the generalized gradi-ent approximations (GGA). The real-time propagationwas taken up to T = 30 fs in 10 as steps and a grid FIG. 4. (a) The structure of the Ag cluster. Ligand atoms:H, white; C, beige; F, green; P, orange; S, yellow. (b) Agatoms in the Ag cluster. Energy (eV) A b s o r ban c e Exp.GLLB-SC shiftedPBE shiftedGLLB-SCPBE
FIG. 5. (a) Photoabsorption spectrum of Ag . Shaded areasrepresent calculated spectra and lines shifted calculated spec-tra. The GLLB-SC spectrum was shifted by 0.14 eV and thePBE spectrum by 0.28 eV. Gaussian broadening with σ = 0 . eV was applied. (b) The experimental (top panel) and calcu-lated (lower panel) ECD spectra of Ag . Default dzp basissets were used for other than Ag atoms. Gaussian broaden-ing with σ = 0 . eV was applied to approximately match thespectral linewidth with the experimental data. The calcu-lated ECD spectra were shifted according to the shifts donefor calculated absorption spectra. spacing of h = 0 . Å was used.For the Ag cluster, we have also calculated the pho-toabsorption spectrum with LCAO RT-TDDFT (Fig-ure 5a). Both GLLB-SC and PBE reproduce the firstpeak of the measured absorption spectrum. However, theGLLB-SC spectrum is red-shifted by 0.14 eV and PBEby 0.28 eV.Table I and Figure 5b present the comparison betweenthe experimentally measured and the calculated ECDspectra. The TDDFT spectra are shifted to higher ener-gies by the same amount as the optical spectra in Figure5a). Both the PBE and GLLB-SC functional capturethe main features of the experimental ECD, which arethe four positive peaks (a, c, e, f) and the two negativepeaks (b, d).We now briefly discuss the differences between the the-oretical spectra and the experimental spectrum. The cal-culated absorption and ECD spectra are shifted to lowerenergies most likely due to the underestimation of theenergy gap between occupied and unoccupied KS statesin the DFT simulations and mismatches in the Ag d-band location. The underestimation is less pronouncedin the GLLB-SC calculations, because GLLB-SC intro-duces an orbital-energy dependent localization of the ex-change hole and describes Ag d-orbitals more accurately. TABLE I. The ECD peak positions in Figure 5 (in unit eV).
Peak Experiment PBE GLLB-SCa b c d e f However, the improved description of the energy gap inGLLB-SC does not remove the mismatch of peaks d ande, suggesting there may be transitions that need a betterdescription. Furthermore, the ECD spectrum was mea-sured in a solvent. The fact that our calculations areperformed for the experimental crystal X-ray structureand without conformational sampling may contribute tothe differences. Using the dzp+NGTO basis set does notremove these differences as demonstrated in Supplemen-tary Figure S4.The calculation of Ag system took 24 hours with 200cores with the PBE exchange-correlation functional and33 hours with the GLLB-SC functional. This is remark-ably fast for TDDFT ECD calculations of such a largesystem. IV. CONCLUSIONS
We present a RT-TDDFT implementation for calculat-ing ECD in GPAW package, which supports both LCAO mode and grid mode. While RT-TDDFT/LCAO is lessaccurate than RT-TDDFT/GRID, our tests have shownthat the LCAO method nevertheless produces matchingspectra in the experimentally relevant energy ranges.The high computational efficiency of the RT-TDDFT/LCAO is enabled by the combination of lo-calized orbitals and the PAW method. We demon-strated the efficiency of our code by computing ECDspectra of a large hybrid nanocluster with thousandatoms, a system whose ECD is challenging to compute byRT-TDDFT/GRID or conventional linear-response for-malisms.Our RT-TDDFT implementation with localized or-bitals and PAW in GPAW opens the door to study thelarge-scale chiral systems with good accuracy and effi-ciency. We expect that our open-source implementationwill be advantageous for studying the chiroptical prop-erty of large systems without excessive computationalcost, which will help to develop many chirality relatedapplications.
ACKNOWLEDGMENTS
This work was supported by the Academy of Finland,Projects 308647, 314298, 279240 and 312556. T.P.R.acknowledges support from the European Union’s Hori-zon 2020 research and innovation programme under theMarie Skłodowska-Curie grant agreement No 838996.M.K. acknowledges funding from Academy of Finlandunder grant No 295602. We acknowledge computationalresources provided by the CSC – IT Center for Science(Finland), the Aalto Science-IT project, and the SwedishNational Infrastructure for Computing (SNIC) at PDC(Stockholm). ∗ mikael.kuisma@jyu.fi † xi.6.chen@aalto.fi L. D. Barron, Space Sci. Rev. , 187 (2008). E. Tokunaga, T. Yamamoto, E. Ito, and N. Shibata, Sci.Rep. , 17131 (2018). A. L. Spek, Acta Crystallogr. B , 659 (2016). W. A. Nugent, T. V. RajanBabu, and M. J. Burk, Science , 479 (1993). T. J. Davis and D. E. Gómez, Phys. Rev. B , 235424(2014). C. Gautier and T. Bürgi, ChemPhysChem , 483 (2009). J. R. Sanchez-Valencia, T. Dienel, O. Gröning, I. Sho-rubalko, A. Mueller, M. Jansen, K. Amsharov, P. Ruffieux,and R. Fasel, Nature , 61 (2014). M. Bode, M. Heide, K. von Bergmann, P. Ferriani,S. Heinze, G. Bihlmayer, A. Kubetzka, O. Pietzsch,S. Blügel, and R. Wiesendanger, Nature , 190 (2007). C. Noguez and I. L. Garzón, Chem. Soc. Rev. , 757(2009). S. Knoppe and T. Bürgi, Acc. Chem. Res. , 1318 (2014). D. Zerrouki, J. Baudry, D. Pine, P. Chaikin, and J. Bi-bette, Nature , 380 (2008). H.-E. Lee, H.-Y. Ahn, J. Mun, Y. Y. Lee, M. Kim, N. H. Cho, K. Chang, W. S. Kim, J. Rho, and K. T. Nam,Nature , 360 (2018). J. T. A. Jones, T. Hasell, X. Wu, J. Bacsa, K. E. Jelfs,M. Schmidtmann, S. Y. Chong, D. J. Adams, A. Trewin,F. Schiffman, F. Cora, B. Slater, A. Steiner, G. M. Day,and A. I. Cooper, Nature , 367 (2011). J. Kypr, I. Kejnovská, D. Ren˘ciuk, and M. Vorlí˘cková,Nucleic Acids Res. , 1713 (2009). L. Chang, O. Baseggio, L. Sementa, D. Cheng, G. Fron-zoni, D. Toffoli, E. Aprà, M. Stener, and A. Fortunelli, J.Chem. Theory Comput. , 3703 (2018). X. Chen, E. Makkonen, D. Golze, and O. Lopez-Acevedo,J. Phys. Chem. Lett. , 4789 (2018). L. A. Espinosa Leal, A. Karpenko, S. Swasey, E. G. Gwinn,V. Rojas-Cervellera, C. Rovira, and O. Lopez-Acevedo, J.Phys. Chem. Lett. , 4061 (2015). T. D. Crawford, Theor. Chem. Acc. , 227 (2006). E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997(1984). C. Diedrich and S. Grimme, J. Phys. Chem. A , 2524(2003). J. Autschbach, T. Ziegler, S. J. A. van Gisbergen, andE. J. Baerends, J. Chem. Phys. , 6930 (2002). M. E. Casida, in
Recent Advances in Density FunctionalMethods, Part I , edited by D. P. Chong (World Scientific,Singapore, 1995) p. 155. M. R. Provorse and C. M. Isborn, Int. J. Quantum Chem. , 739 (2016). M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Phys.Rev. Lett. , 1212 (1996). M. E. Casida, J. Mol. Struct.: THEOCHEM , 3(2009). K. Yabana and G. F. Bertsch, Phys. Rev. B , 4484(1996). K. Yabana and G. F. Bertsch, Phys. Rev. A , 1271(1999). J. J. Goings, P. J. Lestrange, and X. Li, WIREs Comput.Mol Sci. , e1341 (2018). S. Tussupbayev, N. Govind, K. Lopata, and C. J. Cramer,J. Chem. Theory Comput. , 1102 (2015). D. Varsano, L. A. Espinosa-Leal, X. Andrade, M. A. L.Marques, R. di Felice, and A. Rubio, Phys. Chem. Chem.Phys. , 4481 (2009). J. Mattiat and S. Luber, Chem. Phys. , 110464 (2019). J. Enkovaara, C. Rostgaard, J. Mortensen, J. Chen,M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad,V. Haikola, H. Hansen, and et al , J. Phys. Condens. Mat-ter , 253202 (2010). M. Walter, H. Häkkinen, L. Lehtovaara, M. Puska,J. Enkovaara, C. Rostgaard, and J. Mortensen, J. Chem.Phys. , 244101 (2008). P. E. Blöchl, Phys. Rev. B , 17953 (1994). A. H. Larsen, M. Vanin, J. J. Mortensen, K. S. Thygesen,and K. W. Jacobsen, Phys. Rev. B , 195112 (2009). M. Kuisma, A. Sakko, T. P. Rossi, A. H. Larsen, J. Enko-vaara, L. Lehtovaara, and T. T. Rantala, Phys. Rev. B , 115431 (2015). T. P. Rossi, M. Kuisma, M. J. Puska, R. M. Nieminen,and P. Erhart, J. Chem. Theory Comput. , 4779 (2017). H. Yang, J. Yan, Y. Wang, G. Deng, H. Su, X. Zhao, C. Xu,B. K. Teo, and N. Zheng, J. Am. Chem. Soc. , 16113(2017). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). T. P. Rossi, S. Lehtola, A. Sakko, M. J. Puska, and R. M.Nieminen, J. Chem. Phys. , 094114 (2015). C. Liu, J. Kloppenburg, Y. Yao, X. Ren, H. Appel,Y. Kanai, and V. Blum, J. Chem. Phys. , 044105(2020). K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun,V. Gurumoorthi, J. Chase, J. Li, and T. L. Windus, J.Chem. Inf. Model. , 1045 (2007). B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibson, andT. L. Windus, J. Chem. Inf. Model. , 4814 (2019). L. Himanen, A. Geurts, A. S. Foster, and P. Rinke, Adv.Sci. , 1900808 (2019). E. Makkonen, P. Rinke, O. Lopez-Acevedo, and X. Chen,Int. J. Mol. Sci. , 2346 (2018). “National Institute of Standards and Technology, Compu-tational Chemistry Comparison and Benchmark DataBase,Release 20. CISD/6-31G* calculated geometry for C H O(Propylene oxide) A C ,” . M. Kuisma, J. Ojanen, J. Enkovaara, and T. Rantala,Phys. Rev. B82