Method for Calculating Excited Electronic States Using Density Functionals and Direct Orbital Optimization with Real Space Grid or Plane Wave Basis Set
Aleksei V. Ivanov, Gianluca Levi, Elvar ?. Jónsson, Hannes Jónsson
DDirect Optimization Method for VariationalExcited-State Density Functional CalculationsUsing Real Space Grid or Plane Waves
Aleksei V. Ivanov, Gianluca Levi, Elvar Ö. Jónsson, and Hannes Jónsson ∗ Science Institute and Faculty of Physical Sciences, University of Iceland, VR-III, 107Reykjavík, Iceland
E-mail: [email protected]
Abstract
A direct optimization method is presented for density functional calculations ofexcited electronic states using either a real space grid or a plane wave basis set. Themethod is variational, provides atomic forces in the excited states, and can be applied toKohn-Sham (KS) functionals as well as orbital-density dependent functionals (ODD)including explicit self-interaction correction. The implementation for KS functionalsinvolves two nested loops: (1) An inner loop for finding a stationary point in a sub-space spanned by the occupied and a few virtual orbitals corresponding to the excitedstate; (2) an outer loop for minimizing the energy in a tangential direction. For ODDfunctionals, a third loop is used to find the unitary transformation that minimizes theenergy functional among occupied orbitals only. Combined with the maximum overlapmethod, the algorithm converges in challenging cases where conventional self-consistentfield algorithms tend to fail. The benchmark tests presented include two charge-transferexcitations in nitrobenzene and an excitation of CO to degenerate π ∗ orbitals where theimportance of complex orbitals is illustrated. An application of the method to several a r X i v : . [ phy s i c s . c o m p - ph ] F e b etal-to-ligand charge-transfer and metal-centred excited states of an Fe II photosensi-tizer complex is described and the results compared to reported experimental estimates.The method is also used to study the effect of Perdew-Zunger self-interaction correc-tion on valence and Rydberg excited states of several molecules, both singlet and tripletstates. The correction is found to improve the description of molecular bond stretchingbut calculated values of the excitation energy are improved only slightly, by ca. Density functional theory (DFT) is commonly used in computational studies of molecules andmaterials as it can in many cases give reasonable accuracy without too much computationaleffort. Calculations of ground electronic states can be performed even by non-experts thanksto well-established algorithms and software implementations. This does not, however, applyto calculations of excited electronic states, although such states are of great importance inmany rapidly developing fields such as ultrafast spectroscopy, solar energy conversion andphotocatalysis. Even though DFT is formulated theoretically as a ground state theory, itcan also be used to provide useful estimates of excited states. The most commonly usedexcited-state extension of DFT is time-dependent density functional theory (TDDFT).
In practical implementations, TDDFT calculations are carried out using some ground-statedensity functional, linear-response theory and an adiabatic approximation that neglects thetime dependence of the exchange-correlation (XC) kernel. Within these approximations,TDDFT typically provides a fairly good description of low-lying valence excitations, butoften fails to describe higher excitations and conical intersections between the ground andexcited states. Alternative approaches with wide applicability and similar computational effort can bebased on time-independent DFT such as ensemble DFT, excited-state DFT (eDFT) ∆ self-consistent field, ∆ SCF), the quasi-particle energy DFT(QE-DFT) approach , as well as constrained DFT, orthogonality constrained, andconstricted DFT approaches. There, an excited state is found from a stationary solutionof self-consistent equations and the computational robustness strongly depends on the algo-rithm used. In eDFT, the excited state is a solution of the Kohn-Sham (KS) equations fornon-aufbau orbital occupation numbers. Commonly used methods in ground-state calcula-tions are based on some iterative eigensolver such as the Davidson algorithm or the residualminimization method–direct inversion in the iterative subspace (RMM-DIIS) enhancedin various ways to improve robustness and rate of convergence. However, these algo-rithms are not specifically designed for calculations of excited states. The maximum overlapmethod (MOM) can be used to reduce the probability of convergence on the ground statein the iterative calculation, but convergence problems often occur. The basic problem liesin the fact that excited-state calculations do not involve finding the global minimum of theenergy as a function of the electronic degrees of freedom, but rather a more general station-ary point on the high-dimensional electronic energy surface. A method that can in generalconverge on a saddle point on the energy surface rather than a minimum is required. To findan N-th order saddle point, one needs to maximize the energy with respect to N degrees offreedom while minimizing with respect to all the others. The degrees of freedom along whichthe energy needs to be maximized are not known a priori . This makes a search for a saddlepoint significantly more difficult than a search for a minimum. Therefore, while calculationsbased on SCF-MOM-type algorithms can in principle converge on excited states, they are inpractice not reliable.Alternatively, direct optimization (DO) can be used to converge on solutions of the KSequations. While this approach has mainly been used in energy minimization to find aground-state solution , it can be extended to calculations of excited states.
Onepossible approach is to formulate the problem as a search for a minimum of the norm ofthe gradient. But, it is then important to also ensure that the norm of the gradient is3ero at the minimum and the evaluation of the gradient of this objective function increasesthe computational effort significantly. A different approach is based on direct optimizationof the energy using a quasi-Newton method that can develop negative eigenvalues of theHessian consistent with the type of saddle point searched for.
A preconditioner is thenneeded to estimate the degrees of freedom for which the energy needs to be maximizedand thereby ensure convergence on the desired saddle point. When combined with MOM,this approach can perform better than conventional SCF-MOM algorithms.
The numberof degrees of freedom for which the optimization needs to be carried out is, however, animportant consideration and so far the DO-MOM approach has only been formulated andimplemented in the context of the linear combination of atomic orbitals (LCAO) basis setwhere the number of degrees of freedom is relatively small.Real space grid (RSG) and plane wave (PW) basis sets have the advantage that thecomplete basis set limit can be reached by varying systematically a single parameter such asthe mesh spacing or plane wave cutoff. They are also more universal and can easily be appliedto diffuse states such as Rydberg and metallic states where a typical LCAO basis set needsto be supplemented by specially tailored diffuse basis functions. This is demonstrated for aRydberg state of NH in section 3 of this article. In calculations relevant to, for example,ultrafast experiments, the system may evolve through a series of localized and delocalizedstates, making it challenging to design an LCAO basis set that is complete enough for allthe relevant states. It is clearly more convenient to use an RSG or PW basis set in thesecases.In this article, a DO-MOM algorithm for calculating excited electronic states that can beused with both RSG and PW basis sets is presented. The calculations are variational and,since the Hellmann-Feynman theorem is satisfied not only at the minimum but also at anystationary point on the electronic energy surface, they can be used to evaluate the atomicforces, thereby providing a powerful tool for exploring excited-state potential energy surfaces(PESs) in, for example, simulations of the dynamics or minimum energy path calculations. The testsinclude charge-transfer excitations in nitrobenzene that are known to be challenging casesfor conventional algorithms. While the SCF-MOM algorithms show erratic behavior, likelybecause of the presence of several orbitals with similar energy, the DO-MOM calculationconverges in a robust way. Another test involves calculations of an electronic excitation ofthe CO molecule to degenerate states. Again, SCF-MOM shows erratic behavior while theDO-MOM calculation converges smoothly. There, the advantage of using complex orbitalsinstead of real orbitals is, furthermore, demonstrated.Two applications of the DO-MOM method are presented. The first one involves calcula-tions of metal-to-ligand charge-transfer (MLCT) and metal-centered (MC) excitations of the[Fe(bmip) ] (bmip=2,6-bis(3-methyl-imidazole-1-ylidine)-pyridine) complex, a prototypeof a class of Fe-based photosensitizers. This complex has been studied experimentallyusing X-ray emission and scattering with femtosecond resolution and theoretically usingTDDFT. In order to optimize its performance as a photosensitizer, an understandingof the transitions between the various states and the way they are affected by the ligandsis needed. An initial excitation to a singlet MLCT state is believed to be followed in partby a relaxation to a lower energy triplet MLCT state while another part decays to a tripletMC state where it generates a vibrational wavepacket along a metal-ligand bond stretchingcoordinate.
This branching occurs on ultrafast time scale, on the order of 100 fs, andinfluences the performance of the complex as photosensitizer. Dynamics of the molecule inthe lowest-lying, dark singlet MLCT state has been simulated using energy surfaces calcu-lated with TDDFT but a higher energy, bright singlet MLCT state is likely populated in theexperiments.
As a result, direct comparison with the ultrafast branching observed in theexperiment could not be made. Here, the DO-MOM method is used to calculate six excitedstates that are close in energy, including the bright singlet MLCT state along the metal-ligand bond stretching coordinate that is believed to be activated during the photoinduced5ynamics. The calculated excitation energy agrees well with the experimentally observedvalue and, furthermore, an estimate of the vibrational period in the lowest triplet MC stateis also found to be in good agreement with experimental observations. This demonstratesthat the DO-MOM method could, in future work, be used in dynamics simulations to helpinterpret the experiments on this and similar photosensitizer complexes.A second application of the DO-MOM method, which illustrates its applicability to ODDfunctionals, is a study of the effect of Perdew-Zunger self-interaction correction (PZ-SIC) on the estimated excitation energy in 13 transitions to singlet and triplet excited states in9 molecules. The results are compared with theoretical best estimates as well as experi-mental values. Also, O-H bond stretching curves in the water molecule and water dimerare calculated. It is found that PZ-SIC improves the shape of the curves, producing a localminimum analogous to what has been found in high level wave function calculations whilethe uncorrected functional does not.The article is organized as follows. In section 2, the DO-MOM algorithm for variationaldensity functional calculations of excited states is presented. Section 3 shows the resultsof the numerical tests. In sections 4 and 5, the applications of the DO-MOM method toexcited-state energy curves of the Fe II complex and the effect of PZ-SIC on excited states ofmolecules are presented, respectively. Finally, a discussion and conclusions are presented insection 6. In generalized KS-DFT, the energy of an electronic system is given by E [ Ψ ] = T s + Z d r ρ ( r ) v ext ( r ) + U [ ρ ] + E xc , (1)6here T s is the kinetic energy of a system of non-interacting electrons that have the samedensity ρ ( r ) as the interacting electron system T s = − M X i =1 X σ =0 , f iσ (cid:10) ψ iσ (cid:12)(cid:12) ∇ (cid:12)(cid:12) ψ iσ (cid:11) (2)and M X i =1 X σ =0 , f iσ h r | ψ iσ i h ψ iσ | r i = ρ ( r ) . (3)with occupation numbers { f iσ } . The occupation numbers can be chosen to be non-aufbau inorder to represent an excited state. M is the number of bands (orbitals) in the calculationsand σ is the spin index. v ext ( r ) is the external potential and U is the classical Coulombenergy U [ ρ ] = 12 Z Z d r d r ρ ( r ) ρ ( r ) | r − r | . (4) E xc is the XC energy, approximated in practice as a semilocal functional of the density and itsgradient, but can also include an explicit dependence on the orbitals (as in meta-generalizedgradient functionals and hybrid functionals).Excited states are obtained when the total energy is stationary with respect to the orbitals Ψ = ( | ψ i , · · · , | ψ M i ) T with non-aufbau occupation numbers and can correspond to saddlepoints. The electronic energy surface has dimensionality M N b , where N b is the numberof grid points in an RSG or number of PW coefficients. Even for a small molecule, thenumber of grid points can easily become large, on the order of . In order to facilitatethe stationary-point search problem, the orbitals Ψ are expanded in a linear combination ofsome auxiliary orbitals Φ Ψ = U Φ , (5)where U is an M × M unitary matrix. The auxiliary orbitals, Φ , can be chosen to be theground-state orbitals or any set of orbitals that represents an initial guess for the excited7tate. The energy can then be considered as a functional of both U and Φ E [ Ψ ] = E [ U Φ ] = F [ U, Φ ] . (6)Therefore, stationary points of E [ Ψ ] can be found in two steps: first by extremizing F withrespect to the expansion coefficients U and then by minimizing the functional L [ Φ ] = stat U F [ U, Φ ] (7)with respect to Φ : min Φ L [ Φ ] = min Φ stat U F [ U, Φ ] = stat Ψ E [ Ψ ] . (8)The introduction of the additional functional L [ Φ ] reduces the stationary-point search prob-lem into two simpler tasks. First, instead of finding a saddle point in a wave function space,one finds the saddle-point in the space of unitary matrices U of low dimensionality. Thereduction of dimensionality can further be achieved by decreasing the number of the virtualorbitals. For example, the first excited state of ammonia can be obtained by including only4 virtual orbitals and, therefore, the dimensionality of the problem equals 24, within a frozencore approximation. This is a significant simplification as compared to the original problemof finding a saddle point in a M N b dimensional space. Furthermore, an efficient algorithmbased on a recently proposed quasi-Newton method for finding saddle points in the spaceof unitary matrices can readily be used with minor modifications. The second simpler task isthe outer loop minimization where conventional energy function minimization algorithms generalized for a wave function optimization can be employed.The division of the original optimization problem into separate minimizations of differentdegrees of freedom is a standard technique employed for ground-state calculations of metals, self-interaction corrected density functional calculations and ensemble density functionalcalculations. The inner and outer loops are further described below. The optimization of8 is further performed in conjunction with the maximum overlap method (MOM) used todistribute the occupation numbers { f m } , similar to what has previously been done in thecontext of LCAO, and it is described in Appendix B.Additional considerations need to be addressed when a non-unitary invariant functionalis used as in the case of PZ-SIC. Semilocal approximations of the XC energy possess aspurious self-interaction error due to the inability of semilocal functionals to cancel the non-local self-interaction error in the classical Coulomb energy. Perdew and Zunger proposed anorbital-by-orbital correction E SIC [ Ψ ] = E [ Ψ ] − X iσ ( U [ ρ iσ ] + E xc [ ρ iσ , (9)making any approximation of the energy functional self-interaction free for one-electronsystems. The functional E SIC is not invariant with respect to unitary transformations of theoccupied orbitals and, therefore, an additional inner loop needs to be included in order tofind the optimal orbitals in the occupied subspace minimizing the self-interaction correctedenergy. In this case, the functional defined on the occupied subspace becomes unitaryinvariant.
Thus, finding a solution that corresponds to the SIC excited state is achievedin a three-loop optimization: stat Ψ E SIC [ Ψ ] = min Φ stat U min O E SIC [ U O Φ ] = min Φ stat U F SIC [ U, Φ ] (10)where the unitary minimization min is performed among occupied orbitals only. Moredetails on this additional inner loop are given in Appendix A.9 .1 Inner loop: Finding a stationary point of F [ U, Φ ] with respectto U . To find a stationary point of F [ U, Φ ] for a given Φ , an exponential transformation is made,analogous to what has previously been done for energy minimization in wave-function basedcalculations and density functional calculations, and, more recently, for saddle-point searches . The unitary matrix is parametrized as U = e A (11)where A is a skew-Hermitian matrix, A † = − A . The gradient of F with respect to theelements of A is ∂ F ∂A ∗ ij = 2 − δ ij (cid:20)Z e tA Le − tA dt (cid:21) ij , (12)where the elements of L are given by L ij = h φ j | ˆ h i | φ i i − h φ i | ˆ h j | φ j i (13)with ˆ h j defined as ˆ h j | ψ j i = δ E δ h ψ j | = f j (cid:20) − ∇ + ˆ v ext + ˆ v H + ˆ v xc (cid:21) | ψ j i (14)During the optimization, the elements of A are found iteratively using a limited-memoryversion of the symmetric rank one quasi-Newton algorithm. The initial inverse Hessian ispreconditioned with a diagonal matrix with elements K ij = 1 − (cid:15) i − (cid:15) j )( f i − f j ) (15)10here the (cid:15) i are the eigenvalues of the KS Hamiltonian. Since this preconditioner is validonly for the canonical representation of the Hamiltonian, the auxiliary orbitals are updatedto the canonical orbitals every X th iteration of the outer loop if the inner loop reaches amaximum number of iterations Φ ← C Ψ = CU Φ , (16)and set U = I (17)where C is the unitary matrix that transforms the auxiliary orbitals to the canonical orbitals.Further implementation details of the inner loop can be found in Refs. L [ Φ ] . Let M be a manifold in the Hilbert space such that M = { Φ : Z d r φ ∗ i ( r ) φ j ( r ) = δ ij , i, j = 1 . . . M } . (18)The tangent space to this manifold at Φ is defined as V Φ ( G ) = lim ε → ˆ R [ Φ + ε G ] − Φ ε = ˆ R ε [ Φ + ε G ] (cid:12)(cid:12) ε =0 (19)where ˆ R is the orthonormalization operator such that ˆ R [ Φ + ε G ] ∈ M , and G is a vector inthe Hilbert space. For example, ˆ R can be chosen as the Löwdin transformation. Let S X , Y be the overlap matrix between two vectors X , Y from the Hilbert space. Then S Φ + ε G , Φ + ε G = Z d r ( Φ + ε G )( Φ + ε G ) † = I + ε ( S Φ , G + S G , Φ ) + ε S G , G = I + ε W Φ , G + ε S G , G (20)11ith W Φ , G = S Φ , G + S G , Φ . For the Löwdin transformation ˆ R [ Ψ + ε G ] = S − / Ψ + ε G , Ψ + ε G ( Ψ + ε G ) (21)and therefore, the tangent space at Φ is V Φ ( G ) = G − W Φ , G Φ , (22)obtained after substituting eqs (20) and (21) into eq (19) keeping only first order terms withrespect to ε .The gradient of L can be calculated using the chain rule: δ L δ Φ ∗ = U δ E δ Ψ ∗ , (23)where (cid:18) δ E δ Ψ ∗ (cid:19) j = δ E δ h ψ j | (24)After defining the tangent space in equation (19) and the gradient in equations (23) and(24), the minimization of L can be written as Minimization algorithm. • Set k ← , choose initial guess Φ ( k ) ; calculate gradient G ( k ) = δ L /δ Φ ∗ ( k ) using eqs (23)and (24).• Project gradient on the tangent space V ( k ) = V ( G ( k ) ) at Φ ( k ) . ε and calculate residualerror ∆ ( k ) = k V ( k ) k .• While ∆ ( k ) > ε :1. Compute search direction P ( k ) according to the chosen minimization algorithmand apply preconditioning (for example, for gradient descent, P ( k ) = − G ( k ) and12nverse kinetic energy operator as preconditioner ).2. Project the search direction on the tangent space V ( k ) = V ( P ( k ) ) at Φ ( k ) .3. Choose optimal step length α ( k ) along V ( k ) and compute Φ ( k +1) ← Φ ( k ) + α ( k ) V ( k ) (25)4. Orthonormalize the wave functions, Φ ( k +1) ← ˆ R [ Φ ( k +1) ] .5. Compute new gradient G ( k +1) and project it onthe tangent space V ( k +1) = V ( G ( k +1) ) at Φ ( k +1) .6. Calculate residual ∆ ( k +1) = k V ( k +1) k .7. k ← k + 1 .• End.The search direction P ( k ) can be chosen using a conjugate gradient or a limited memoryquasi-Newton algorithm. Here, the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm as described in Ref. is used. For minimization, this algorithm is knownto give fast and robust convergence. The DO-MOM algorithm has been implemented in a development branch of the Grid-basedprojector augmented wave (GPAW) software and can be used with either a finite-differenceRSG or PW basis set. The calculations are carried out using the frozen core approximationand the projector-augmented wave method. The iterative SCF algorithms used here arebased on either the Davidson algorithm or the RMM-DIIS algorithm as implemented inGPAW. Both versions of the SCF algorithm make use of Pulay density mixing and MOM (see Appendix B). Default values of the convergence parameters are used. The Pulay density13ixing uses densities from three previous steps and the coefficient used in the linear mixingof the density with the density residual vector is 0.15. No damping of short-wavelengthdensity changes is used. In the DO-MOM calculations, Pulay density mixing is not used.Instead, the density is calculated from the orbitals obtained at each iteration. The gradientof the energy projected on the tangent space in the outer loop in DO, or the residual vectorin SCF, is preconditioned with the inverse kinetic energy operator.
In the outer loop of the DO-MOM algorithm, the search direction is calculated accordingto the L-BFGS algorithm, using only the previous step to estimate the Hessian matrix. Fora quasi-Newton algorithm a step length of 1 is a natural choice. However, the first step inthe optimization corresponds to the gradient descent algorithm and the following maximumstep length update is used in order to avoid too large changes in the orbitals: if the norm ofthe search direction ∆ = P † P (26)is larger than α max ( ∆ > α max ) then P ← α max ∆ P . (27)The value α max = 0 . is found to give reliable convergence. For the inner loop optimization,the limited-memory symmetric rank one update is used.The calculations were carried out in the following way if not stated otherwise: Themolecule is placed in a rectangular box with at least 7 Å vacuum space in all directions fromthe nuclei to the boundary of the box. Open boundary conditions are used. A grid meshspacing of 0.2 Å is employed. All the calculations are spin-unrestricted and use the PBEfunctional. The advantage of RSG over LCAO is illustrated in a calculation of the 3 p z Rydbergstate of NH for excitation from the HOMO, see Fig. 1. An LCAO calculation using theaug-cc-pVDZ basis set is clearly not sufficient to reproduce the RSG calculation while the14xpanded d-aug-cc-pVDZ basis set, which includes additional diffuse functions, gives closeagreement. w a v e f un c t i o n [ A n g / ] FD (h=0.2 Å )LCAO (aug-cc-vPDZ)LCAO (d-aug-cc-vPDZ) Figure 1: Comparison of three basis sets in a PBE calculation of the 3 p z excited state ofNH for excitation from the HOMO. The 3 p z orbital is shown along the rotational axis of themolecule ( z -axis). The finite-difference real-space grid calculation with grid spacing of 0.2 Å(red) is closely reproduced by an LCAO calculation when the expanded d-aug-cc-pVDZ basisset including an extra set of diffuse functions is used (blue dashed), but the aug-cc-pVDZbasis set is clearly not sufficient (black dashed).The energy of an open-shell singlet state is calculated using the spin purification: E s = 2 E ( ↑↓ ) − E ( ↑↑ ) , (28)where E ( ↑↑ ) is the energy of the triplet state and E ( ↑↓ ) is the energy of the mixed spin state.Both states are calculated independently and variationally. The singlet excited-state energywill hereafter be referred to as the energy calculated according to eq (28). This has beenfound to give a better estimate of the singlet excited-state energy in eDFT calculations usingsemilocal KS functionals compared to the estimate obtained from the mixed spin state. .1 Test I: Charge-transfer excitations in nitrobenzene Charge-transfer excitations in nitrobenzene are known to be challenging cases for conven-tional algorithms and are often used as benchmark tests.
The A ( π → π ∗ ) excita-tion transfers an electron from the benzene ring to the nitro group while in the A (n π → π ) excitation the transfer is in the opposite direction.A ground-state calculation including 9 virtual orbitals was first performed to obtain theinitial orbitals. The A ( π → π ∗ ) excited state was then calculated by promoting an electronfrom the HOMO-2 to the LUMO orbital, while the A (n π → π ) state was calculated bypromoting an electron from the HOMO-4 to the LUMO+1 orbital.An analysis of the performance of the DO-MOM calculation and comparison with thetwo versions of the SCF-based methods is presented in Fig. 2. Both the Davidson andRMM-DIIS implementations of SCF-MOM quickly approach the excitation energy of thetarget solution but then show erratic behaviour. This is attributed to the presence of severalorbitals with energy close to that of the orbital from which excitation occurs. The energydifference between HOMO-4 and HOMO-1 is only 0.56 eV and a change in the ordering of theorbitals occurs during the optimization of the excited state.
It is known that for orbitalsthat are energetically close, SCF algorithms have a difficulty converging unless smearingof occupation numbers is used or the parameters in the Pulay mixing are fine tuned. Incontrast, the DO-MOM algorithm shows robust convergence. Tight convergence is obtainedwithin 30 to 45 outer loop iterations as shown in Fig. 2(b) and (e). Initially, during eachouter loop iteration, several inner loop iterations are performed as shown in Fig. 2(c) and (f).Towards the minimum of the energy functional L , no inner-loop iterations are performed,only outer-loop iterations. 16 E n e r g y ( e V ) (a) SCF-Dav-MOMSCF-R-MOMDO-MOM 0 10 20 30 40Iteration No. (Outer Loop)876543210 L o g [ R e s i d u a l ] (b) SCF-Dav-MOMSCF-R-MOMDO-MOM 0 10 20 30Iteration No. (Outer Loop)0246810 N u m b e r o f I t e r a t i o n s ( I nn e r L oo p ) (c) DO-MOM20 40 60Iteration No. (Outer Loop)7.167.187.207.227.24 E n e r g y ( e V ) (d) SCF-Dav-MOMSCF-R-MOMDO-MOM 0 20 40 60Iteration No. (Outer Loop)86420 L o g [ R e s i d u a l ] (e) SCF-Dav-MOMSCF-R-MOMDO-MOM 0 10 20 30 40Iteration No. (Outer Loop)0246810 N u m b e r o f I t e r a t i o n s ( I nn e r L oo p ) (f) DO-MOM
Figure 2: Comparison of the performance of various algorithms in calculations of the A ( π → π ∗ ) excitation (a,b,c) and A (n π → π ) excitation (d,e,f) of nitrobenzene. TheSCF-MOM method, based either on Davidson (SCF-Dav-MOM) or RMM-DIIS (SFC-R-MOM) algorithms, does not converge well, while DO-MOM gives robust convergence. (a)and (d) show the convergence of the energy; (b) and (e) show the norm of the residual,and (c) and (f) show the number of inner-loop iterations for each outer-loop iteration in theDO-MOM calculation. 17 .2 Test II: Excitation to degenerate orbitals Degenerate electronic states need to be represented by multi-determinant wave functions.When a single determinant is used in a KS-DFT calculation two problems occur. The firstis a technical problem, as the SCF algorithm has difficulty converging unless large enoughsmearing of the occupation numbers is used. The second is a conceptual problem, in that dif-ferent single determinants which should in principle be degenerate can give different electrondensities and, as a result, different total energy. This occurs, for example, in calculations ofopen-shell atoms.
In the case of degenerate excited states, an additional problem arises.In order to unambiguously assign an excited state, the wave function needs to have thesymmetry of the excited state. With real-valued orbitals, which are most commonly em-ployed in electronic structure calculations, this requirement is not necessarily satisfied dueto symmetry breaking, as is demonstrated below.Consider the lowest valence excited state in carbon monoxide. Using a single determinant,this excited state can be described by the promotion of an electron from a σ orbital (ground-state HOMO) to one of the two lowest degenerate π ∗ orbitals, ( π ∗ x or π ∗ y , the ground-stateLUMO, in the case of real wave functions). For such an excitation, the SCF method withthe Davidson algorithm does not converge with integer occupation numbers. The energyoscillates around the excited-state solution, as shown in Fig. 3(a). The DO-MOM algorithm,however, gives smooth convergence to the excited-state solution.Since the chosen orbitals are real, they are not eigenstates of the z -component of theangular momentum operator, and the angular momentum of the single-determinant wavefunction around the internuclear axis ( z -axis) is not defined. In addition, the resultingelectron density lacks uniaxial symmetry. It has instead an elliptic shape in the x - y planewith orientation depending on which orbital is occupied, π ∗ x or π ∗ y , see Fig. 3(b). This isinconsistent with the symmetry of the molecule. In the DO-MOM calculation, the orbitalscan be chosen to be complex valued functions without any modifications of the algorithm.If the LUMO is chosen as a complex π ∗ +1 or π ∗− orbital, where +1 or -1 is the eigenvalue18f the z -component angular momentum operator, the single-determinant excited-state wavefunction has a well-defined angular momentum and can unambiguously be identified as a Π state. The solution DO-MOM converges to using complex orbitals has 0.15 eV higher energycompared to the real valued solution, but it is more accurate since the total density then hasrotational symmetry around the internuclear axis [see Fig. 3(b)]. Thus, the use of complexorbitals not only allows one to properly represent the total angular momentum of the excitedstate, but it also provides a density with the correct symmetry. The spin symmetry is stillbroken, however, in the unrestricted calculation.The importance of using complex orbitals in order to provide correct description ofthe ground state has been emphasized in calculations of atoms and molecules using self-interaction corrected functionals as well as within restricted Hartree-Fock theory andKS formalism. In particular, in the work of Lee et. al it was shown that real orbitalsbreak the cylyndrical symmetry of the density in the singlet ground state of O while complexorbitals restore such symmetry within the restricted KS formalism. Here, it is shown thata similar situation occurs in the excited states of open-shell singlets and that the symmetrycan be restored in the spin-unrestricted formalism using complex orbitals.19 E n e r g y ( e V ) (a) SCF-Dav-MOMDO-MOM 0 3 6 9 12 15 18 21Iteration No. (Outer Loop)6.506.556.606.656.706.75 E n e r g y ( e V ) (b) Real OrbitalsComplex orbitals
Figure 3: (a) Comparison of the performance of the SCF-MOM method implemented withDavidson and the DO-MOM method in calculations of the lowest n → π ∗ excitation of acarbon monoxide molecule. (b) Comparison of the DO-MOM calculations with real andcomplex orbitals. The inset shows the total electron density obtained with real orbitals(red) and complex orbitals (blue). The internuclear axis is perpendicular to the plane of theimage. The density obtained using real orbitals lacks the uniaxial symmetry. II carbene photo-sensitizer The first application of the DO-MOM method involves calculations of four MLCT and twoMC excited states of the [Fe(bmip) ] complex that consists of 63 atoms (see Fig. 4). Thecalculations are carried out with the BLYP functional. ] complex together with the orbitals involvedin the transition to the bright MLCT state (C atoms green, N atoms blue, H atoms white,Fe atom orange). (a) HOMO orbital of the ground state localized on the Fe atom. (b)LUMO+2 orbital of the ground state delocalized over the ligands. Promotion of the electronfrom HOMO to LUMO+2 corresponds to a charge-transfer excitation (denoted as MLCT state in the text). The isosurfaces are shown with grey and orange colors and correspond toisovalues of ± . Å − / .Fig. 5 shows the energy of the various states calculated along the metal-ligand bondstretching coordinate (the Q breathing normal mode according to Ref. ) that is believedto account for the nuclear dynamics following photoexcitation. The singlet state labelled MLCT in Fig. 5, corresponding to a HOMO-to-LUMO transition (see Fig. 4), has anexcitation energy of 2.58 eV, only 0.13 eV lower than the position of the maximum of theexperimental UV/VIS absorption spectrum of the complex dissolved in acetonitrile. Indeed,this state has the same character as the state with largest oscillator strength in TDDFTcalculations; thus, confirming that MLCT corresponds to the bright MLCT state witha calculated excitation energy in good agreement with experiment. The triplet with sameorbital occupancy (labelled MLCT in Fig. 5), the lower-lying singlet ( MLCT ) and tripletMLCT ( MLCT ) states arising from the HOMO-to-LUMO excitation are also shown, aswell as the lowest triplet MC state and the corresponding singlet with same character (arisingfrom HOMO-1 to LUMO+4 excitation). The orbitals involved in the transitions to these21tates are shown in the Supporting Information.
10 5 0 5 10Metal-Ligand-Stretching Normal Coordinate2.02.53.03.54.04.55.0 E n e r g y ( e V ) MLCT MLCT MLCT MLCT MC MC Figure 5: Potential energy curves along a metal-ligand normal coordinate (the breathingmode Q as defined in Ref. ) that is expected to be responsible for the nuclear dynamicsduring the relaxation from the initially photoexcited MLCT state to the MC state. MLCT:metal-to-ligand charge-transfer excitation, MC: metal-centered excitation.The combined ultrafast X-ray emission and scattering experiments have detected vibra-tional wavepacket dynamics along the metal-ligand stretching coordinate in a MC statewith a period of 278 fs. The curvature of the energy curve for the MC state calculatedwith DO-MOM is used to estimate the vibrational period, obtaining a value of 280 fs (seeSupporting Information). The shape of the PES predicted by DO-MOM with BLYP, there-fore, appears to agree well with the experiment in this respect, lending support for the use ofDO-MOM in future dynamics simulations to study this and other photosensitizer complexes.22
Application II: Assessment of Perdew-Zungerself-interaction correction
Calculations are carried out for 9 valence and 9 Rydberg excitations in 13 molecules involvingboth singlet and triplet excited states. The results of the calculations are compared totheoretical best estimates from Ref., which include corrections for basis set limitations and"all-electron" effects. The "all-electron" relaxation effects are estimated to be small, around . − . eV, in the present cases. On the other hand, the basis set correction can besignificant, and is important for a consistent comparison with the finite-difference RSG resultsobtained in the present work. The molecules are placed in a rectangular box with at least 9Å vacuum space in in all directions from any of the nuclei to the edge of the simulation box.This is found to be large enough to correctly describe even the diffuse Rydberg orbitals. Theatomic coordinates of the molecules are those given in Ref. The number of virtual orbitals isset to 8. The excitation energy is calculated with respect to the energy of the singlet groundstate. The excited states are calculated using the lowest energy single Slater determinant fora given transition without enforcing point-group symmetry constraints on the total density.The SIC calculations are carried out using the PBE functional with full PZ-SIC (PBE-SIC) and with the PBE functional where the correction is scaled by a half (PBE-SIC/2) assuch scaling has previously been found to provide better estimates of atomization energy ofmolecules and band gaps of solids.
The DO-MOM calculated values of the excitation energy are given in Table 1 for thetriplet states and in Table 2 for the singlet states. Since mixed spin states are often used inpractice as an approximation to a singlet energy surface, such calculations are also performedand the results are given in Table S1 in the Supporting Information.The calculations using the PBE functional give a mean error (ME) of -0.27 eV and a rootmean square error (RMSE) of 0.31 eV with respect to the theoretical best estimates for theexcitations to triplet states while a larger error is obtained for excitations to singlet states,23E of -0.46 eV and RMSE of 0.54 eV. If the spin purification for singlet states is not applied,the error in the excitation energy is significantly larger, with the RMSE being 0.95 eV.The calculations using self-interaction correction, i.e. the PBE-SIC functional, giveslightly better values, the magnitude of the ME with respect to the theoretical best esti-mates being 0.2 eV smaller for the singlet state excitations and 0.1 eV smaller for the tripletstate excitations, compared to PBE. The MAE is not improved as much. As has beenshown previously in ground state calculations, it is important to use complex orbitals inSIC calculations.
This is also found to be the case here in the calculated values of theexcitation energy. If real orbitals are used, the calculated values of the excitation energybecome worse than those obtained with the PBE functional (the RMSE being 0.44 eV fortriplet excitations, see Table S2 in Supporting Information).For the triplet excitations, the scaled self-interaction correction, PBE-SIC/2 functional,gives smaller improvement while for the singlet excitations the MAE and RMSE with respectto the theoretical best estimates is a bit smaller than for full correction PBE-SIC. For mixedspin states, PBE-SIC performs better (see Table S1 in Supporting Information).The small effect PZ-SIC has on the excitation energy is a result of cancellation of theestimated self-interaction energy in the ground and excited states. This orbital-by-orbitalestimate of the self-interaction energy, which is most appropriate for a single electron system,turns out to be of similar magnitude for the ground and excited states. This is a surprisingresult since the classical self-Coulomb energy of a diffuse, Rydberg orbital is known to besmaller than that of a more localized ground-state orbital. Table 3 shows an analysis of thisfor the water molecule. A near cancellation of the total self-interaction energy still occursbecause there is a simultaneous change in the self-XC term that offsets the difference in theclassical self-Coulomb energy. 24able 1: Energy of excitations to triplet states calculated with the DO-MOM method andcomparison with theoretical best estimates as well as experimental values. The calculationsmake use of a generalized gradient approximation Kohn-Sham functional (PBE), with scaledself-interaction correction (SIC/2) and with full self-interaction correction (SIC). The meanerror (ME), mean absolute error (MAE) and root mean square error (RMSE) are given withrespect to theoretical best estimates and with respect to experimental values at the bottomof the table. molecule excitation PBE SIC/2 SIC TBE a EXP b acetaldehyde 1 A ( n → π ∗ ; V) 3.65 3.75 3.79 3.98 3.97acetylene ∆ u ( π → π ∗ ; V) 6.33 5.89 6.04 6.40 6.0ammonia 2 A (n → Π (n → π ∗ ; V) 5.91 5.84 5.66 6.28 6.32diazomethane 1 A ( π → π ∗ ; V) 2.76 2.38 1.88 2.80ethylene 1 B u (n → B u ( π → π ∗ ; V) 4.46 4.63 4.75 4.54 4.36formaldehyde 1 B (n → A"(n → π ∗ ; V) 5.14 5.23 5.27 5.37 5.2hydrogen sulfide 1 A (n → B ( π → A"(n → π ∗ ; V) 4.20 4.35 4.41 4.64thioformaldehyde 1 A (n → π ∗ ; V) 1.71 1.81 1.88 1.941 B (n → A ( π → π ∗ ; V) 3.36 3.33 3.28 3.44 3.28water molecule 1 B (n → A (n → A (n → a Theoretical best estimates as given in Ref. b Experimental values listed in Ref. (see references therein). molecule excitation PBE SIC/2 SIC TBE a EXP b acetaldehyde A ( n → π ∗ ; V) 3.94 3.74 3.59 4.31 4.27acetylene ∆ u ( π → π ∗ ; V) 6.69 7.72 7.76 7.10 7.2ammonia 2 A (n → Π (n → π ∗ ; V) 7.48 7.96 9.36 8.48 8.51diazomethane 1 A ( π → π ∗ ; V) 2.94 2.56 2.14 3.13 3.14ethylene 1 B u (n → B u ( π → π ∗ ; V) 6.72 7.17 7.64 7.92 7.6formaldehyde 1 B (n → A"(n → π ∗ ; V) 5.38 5.17 5.01 5.63 5.8hydrogen sulfide 1 A (n → B ( π → A"(n → π ∗ ; V) 4.65 4.77 4.89 5.21thioformaldehyde 1 A (n → π ∗ ; V) 1.91 1.74 1.57 2.20 2.031 B (n → A ( π → π ∗ ; V) 5.36 6.02 6.66 6.34 6.2water molecule 1 B (n → A (n → A (n → a Corrected theoretical best estimates as given in Ref. b Experimental values listed in Ref. (see references therein). s orbital has smaller classicalself-Coulomb energy than the valence ground-state (GS) orbitals. However, the estimate ofthe total self-interaction energy is of similar magnitude for all the orbitals. orbital Coulomb XC Coulomb + XCGS a a s orbital b a Two orbitals from a ground-state calculation with PBE-SIC/2 (see Figure S5 in the SupplementaryInformation), giving different SIC estimates. b Orbital obtained from a DO-MOM excited-state calculation.
In order to further investigate the effect of PZ-SIC on excited states, the change inenergy as an O-H bond is stretched in a water molecule and a dimer of water molecules iscalculated for the ground and lowest singlet excited states. The spin-purification correctionis not applied in this case. First, the ground-state geometry is optimized until the maximumof the force on the atoms has magnitude below 0.01 eV/Å. Then, the hydrogen-bonded O-H bond is stretched by changing the position of the hydrogen atom in increments of 0.1Å, while keeping the positions of all other atoms frozen. The results of these calculationswith both PBE and PBE-SIC/2 are presented in Fig. 6 and compared with the results ofCR-EOMCCSD(T),ID/aug-cc-pVTZ calculations from Ref. The energy curves obtainedwith PBE and PBE-SIC/2 lie close to each other. However, the PBE-SIC/2 energy curvereproduces better the S-shape of the high level reference curves for both the monomer andthe dimer. For the water dimer in the lowest excited state, PBE-SIC/2 reproduces the localminimum at short bond length while PBE predicts a barrierless path towards the secondconstrained energy minimum near 1.8 Å O-H distance.27 .8 1.0 1.2 1.4 1.6 1.8 2.0OH distance ( Å )0246810 E n e r g y ( e V ) (a) Ref. a PBEPBE-SIC/2 0.8 1.0 1.2 1.4 1.6 1.8 2.0OH distance ( Å )0246810 (b) WaterDonor WaterAcceptor
Figure 6: Potential energy curves calculated for ground and lowest excited singlet statesusing mixed-spin determinant along the hydrogen-bonded O-H bond with all other moleculargeometry parameters fixed. (a) Water monomer, (b) Water dimer. The insets show theatomic configurations at the energy minimum and the endpoint of the scan in the O-H bondlength. a CR-EOMCCD(T),ID/aug-cc-pVTZ calculations from Ref. The RSG or PW representations have an advantage over LCAO in that full basis set limitcan be reached by systematically changing only a single parameter. However, such calcu-lations involve larger computational cost than LCAO and are limited to the frozen coreapproximation. A useful strategy for reducing the computational cost is to obtain initialorbitals from an LCAO calculation of the excited state and then switch to RSG or PWmode. All three types of representations can be available in the same software, as is the casewith our implementation in the GPAW software, making such a hybrid approach relativelystraightforward. 28he search of saddle points on multidimensional electronic energy surfaces is a consid-erably more difficult problem than minimization due to the fact that the energy must bemaximized along a few degrees of freedom that are not known a priori . Even if the orderof a saddle point is known, the success of the extremization of the energy depends on howclose the initial guess for the orbitals is to the target solution. The ground-state orbitals canrepresent a good enough initial guess for the iterative convergence on an excited state, butin some cases another excited-state solution provides a better choice. This can be especiallyimportant when the lowest energy excited state does not correspond to a HOMO-LUMOexcitation. For example, in the water dimer the lowest excited state corresponds to anelectron-hole pair localized on the hydrogen bond donor molecule but the HOMO-LUMOexcitation corresponds to the transfer of an electron from the donor to the acceptor moleculewith the hole localized on the donor and the excited electron localized on the acceptor. Inthe calculations of the potential energy curve for the lowest excited state in the water dimerpresented in Sec. 5, a two-step procedure was used: First, the HOMO-LUMO excitation wascomputed and the obtained solution then used as an initial guess in the calculation on thelowest excited state.In conclusion, we have developed a formulation of the DO-MOM algorithm for densityfunctional calculations of excited states that can be applied to RSG and PW representations,where it is not possible to include all virtual orbitals, unlike with LCAO basis sets. Thisimplementation of the DO-MOM algorithm is found to be robust and able to converge on ex-cited states that are challenging for commonly used SCF-MOM algorithms. The importanceof complex-valued orbitals in calculations of excited states with a KS-DFT functional is alsodemonstrated. In the case of the lowest excited state of CO, it was shown that real-valuedorbitals break the uniaxial symmetry of the electron density, while complex-valued orbitalsrestore the symmetry and correspond to higher energy solutions. The fact that complexorbitals provide a better description of the electronic system is in par with findings fromground-state calculations. ] complex shows that: (1) the approach is robust enough to allow the calculation of severalexcited states that are close in energy, including regions where the potential energy curvescross; and (2) the predicted excited-state properties are in good agreement with experimentalvalues. As for any single-determinant method, the quality of the excited states obtainedwith DO-MOM is expected to degrade when the states have a multi-determinantal characterand this can represent a limitation for the application of the methodology to dynamicssimulations. However, when explicit solvent effects are included in the simulations, thesymmetry breaking induced by the solvent can lift degeneracy, thereby reducing the presenceof multi-determinantal states. So far, excited-state dynamics simulations of systems as largeas transition metal complexes using a variational eDFT approach have been limited to a singleBorn-Oppenheimer surface. The present results represent a preliminary indication thatthe DO-MOM method is viable for nonadiabatic molecular dynamics simulations includingmultiple excited states. Current efforts in this direction include evaluating nonadiabaticcouplings between the excited states obtained variationally, using numerical or analyticalapproaches such as the one recently presented in Ref. A simulation of the excited-staterelaxation of the [Fe(bmip) ] complex starting from an excitation to the bright MLCTstate could help explain the experimentally observed ultrafast MLCT/ MC branching thataffects its performance as a photosensitizer.An assessment of self-interaction correction, the PZ-SIC, applied to the PBE functionalhas been performed for excited states. For valence and Rydberg excitations, the correctionimproves the value obtained for the excitation energy, but only by ca. better than results obtained with the PBE functional. Further improvementin excited-state energy calculations may require going beyond the self-interaction correctionbased on orbital densities, for example by including the complementary density in the errorestimate. An important advantage of the RSG and PW implementation presented here is thepossibility of doing calculations for extended systems that are subject to periodic boundaryconditions. This includes, for example, excited states of defects in crystals and solvent effectson electronic excitations of molecular complexes. This will be the topic of future studies withthe DO-MOM method.
Acknowledgement
This work was supported by the University of Iceland Research Fund and the IcelandicResearch Fund. AVI is supported by a doctoral fellowship from the University of Iceland.GL thanks the Icelandic Research Fund for financial support under grant number 196070.The authors are grateful to Mátyás Pápai for discussions on the Fe photosensitizer complex.The calculations were carried out at the Icelandic Research High Performance Computing(IRHPC) facility at the University of Iceland.31 ppendix
A Inner-Loop Minimizaion for ODD functional
The unitary matrix that minimizes an ODD functional, such as the PZ-SIC among occupiedorbitals is found using an exponential transformation: | ψ i i = X ij ( e − A ) ij | φ j i (A.1)where A is a skew-hermitian matrix, A † = − A . Let λ ij = D φ i (cid:12)(cid:12)(cid:12) ˆ h j + ˆ v j (cid:12)(cid:12)(cid:12) φ j E (A.2)with ˆ h j defined in Eq. (14) and ˆ v j = − δ ( U [ ρ j ] + E xc [ ρ j , δρ j (A.3)The gradient of the energy with respect to matrix elements of A is: ∂ E ODD ∂A ij = λ ij − λ ∗ ji + o ( k A k ) (A.4)This expression must be zero at the minimum and for equally occupied orbitals this leads tothe Pederson (or localization) conditions: h φ i | ˆ v j − ˆ v i | φ j i = 0 (A.5)The L-BFGS algorithm with inexact line search is used to find the optimal matrix A and, thereby, the optimal unitary matrix O = e − A . Details of the implementation of theexponential transformation are given in Ref. Maximum Overlap Method within the Projector-Augmented WaveApproach
Let {| ψ n i} be the N orbitals used as initial guess for an excited-state calculation and {| ψ ( k ) m i} the M orbitals at the k th iteration of the optimization. Within the PAW approach, theelements of the overlap matrix between orbitals {| ψ n i} and {| ψ ( k ) m i} are: S ( k ) nm = (cid:10) ψ n (cid:12)(cid:12) ψ ( k ) m (cid:11) = D ˜ ψ n (cid:12)(cid:12)(cid:12) ˜ ψ ( k ) m E + X a,i ,i D ˜ φ n (cid:12)(cid:12)(cid:12) ˜ p ai E (cid:16)(cid:10) φ ai (cid:12)(cid:12) φ ai (cid:11) − D ˜ φ ai (cid:12)(cid:12)(cid:12) ˜ φ ai E(cid:17) D ˜ p ai (cid:12)(cid:12)(cid:12) ˜ ψ ( k ) m E (B.1)where | ˜ ψ n i are pseudo orbitals, | ˜ p ai i ( | ˜ p ai i ) are projector functions localized on atom a ,while | φ ai i ( | φ ai i ) and | ˜ φ ai i ( | ˜ φ ai i ) are partial and pseudo partial waves localized on atom a , respectively. After computing the overlap matrix at iteration k , the occupation numbers { f m } of the orbitals {| ψ ( k ) m i} are chosen as following. An occupation number of 1 is given tothe first N orbitals with the highest numerical weights, evaluated from the projection ontothe manifold {| ψ n i} : ω ( k ) m = N X n =1 | S ( k ) nm | ! / (B.2)The remaining M − N orbitals are left unoccupied (occupation number of 0). References (1) Runge, E.; Gross, E. K. U. Density-functional theory for time-dependent systems.
Phys-ical Review Letters , , 997–1000.(2) Casida, M. E. In Recent Advances in Density Functional Methods ; Chong, D. P., Ed.;World Scientific, 1995; pp 155–192.(3) Dreuw, A.; Head-Gordon, M. Single-Reference ab Initio Methods for the Calculation ofExcited States of Large Molecules.
Chem. Rev. , , 4009–4037.334) Levine, B. G.; Ko, C.; Quenneville, J.; Martínez, T. J. Conical intersections and doubleexcitations in time-dependent density functional theory. Mol. Phys. , , 1039–1051.(5) Maitra, N. T.; Zhang, F.; Cave, R. J.; Burke, K. Double excitations within time-dependent density functional theory linear response. J. Chem. Phys. , , 5932–5937.(6) Tozer, D. J.; Handy, N. C. On the determination of excitation energies using densityfunctional theory. Phys. Chem. Chem. Phys. , , 2117–2121.(7) Oliveira, L. N.; Gross, E. K. U.; Kohn, W. Ensemble-Density functional theory forexcited states. Int. J. Quantum Chem. , , 707–716.(8) Yang, Z.-h.; Pribram-Jones, A.; Burke, K.; Ullrich, C. A. Direct Extraction of Excita-tion Energies from Ensemble Density-Functional Theory. Phys. Rev. Lett. , ,33003.(9) Deur, K.; Mazouin, L.; Fromager, E. Exact ensemble density functional theory forexcited states in a model system: Investigating the weight dependence of the correlationenergy. Phys. Rev. B , , 35120.(10) Hellman, A.; Razaznejad, B.; Lundqvist, B. I. Potential-energy surfaces for excitedstates in extended systems. J. Chem. Phys. , , 4593–4602.(11) Gavnholt, J.; Olsen, T.; Engelund, M.; Schiøtz, J. ∆ self-consistent field method toobtain potential energy surfaces of excited molecules on surfaces. Phys. Rev. B , , 75441.(12) Cheng, C.-L.; Wu, Q.; Van Voorhis, T. Rydberg energies using excited state densityfunctional theory. J. Chem. Phys. , , 124112.3413) Gilbert, A. T.; Besley, N. A.; Gill, P. M. Self-consistent field calculations of excitedstates using the maximum overlap method (MOM). J. Phys. Chem. A , ,13164–13171.(14) Kowalczyk, T.; Yost, S. R.; Voorhis, T. V. Assessment of the ∆ SCF density functionaltheory approach for electronic excitations in organic dyes.
J. Chem. Phys. , .(15) Levi, G.; Ivanov, A. V.; Jónsson, H. Variational Density Functional Calculations ofExcited States via Direct Optimization. Journal of Chemical Theory and Computation , , 6968–6982.(16) Hait, D.; Head-Gordon, M. Excited state orbital optimization via minimizing the squareof the gradient: General approach and application to singly and doubly excited statesvia density functional theory. J. Chem. Theory Comput. , , 1699–1710.(17) Carter-Fenk, K.; Herbert, J. M. State-Targeted Energy Projection: A Simple and Ro-bust Approach to Orbital Relaxation of Non-Aufbau Self-Consistent Field Solutions. J.Chem. Theory Comput. , , 5067–5082.(18) Mei, Y.; Yang, W. Excited-State Potential Energy Surfaces, Conical Intersections, andAnalytical Gradients from Ground-State Density Functional Theory. J. Phys. Chem.Lett. , , 2538–2545.(19) Ramos, P.; Pavanello, M. Low-lying excited states by constrained DFT. J. Chem. Phys. , , 144103.(20) Roychoudhury, S.; Sanvito, S.; O’Regan, D. D. Neutral excitation density-functionaltheory: an efficient and variational first-principles method for simulating neutral exci-tations in molecules. Sci. Rep. , , 8947.(21) Karpinski, N.; Ramos, P.; Pavanello, M. Capturing multireference excited states byconstrained-density-functional theory. Phys. Rev. A , , 32510.3522) Evangelista, F. A.; Shushkov, P.; Tully, J. C. Orthogonality Constrained Density Func-tional Theory for Electronic Excited States. J. Phys. Chem. A , , 7378–7392.(23) Ziegler, T.; Krykunov, M.; Cullen, J. The implementation of a self-consistent constrictedvariational density functional theory for the description of excited states. J. Chem.Phys. , , 124107.(24) Park, Y. C.; Senn, F.; Krykunov, M.; Ziegler, T. Self-Consistent Constricted Varia-tional Theory RSCF-CV( inf )-DFT and Its Restrictions to Obtain a Numerically Stable ∆ sCF-DFT-like Method: Theory and Calculations for Triplet States. J. Chem. TheoryComput. , , 5438–5452.(25) Davidson, E. The iterative calculation of a few of the lowest eigenvalues and corre-sponding eigenvectors of large real-symmetric matrices. J. Comput. Phys. , ,87.(26) Pulay, P. Convergence acceleration of iterative sequences. the case of scf iteration. Chem. Phys. Lett. , , 393–398.(27) Pulay, P. Improved SCF convergence acceleration. J. Comput. Chem. , , 556–560.(28) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calcu-lations using a plane-wave basis set. Phys. Rev. B , , 11169–11186.(29) Furthmüller, J.; Kresse, G. Efficiency of ab-initio total energy calculations for metalsand semiconductors using a plane-wave basis set. , , 15–50.(30) Garza, A. J.; Scuseria, G. E. Comparison of self-consistent field convergence accelerationtechniques. J. Chem. Phys. , .(31) Gillan, M. J. Calculation of the vacancy formation energy in aluminium. J. Phys. Con-dens. Matter , , 689–711. 3632) Payne, M.; Teter, M.; Allan, D.; Arias, T.; Joannopoulos, J. Iterative minimizationtechniques for ab initio total-energy calculations: Molecular dynamics and conjugategradients. Reviews of Modern Physics , , 1045–1097.(33) Hutter, J.; Parrinello, M.; Vogel, S. Exponential transformation of molecular orbitals. The Journal of Chemical Physics , , 3862–3865.(34) Marzari, N.; Vanderbilt, D.; Payne, M. C. Ensemble Density-Functional Theory for AbInitio Molecular Dynamics of Metals and Finite-Temperature Insulators. Phys. Rev.Lett. , , 1337–1340.(35) Ismail-Beigi, S.; Arias, T. A. New algebraic formulation of density functional calcula-tion. Comput. Phys. Commun. , , 1–45.(36) Van Voorhis, T.; Head-Gordon, M. A geometric approach to direct minimization. Molec-ular Physics , , 1713–1721.(37) VandeVondele, J.; Hutter, J. An efficient orbital transformation method for electronicstructure calculations. Journal of Chemical Physics , , 4365–4369.(38) Weber, V.; Vandevondele, J.; Hutter, J.; Niklasson, A. Direct energy functional mini-mization under orthogonality constraints. Journal of Chemical Physics , .(39) Freysoldt, C.; Boeck, S.; Neugebauer, J. Direct minimization technique for metals indensity functional theory. Phys. Rev. B , , 241103.(40) Ye, H.-Z.; Welborn, M.; Ricke, N. D.; Van Voorhis, T. σ -SCF: A direct energy-targetingmethod to mean-field excited states. J. Chem. Phys. , , 214104.(41) Levi, G.; Ivanov, A. V.; Jónsson, H. Variational calculations of excited states via directoptimization of the orbitals in DFT. Faraday Discuss. , , 448–466.3742) Ásgeirsson, V.; Jónsson, H. In Handbook of Materials Modeling: Methods: Theory andModeling ; Andreoni, W., Yip, S., Eds.; Springer International Publishing: Cham, 2020;pp 689–714.(43) Liu, Y.; Harlang, T.; Canton, S. E.; Chábera, P.; Suárez-Alcántara, K.; Fleckhaus, A.;Vithanage, D. A.; Göransson, E.; Corani, A.; Lomoth, R.; Sundström, V.; Wärn-mark, K. Towards longer-lived metal-to-ligand charge transfer states of iron(ii) com-plexes: an N-heterocyclic carbene approach.
Chem. Commun. , , 6412.(44) Harlang, T. C.; Liu, Y.; Gordivska, O.; Fredin, L. A.; Ponseca, C. S.; Huang, P.;Chábera, P.; Kjaer, K. S.; Mateos, H.; Uhlig, J.; Lomoth, R.; Wallenberg, R.;Styring, S.; Persson, P.; Sundström, V.; Wärnmark, K. Iron sensitizer converts light toelectrons with 92% yield. Nature Chemistry , , 883–889.(45) Lindh, L.; Chábera, P.; Rosemann, N. W.; Uhlig, J.; Wärnmark, K.; Yartsev, A.; Sund-ström, V.; Persson, P. Photophysics and Photochemistry of Iron Carbene Complexesfor Solar Energy Conversion and Photocatalysis. Catalysts , , 315.(46) Kunnus, K.; Vacher, M.; Harlang, T. C. B.; Kjær, K. S.; Haldrup, K.; Biasin, E.; vanDriel, T. B.; Pápa, M.; Chabera, P.; Liu, Y.; Tatsuno, H.; Timm, C.; Källman, E.; Del-cey, M.; Hartsock, R. W.; Reinhard, M. E.; Koroidov, S.; Laursen, M. G.; Hansen, F. B.;Vester, P.; Christensen, M.; Sandberg, L.; Németh, Z.; Szemes, D. S.; Bajnóczi, É.;Alonso-Mori, R.; Glownia, J. M.; Nelson, S.; Sikorski, M.; Sokaras, D.; Lemke, H. T.;Canton, S.; Møller, K. B.; Nielsen, M. M.; Vankó, G.; Wärnmark, K.; Sundström, V.;Persson, P.; Lundberg, M.; Uhlig, J.; Gaffney, K. J. Origin of vibrational wavepacketdynamics in Fe carbene photosensitizer determined with femtosecond X-ray emissionand scattering. Nature Communications , , 1–11.(47) Pápai, M.; Vankó, G.; Rozgonyi, T.; Penfold, T. J. High-Efficiency Iron Photosensitizer38xplained with Quantum Wavepacket Dynamics. J. Phys. Chem. Lett. , , 2009–2014, and private communication with the first author.(48) Pápai, M.; Rozgonyi, T.; Penfold, T. J.; Nielsen, M. M.; Møller, K. B. Simulation ofultrafast excited-state dynamics and elastic x-ray scattering by quantum wavepacketdynamics. The Journal of Chemical Physics , , 104307.(49) Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approxima-tions for many-electron systems. Phys. Rev. B , , 5048–5079.(50) Nocedal, J.; Wright, S. Numerical Optimization ; Springer, New York, 2006.(51) Stengel, M.; Spaldin, N. A. Self-interaction correction with Wannier functions.
Phys.Rev. B , , 155106.(52) Klüpfel, P.; Klüpfel, S.; Tsemekhman, K.; Jónsson, H. Optimization of functionals oforthonormal functions in the absence of unitary invariance. Lect. Notes Comput. Sci.(including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics) , , 23–33.(53) Lehtola, S.; Head-Gordon, M.; Jónsson, H. Complex orbitals, multiple local minima,and symmetry breaking in perdew-zunger self-interaction corrected density functionaltheory calculations. J. Chem. Theory Comput. , , 3195–3207.(54) Lehtola, S.; Parkhill, J.; Head-Gordon, M. Orbital optimisation in the perfect pairinghierarchy: applications to full-valence calculations on linear polyacenes. Mol. Phys. , , 547–560.(55) Borghi, G.; Park, C. H.; Nguyen, N. L.; Ferretti, A.; Marzari, N. Variational minimiza-tion of orbital-density-dependent functionals. Phys. Rev. B , .(56) Marzari, N.; Vanderbilt, D.; Payne, M. C. Ensemble density-functional theory for ab39nitio molecular dynamics of metals and finite-temperature insulators. Phys. Rev. Lett. , , 1337–1340.(57) Lehtola, S.; Jónsson, H. Variational, Self-Consistent Implementation of thePerdew–Zunger Self-Interaction Correction with Complex Optimal Orbitals. J. Chem.Theory Comput. , , 5324–5337.(58) Fern Rico, J. Á.; Paniagua, M.; Fern Alonso, J. I.; Fantucci, P. Restricted Hartree–Fockapproximation. II. Computational aspects of the direct minimization procedure. Journalof Computational Chemistry , , 41–47.(59) Rico, J. F.; De La Vega, J. M.; Alonso, J. I.; Fantucci, P. Restricted Hartree–Fockapproximation. I. Techniques for the energy minimization. Journal of ComputationalChemistry , , 33–40.(60) Douady, J.; Ellinger, Y.; Subra, R.; Levy, B. Exponential transformation of molecu-lar orbitals: A quadratically convergent SCF procedure. I. General formulation andapplication to closed-shell ground states. The Journal of Chemical Physics , ,1452–1462.(61) Head-Gordon, M.; Pople, J. Optimization of wave function and geometry in the finitebasis Hartree-Fock method. Journal of Physical Chemistry , , 3063–3069.(62) Ivanov, A. V.; Jónsson E.; Vegge, T.; Jónsson, H. Direct Energy Minimization Based onExponential Transformation in Density Functional Calculations of Finite and ExtendedSystems. arXiv:2101.12597 [physics.comp-ph] ,(63) Briggs, E. L.; Sullivan, D. J.; Bernholc, J. Large-scale electronic-structure calculationswith multigrid acceleration. Phys. Rev. B , , R5471–R5474.(64) Ref. 50, p. 121.(65) Ref. 50, p. 177. 4066) Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Dulak, M.; Ferrighi, L.; Gavn-holt, J.; Glinsvad, C.; Haikola, V.; Hansen, H. A.; Kristoffersen, H. H.; Kuisma, M.;Larsen, A. H.; Lehtovaara, L.; Ljungberg, M.; Lopez-Acevedo, O.; Moses, P. G.; Oja-nen, J.; Olsen, T.; Petzold, V.; Romero, N. A.; Stausholm-Møller, J.; Strange, M.;Tritsaris, G. A.; Vanin, M.; Walter, M.; Hammer, B.; Häkkinen, H.; Madsen, G.K. H.; Nieminen, R. M.; Nørskov, J. K.; Puska, M.; Rantala, T. T.; Schiøtz, J.; Thyge-sen, K. S.; Jacobsen, K. W. Electronic structure calculations with GPAW: a real-spaceimplementation of the projector augmented-wave method. J. Phys.: Condens. Matter , , 253202.(67) Mortensen, J.; Hansen, L.; Jacobsen, K. W. Real-space grid implementation of theprojector augmented wave method. Phys. Rev. B , , 035109.(68) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B , , 17953.(69) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation MadeSimple. Phys. Rev. Lett. , , 3865–3868.(70) Ziegler, T.; Rauk, A.; Baerends, E. J. On the calculation of multiplet energies by thehartree-fock-slater method. Theor. Chim. Acta , , 261–271.(71) Mewes, J. M.; Jovanović, V.; Marian, C. M.; Dreuw, A. On the molecular mechanism ofnon-radiative decay of nitrobenzene and the unforeseen challenges this simple moleculeholds for electronic structure theory. Phys. Chem. Chem. Phys. , , 12393–12406.(72) Becke, A. D. Current density in exchange-correlation functionals: Application to atomicstates. The Journal of Chemical Physics , , 6935–6938.(73) Johnson, E. R.; Dickson, R. M.; Becke, A. D. Density functionals and transition-metalatoms. The Journal of Chemical Physics , , 184104.4174) Klüpfel, S.; Klüpfel, P.; Jónsson, H. Importance of complex orbitals in calculating theself-interaction-corrected ground state of atoms. Phys. Rev. A , , 050501.(75) Klüpfel, S.; Klüpfel, P.; Jónsson, H. The effect of the Perdew-Zunger self-interactioncorrection to density functionals on the energetics of small molecules. J. Chem. Phys. , , 124102.(76) Lehtola, S.; Jónsson, E. Ö.; Jónsson, H. Effect of Complex-Valued Optimal Orbitals onAtomization Energies with the Perdew–Zunger Self-Interaction Correction to DensityFunctional Theory. J. Chem. Theory Comput. , , 4296–4302.(77) Small, D. W.; Sundstrom, E. J.; Head-Gordon, M. Restricted Hartree Fock usingcomplex-valued orbitals: A long-known but neglected tool in electronic structure theory. J. Chem. Phys. , , 24104.(78) Lee, J.; Bertels, L. W.; Small, D. W.; Head-Gordon, M. Kohn-Sham Density FunctionalTheory with Complex, Spin-Restricted Orbitals: Accessing a New Class of Densitieswithout the Symmetry Dilemma. Phys. Rev. Lett. , , 113001.(79) Becke, A. D. Density-functional exchange-energy approximation with correct asymp-totic behavior. Phys. Rev. A , , 3098–3100.(80) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energyformula into a functional of the electron density. Phys. Rev. B , , 785–789.(81) Loos, P.-F.; Scemama, A.; Blondel, A.; Garniron, Y.; Caffarel, M.; Jacquemin, D. AMountaineering Strategy to Excited States: Highly Accurate Reference Energies andBenchmarks. J. Chem. Theory Comput. , , 4360–4379.(82) Jonsson, H. Simulation of surface processes. Proc. Natl. Acad. Sci. , , 944–949.(83) Chipman, D. M. Stretching of hydrogen-bonded OH in the lowest singlet excited elec-tronic state of water dimer. J. Chem. Phys. , , 44305.4284) Levi, G.; Biasin, E.; Dohn, A. O.; Jónsson, H. On the interplay of solvent and confor-mational effects in simulated excited-state dynamics of a copper phenanthroline pho-tosensitizer. Phys. Chem. Chem. Phys. , , 748–757.(85) Levi, G.; Papai, M.; Henriksen, N. E.; Dohn, A. O.; Møller, K. B. Solution structureand ultrafast vibrational relaxation of the PtPOP complex revealed by ∆ SCF-QM/MMDirect Dynamics simulations.
J. Phys. Chem. C , , 7100–7119.(86) Dohn, A. O.; Jónsson, E. O.; Kjær, K. S.; B. van Driel, T.; Nielsen, M. M.; Jacob-sen, K. W.; Henriksen, N. E.; Møller, K. B. Direct Dynamics Studies of a BinuclearMetal Complex in Solution: The Interplay Between Vibrational Relaxation, Coherence,and Solvent Effects. J. Phys. Chem. Lett. , , 2414–2418.(87) Ramos, P.; Pavanello, M. Nonadiabatic couplings from a variational excited statemethod based on constrained DFT. The Journal of Chemical Physics , ,014110.(88) Lundin, U.; Eriksson, O. Novel method of self-interaction corrections in density func-tional calculations. Int. J. Quantum Chem. , , 247–252.(89) Pederson, M. R.; Heaton, R. A.; Lin, C. C. Local-density Hartree-Fock theory of elec-tronic states of molecules with self-interaction correction. J. Chem. Phys. , ,1972–1975. 43 raphical TOC Entry En e r g y Metal-ligand stretching coordinate
HOMO LUMO+2 charge-transfer excitation HOMO LUMO+2 upporting Information:Direct Optimization Method for VariationalExcited-State Density Functional Calculationsusing Real-Space Grid or Plane Waves Aleksei V. Ivanov, Gianluca Levi, Elvar ¨O. J´onsson, and Hannes J´onsson ∗ Science Institute and Faculty of Physical Sciences, University of Iceland, VR-III, 107Reykjav´ık, Iceland
E-mail: [email protected]
S-1 a r X i v : . [ phy s i c s . c o m p - ph ] F e b ontents ] S-32 Curvature of the MC potential energy curve of [Fe(bmip) ] S-63 Excitation Energy obtained from the Mixed-Spin Determinant S-84 The Perdew Zunger Self-Interaction Correction Restricted to Real Orbitals(rSIC) S-95 Ground-State SIC Complex Orbitals of the Water Monomer S-10References S-11
S-2
Ground-state molecular orbitals of [Fe(bmip) ] ] complex involved in the electronictransitions investigated in the present work. The isofurfaces are represented with gray andorange colors and correspond to isovalues of ± .
16 ˚A − / .Figure S1: HOMO-1S-3igure S2: LUMOS-4igure S3: LUMO+4S-5 Curvature of the MC potential energy curve of [Fe(bmip) ] MC state of [Fe(bmip) ] calculated with DO-MOM (BLYP) in the present work and linear response TDDFT (B3LYP*) from Ref. S1 arefitted using a fourth order polynomial: E ( Q ) = − d e + k Q − r e ) − α ( Q − r e ) + β ( Q − r e ) , (1)where k is the curvature. Below are the results of the fitting of the curve calculated withDO-MOM (BLYP): d e = − . k = 1 . × − r e = 5 . α = 1 . × − β = − . × − and of the curve calculated with TDDFT (B3LYP*): d e = − . k = 1 . × − r e = 5 . α = 9 . × − β = 2 . × − The fitting curves are shown in Fig. S4 together with the points obtained from the DO-MOM(BLYP) and TDDFT (B3LYP*) calculations. The curvatures of the fitting curves differ byaround 5%. The oscillation period estimated on the basis of TDDFT is 285 fs, S1 and thus,S-6e obtain an estimation of an oscillation period on the basis of eDFT as T eDF T = T T DDF T r k T DDF T k eDF T ≈
280 (fs)
10 5 0 5 10Metal-Ligand-Stretching Normal Coordinate2.02.53.03.54.04.5 E n e r g y ( e V ) DO-MOM-BLYP (FD=0.2)TDDTF-B3LYP (TZVPD) a Figure S4: Potential energy curves along a metal-ligand normal coordinate (the breathingmode Q as defined in Ref. S1 ) of the lowest MC state. The DO-MOM (BLYP) calculationswere performed using a finite-difference real-space basis set as described in the main text ofthe article. a Calculations from Ref. S1 S-7
Excitation Energy obtained from the Mixed-Spin De-terminant
Table S1: Mixed-Spin States species excitation PBE a SIC/2 b SIC b TBE c EXP d acetaldehyde 1 A ( n → π ∗ ; V) 3.79 3.74 3.69 4.31 4.27acetylene 1 ∆ u ( π → π ∗ ; V) 6.51 6.81 6.90 7.10 7.2ammonia 2 A (n → Π(n → π ∗ ; V) 6.70 6.90 7.51 8.48 8.51diazomethane 1 A ( π → π ∗ ; V) 2.85 2.47 2.01 3.13 3.14ethylene 1 B u (n → B u ( π → π ∗ ; V) 5.59 5.90 6.20 7.92 7.6formaldehyde 1 B (n → A”(n → π ∗ ; V) 5.26 5.20 5.14 5.63 5.8hydrogen sulfide 1 A (n → B ( π → A”(n → π ∗ ; V) 4.43 4.56 4.65 5.21thioformaldehyde 1 A (n → π ∗ ; V) 1.81 1.78 1.72 2.20 2.031 B (n → A ( π → π ∗ ; V) 4.36 4.67 4.97 6.34 6.2water 1 B (n → A (n → A (n → a Employing real-valued orbitals; b Initial guess for the wave functions is PBE real-valued orbitals followed by complexRudenberg-Edmiston localization; c Corrected theoretical best estimates as given in Ref.; S2 d Listed in Ref. S2 (see referencestherein); S-8
The Perdew Zunger Self-Interaction Correction Re-stricted to Real Orbitals (rSIC)
Table S2: Excitation energy from the PBE-rSIC/2 functional species excitation Mixed Spin Triplet Singletacetaldehyde 1 A ( n → π ∗ ; V) 3.71 3.73 3.69acetylene 1 ∆ u ( π → π ∗ ; V) 6.56 5.60 7.52ammonia 2 A (n → Π(n → π ∗ ; V) 6.94 5.79 8.09diazomethane 1 A ( π → π ∗ ; V) 2.26 2.15 2.36ethylene 1 B u (n → B u ( π → π ∗ ; V) 5.90 4.62 7.19formaldehyde 1 B (n → A”(n → π ∗ ; V) 5.01 5.06 4.96hydrogen sulfide 1 A (n → B ( π → A”(n → π ∗ ; V) 4.66 4.39 4.92thioformaldehyde 1 A (n → π ∗ ; V) 1.75 1.80 1.691 B (n → A ( π → π ∗ ; V) 4.72 3.15 6.29water 1 B (n → A (n → A (n → a ) -0.73 -0.38 -0.40ME (EXP b ) -0.59 -0.23 -0.24MAE (TBE a ) 0.73 0.39 0.46MAE (EXP a ) 0.59 0.26 0.34RMSE (TBE a ) 0.95 0.44 0.50RMSE (EXP a ) 0.88 0.31 0.41 a Deviation from corrected theoretical best estimates as given in Ref.; S2 b Deviation from experimental values collected inRef. S2 (see references therein); S-9
Ground-State SIC Complex Orbitals of the WaterMonomer
Figure S5: Real and imaginary parts of optimal ground-state SIC complex orbitals of thewater monomer. The isofurfaces are represented with blue and red colors and correspondto isovalues of ± .
05 ˚A − / . There are two degenerate pairs and only one orbital of eachdegenerate pair is shown. One optimal orbital (GS as denoted in Table 3 of the main text)is made of real part (a) and imaginary part (b). The other optimal orbital (GS as denotedin Table 3 of the main text) is made of real part (c) and imaginary part (d). Each figure(a), (b), (c), and (d) contains its own isosurfaces from two different views as shown by theaxis in insets. S-10 eferences (S1) P´apai, M.; Rozgonyi, T.; Penfold, T. J.; Nielsen, M. M.; Møller, K. B. Simulation ofultrafast excited-state dynamics and elastic x-ray scattering by quantum wavepacketdynamics. The Journal of Chemical Physics , , 104307.(S2) Loos, P.-F.; Scemama, A.; Blondel, A.; Garniron, Y.; Caffarel, M.; Jacquemin, D. AMountaineering Strategy to Excited States: Highly Accurate Reference Energies andBenchmarks. J. Chem. Theory Comput. ,14