Calculating Excited State Properties Using Kohn-Sham Density Functional Theory

The accuracy of excited states calculated with Kohn-Sham density functional theory using the maximum overlap method has been assessed for the calculation of adiabatic excitation energies, excited state structures, and excited state harmonic and anharmonic vibrational frequencies for open-shell singlet excited states. The computed Kohn-Sham adiabatic excitation energies are improved significantly by post self-consistent field spin-purification, but remain too low compared with experiment with a larger error than time-dependent density functional theory. Excited state structures and vibrational frequencies are also improved by spin-purification. The structures show a comparable accuracy to time-dependent density functional theory, while the harmonic vibrational frequencies are found to be more accurate for the majority of vibrational modes. The computed harmonic vibrational frequencies are also further improved by perturbative anharmonic corrections, suggesting a good description of the potential energy surface. Overall, excited state Kohn-Sham density functional theory is shown to provide an efficient method for the calculation of excited state structures and vibrational frequencies in open-shell singlet systems and provides a promising technique that can be applied to study large systems.


Introduction
Insight into the properties of electronic excited states and the reaction mechanisms that molecules undergo following photolysis is key to provide a full understanding of a wide number of processes involving the excited states of organic, transition metal and DNA/RNA related systems. 1 Theoretical calculations of molecular properties are well established, and the application of theoretical methods to describe the electronic ground state of molecules is often considered to be routine, although it can still be very demanding.The description of excited states presents a greater challenge due to problems such as the multiconfigurational nature of electronic excited states, the variety of different types of excited state, and the complexity that open shell systems with unpaired electrons can introduce.
Ideally it would be possible to treat excited states in a similar manner to the calculation of ground states.Excited states are commonly calculated using wavefunction based approaches such as complete active space self-consistent-field (CASSCF) coupled with multiconfigurational perturbation theory (CASPT2) 2 and multireference configuration interaction (MRCI) 3 .However, these wavefunction based approaches quickly become prohibitively expensive as the size of the system increases, particularly when optimizing a molecular structure within an excited state.Furthermore, the calculation of even more computationally demanding properties, such as accurate vibrational frequencies, [4][5][6][7] can also be of importance for the assignment of experimental vibrational spectra. 8,9 his motivates the search for alternative, computationally less expensive methods, to explore potential energy surfaces within electronically excited states.Kohn-Sham density functional theory (DFT) is computationally less expensive and can be applied to large systems.1][12][13][14][15][16][17][18][19][20] Results for adiabatic excitation energies, excited state structures, and vibrational frequencies for low-lying excited states of small molecules have previously been found to be in good agreement with gas-phase experimental data when using TDDFT. 10,12 owever, it is also possible to study electronically excited states directly using Kohn-Sham DFT without the need for TDDFT.2][23][24][25][26] One such approach, which is used in this work, is the maximum overlap method (MOM) that maintains the excited state within the Kohn-Sham calculation through an overlap analysis of the orbitals between iterations of the self-consistent field procedure. 24Recently, a constricted variational DFT approach was introduced that overcomes problems associated with excited state Kohn-Sham DFT methods, such as cases of unavoidable variational collapse, 27,28 and the relationship between excited state Kohn-Sham DFT and constrained (constricted) DFT has been described. 29Van Voorhis and co-workers have explored the theoretical foundation for such approaches to computing excited states and have shown that the states obtained through solution of the Kohn-Sham equations for a non-Aufbau occupation of the molecular orbitals correspond to stationary densities of the interacting systems within the adiabatic approximation to TDDFT. 30,31 here are several reasons why computing excited states using excited state Kohn-Sham DFT (also termed eDFT) 30 is an attractive alternative to TDDFT.
TDDFT is known to be inaccurate for particular types of electronic excitation when using common approximate forms of the exchange-correlation functional.Charge transfer, 32,33 Rydberg, 34,35 and core excitations 36 are all treated poorly, and while the use of range corrected functionals can lead to much better agreement with experiment, 36,37 within eDFT these types of transition require no special attention. 30eDFT is also computationally less expensive than TDDFT.Although the cost of TDDFT and eDFT scale similarly with system size, the prefactor for TDDFT is 2-3 times larger. 19Furthermore, unlike TDDFT, all of the functionality available for ground state DFT calculations, such as gradients, secondderivatives, and solvent modelling, are available for eDFT calculations in existing codes with little additional effort.
][40][41][42][43][44][45][46] For excited doublet and triplet states, directly applying eDFT can provide accurate excitation energies. 24,43 owever, computed excitation energies for some excited (open-shell) singlet states are found to be significantly underestimated. 24is failure can be attributed to the description of an open-shell singlet state with a single determinant, and the Zeigler post self-consistent field (SCF) spin-purification formula 47 provides a simple correction for this deficiency where E is the energy of the true singlet state, E S is the energy of the spin mixed (single determinant) state and E T is the energy of the corresponding triplet state.For some states, such as core excited states and Rydberg states 24,42 where E S ≈ E T the correction will be small.
However, for typical valence excited states E T < E S , resulting in an increase in the excitation energy.
While the accuracy of vertical excitation energies using eDFT has been studied, the reliability of eDFT to reproduce the shape of excited state potential energy surfaces has not.This is addressed in this paper, wherein we assess the accuracy of spin-mixed and spin-purified eDFT calculations for adiabatic excitation energies, excited state structures and vibrational frequencies compared to TDDFT for a set of small molecules where experimental gas-phase data are available.Furthermore, anharmonic frequencies provide a good indicator of potential energy surface quality, 4 and anharmonic frequencies based on second-order vibrational perturbation theory (VPT2) 48,49 are also assessed.

Computational Details
Constructing an extensive data set that incorporates excited state structures and vibrational frequencies is hindered by the relative scarcity of experimental gas phase data for excited states compared to the ground state, and by the variety of molecular excited states.The test set of molecules used here has been selected from those used to examine TDDFT excited state properties by Furche et al, 10,12 and includes only the lowest singlet excitations due to their importance in experimental photochemistry, and because they constitute a problematic case for eDFT calculations.For diatomic molecules, the data uses inferred harmonic frequencies for infinitesimal amplitudes of vibration, and anharmonically corrected bond lengths determined by Huber and Herzberg. 50cited state energies, optimized structures and harmonic vibrational frequencies were determined using TDDFT, spin-mixed excited state DFT, which we denote eDFT SM , and spin-purified excited state DFT, which we denote eDFT SP , for the B3LYP, 51,52 B97-1, 53 and EDF1 54 exchange-correlation functionals with the 6-311+G(d,p) basis set.The B3LYP, B97-1, and EDF1 exchange-correlation functionals were chosen for this comparison, as they are known reproduce ground state experimental frequencies well using scaled harmonic, 55,56 anharmonic, 4,57 and harmonic techniques, 4, 58 respectively.All calculations were performed within a locally modified version of the Q-Chem software package. 59The excited state Kohn-Sham calculations were performed using unrestricted Kohn-Sham DFT with the maximum overlap method invoked to prevent the variational collapse to the ground state to give the spin-mixed excited state. 24The molecular orbitals of the ground state were used as an initial guess for the excited state calculation, with a β-electron moved from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO).In the case of HCN the LUMO+1 was occupied.This gives the spin-mixed excited state, which we denote eDFT SM .Subsequently, energies, analytic gradients and Hessian matrices modified according to equation (1) were calculated to give the spin-purified state.These were calculated by combining energies, gradients and Hessian matrices obtained from separated SCF calculations on the singlet and triplet states.
Adiabatic excitation energies are calculated as the difference between the ground and excited state energies for each method.Zero point energy (ZPE) corrections were made (where stated) by subtracting the difference in harmonic vibrational ZPE between the ground and excited state from the uncorrected excitation energies.Analytical derivatives were unavailable for TDDFT with the B97-1 and EDF1 functionals, and for the aug-cc-pVQZ basis set.In these cases, numerical second derivatives of the energy with respect to nuclear displacements where evaluated from finite differences of the energies with a step size of 0.00189 a 0 .Anharmonic vibrational frequencies were determined using VPT2 and the transition-optimised shifted Hermite (TOSH) method. 49Third and fourth order derivatives of the energy with respect to nuclear displacements were calculated using finite difference schemes using analytical energies, gradients and Hessians using the default step size of 0.1 a 0 , 49 the resulting derivatives were used to compute a two-mode coupled quartic force field. 49en for small systems, this requires a large number of energy evaluations.Within the MOM framework, it remains possible that variational collapse to the ground state can occur, and this was found to be the case for the larger molecules studied here when evaluating the anharmonic frequencies.In order to ensure that the correct excited state energy is computed for all energy evaluations, if variational collapse is detected the energy evaluation is repeated with an alternative set of starting orbitals (the ground state orbitals for the minimum energy structure) or by converging a Hartree-Fock calculation for the excited state and using these orbitals as a starting point for the DFT calculation.Linear scaling of the computed harmonic frequencies was performed for comparison with anharmonic experimental data using values established for the 6-311+G(d,p) basis set.Following the work of Merrick et al, 56 B3LYP harmonic frequencies were scaled by 0.9688, B97-1 frequencies by 0.9684, and EDF1 frequencies by 0.9858, and we assume here that scaling factors derived for the ground state can be applied to excited states without further modification.Normal mode displacement diagrams have been generated using the Visual Molecular Dynamics software package. 60

Results and Discussion
The main focus of this work is to explore the performance of the eDFT SM and eDFT SP methods for open-shell singlet excitation energies, excited state structures and vibrational frequencies.However, it is necessary to use an appropriate basis set that ensures sufficient convergence with respect to basis set size.For this reason, excitation energies, geometries, and frequencies were evaluated for a selection of common basis sets using the 1 1 A'' state of H 2 CO as a test case, and the results are shown in Table 1.The results for the quadruple-ζ augcc-pVQZ basis set are sufficiently close to the basis set limit. 10 This basis set is too large to apply to moderately sized molecules and it is necessary to identify a smaller basis set that can adequately reproduce the results of the larger basis set at a lower computational cost.For the 1 1 A'' state of H 2 CO, the 6-31G(d) basis set shows large errors compared to aug-cc-pVQZ and augmentation with diffuse functions on non-hydrogen atoms is necessary for an accurate description of both energies and bond lengths.We have chosen the 6-311+G(d,p) basis set, in order to provide reasonable accuracy at a practical cost.Compared to aug-cc-pVQZ, 6-311+G(d,p) has errors of less than 0.01 eV for adiabatic excitation energies, with energies differing by as little as 0.001 eV for eDFT SM , 0.007 eV for eDFT SP , and 0.026 eV for TDDFT.Vibrational frequencies for the third harmonic mode, containing a significant amount of both carbon-oxygen and carbon-hydrogen motion, differs by 4 and 5 cm -1 from the aug-cc-pVQZ result for the eDFT SM and eDFT SP methods, and by less than 1 cm -1 for TDDFT.The bond length error is roughly 0.4 pm for the carbonyl bond and 0.3 pm for the carbon-hydrogen bonds for all three methods, and while bond lengths do show some improvement for the larger basis sets, this uncertainty is at a level comparable to that of experiment.

Excitation Energies
The calculation of accurate adiabatic excitation energies represents a greater challenge than for vertical excitation energies, since both a good description of the ground and excited state equilibrium geometries is required, in addition to the relative energy difference between them.Accurate ZPE corrections are further reliant on predicting the vibrational frequencies of both ground state and the excited state.Computed adiabatic excitation energies are shown in Table 2, in addition to the values from experiment.The excitation energies from TDDFT are closest to the experimental data and the mean absolute deviation (MAD) of the calculated values from experiment shows little variation between the different functionals used, with only a 0.04 eV difference in MAD, and a maximum difference between the two hybrid functionals of only 0.08 eV.The transitions studied here all involve excitation to low-lying valence states and are described well by standard hybrid and generalized gradient approximation functionals. 61 contrast, for both variants of eDFT, there is a much greater dependence on the functional used.This is quite surprising given that B97-1 and B3LYP both contain a similar fraction of Hartree-Fock exchange.It is known that excitation energies for valence excited states, like the ones considered here, have a large error when described by a single determinant.As expected, the excitation energies for eDFT SM are much too low compared with experiment with MAD values of 0.82 eV for B3LYP, 0.77 eV for B97-1, and 0.93 eV for EDF1, and showing errors in excess of 1 eV for BF, BH, CO, and SiF 2 .Application of the post SCF spin-purification leads to a significant reduction in the error with MAD values of 0.53, 0.38, and 0.58 eV for B3LYP, B97-1, and EDF1, respectively.Closer observation of the data shows that even for B97-1, which shows the best agreement with experiment, there is a consistent underestimation for the majority of excitation energies.This suggests that better agreement with experiment should be possible with a small modification to the functional.
However, it is interesting to note that some of the largest errors in the excitation energies occur for molecules with triple bonds.Such bonding can be described poorly by DFT and as such the source of the error for these excitation energies may lie with the ground state rather than the excited state.

Excited State Structures and Harmonic Frequencies
While excitation energies primarily examine the relative heights of the ground and excited state potential energy surfaces, the calculation of excited state structures and frequencies, depends on predicting the position of the energy minimum and the curvature of the potential energy surface at the minimum correctly.The predicted structures for the excited states are shown in Table 3.Previous work has shown that structures predicted from TDDFT are in good agreement with experiment. 10,12 he results in Table 3 show that the structures for eDFT SM and eDFT SP are very close to the values from TDDFT.For the HCP bond angle and Li 2 bond length, there is a significant error in eDFT SM that is partially corrected by spin-purification.In general, spin-purification leads to a small increase in the predicted bond lengths.There is some discrepancy with experiment for the dihedral angle (φ) in H 2 CS.For this angle, eDFT SM predicts a value that is closer to experiment than eDFT SP .However, the value for eDFT SP is consistent with the predictions of TDDFT raising some doubt over the experimental value.In contrast to excitation energies, the predicted structures are not strongly dependent on the choice of functional.
Unscaled harmonic frequencies are shown in Table 4. Calculated harmonic frequencies tend to overestimate their experimentally observed values for all three functionals, and indicate EDF1 to perform best due to systematic error cancellation.Frequencies calculated using the eDFT SM and TDDFT methods are in poor agreement with experiment, and have MAD values around 70 cm -1 for B3LYP and B97-1, with the eDFT MAD falling to around 55 cm -1 on spin-purification.As expected, on scaling of the harmonic frequencies, shown in Table 5, the agreement with experiment is significantly improved.The scaled MAD decreases by ~10-20 cm -1 compared to the harmonic value, although this improvement is more modest for the EDF1 functional.The results show that scaled harmonic frequencies for eDFT SM and eDFT SP are closer to experiment than those from TDDFT.The results when incorporating spinpurification have the best agreement with experiment, with a MAD of 35 cm -1 and 36 cm -1 for B3LYP and B97-1, respectively.
Comparison between the calculated and experimental harmonic frequencies for the diatomic molecules, shown in Table 6, shows a similar trend to the scaled polyatomic frequencies, with eDFT SM performing worst, followed by TDDFT, and with eDFT SP providing the best direct match with experiment, albeit with significantly larger errors.It is likely that functionals designed to predict harmonic frequencies correctly would lead to an improvement in the values. 62Overall comparison between experimental, harmonic and scaled frequencies indicate that the excited state methods have the following order of accuracy, eDFT SP > TDDFT > eDFT SM , for calculating the vibrational frequencies of open-shell singlet states.

Excited State Normal Mode Analysis
Comparison of the harmonic normal modes of molecules containing four atoms (see Figures 1-3 in the Supplementary Material), shows that while the frequencies vary significantly, the normal mode displacement vectors calculated for the excited state vibrations of H 2 CO, H 2 CS and C 2 H 2 using eDFT SM closely mirror the displacements calculated by TDDFT.The C-O stretching ν 3 mode for the excited state of H 2 CO is one notable exception, and is poorly treated by both methods.Here, the difference between the eDFT and TDDFT hydrogen atom displacement leads to reordering between the ν 3 and ν 4 modes using TDDFT compared with the eDFT order.While there is no assignment for the experimental atomic displacements in the excited state, the eDFT frequencies for the ν 3 mode are in good agreement with the experimental mode at 1183 cm -1 , with an error of less than 35 cm -1 for eDFT SM compared to over 100 cm -1 for TDDFT.Reordering of the ν 3 and ν 4 modes for TDDFT modes leads to a 58 cm -1 reduction (from 66 to 8 cm -1 ) in the error between the TDDFT result and experimental mode at 1293 cm -1 , with a corresponding increase in error for the 1183 cm -1 experimental mode, moving from 118 to 176 cm -1 .The corrected ordering of these two TDDFT modes has been used in Tables 1, 4 and 5. Despite resulting in significantly improved harmonic and scaled harmonic frequencies, the normal modes calculated using the eDFT SP method tend to involve the same atoms and displacement magnitudes as modes calculated using TDDFT and eDFT SM , but with rotated directions.

Anharmonic Vibrational Frequencies
Anharmonic frequencies evaluated using either the VPT2 or TOSH methods together with eDFT are shown in Table 7.For the majority of the vibrational modes, the computed frequencies lie close to their experimental values, suggesting that eDFT still provides a good description of the anharmonic potential energy surface around the equilibrium geometry.This is true for both the spin-mixed and spin-purified frequencies, and indicates that the eDFT SP formula in equation ( 1) remains accurate for greater sampling of potential energy surface despite providing a different description for some of the normal modes.These calculations show small deviations from experiment of below 30 cm -1 for most vibrational modes, which is reflected by a median absolute error of 23 cm -1 .There is a closer agreement with experiment for the spin-purified compared to the spin-mixed values, with MADs for the TOSH method of 41 and 36 cm -1 , respectively.Improvements are also seen even in cases where the harmonic surface is poorly treated within eDFT, such as with the excited states of SiF 2 and HCP, and the lower frequency mode of HCP is found to improve by ~100 cm -1 .The asymmetric C-H stretching mode of the H 2 CO excited state is also poorly treated.Overall, the results indicate that eDFT can provide an accurate description of both anharmonic as well as the harmonic frequencies, with the possibility of providing accuracy across a wider region of the potential energy surface.

Conclusions
In this paper the calculation of excited state properties for open-shell singlet states within excited state DFT has been investigated.Application of the spin-purification formula (eqn. (1)) leads to improved results for adiabatic excitation energies, excited state structures and vibrational frequencies compared to the spin mixed state.Spin-purification has the greatest effect on the excitation energies, which are much closer to experiment but remain consistently too low and less accurate than TDDFT.The computed excitation energies also show a greater sensitivity to the choice of exchange-correlation functional than TDDFT.The predicted structures from eDFT are of a similar accuracy to TDDFT, and spin-purification does correct the few significantly larger errors found with eDFT SM .For calculated scaled and harmonic frequencies, the values from eDFT SP are found to be closer to experiment than for TDDFT, although spin-purification can lead to different normal modes.Anharmonic eDFT frequencies provide accurate vibrational frequencies for the majority modes for both eDFT SM and eDFT SP .Consequently, eDFT represents a computationally inexpensive approach that provides an alternative to TDDFT that can be applied to study the properties and dynamics of electronically excited states of large molecules.

Table 1 .
Basis set dependence for some properties of the 1 1 A'' state of H 2 CO, calculated using B3LYP at the TDDFT, eDFT SM , and eDFT SP levels of theory.Adiabatic excitation energies (ΔE) are given in eV, bond lengths are in pm, and vibrational frequencies for the third harmonic mode (ν 3 ) are in cm -1 .

Table 2 .
Calculated ZPE corrected adiabatic excitation energies compared with experiment (in eV).

Table 3 .
Calculated excited state structural parameters compared with experiment.Bond lengths are in pm, and angles are in degrees.Experimental diatomic bond lengths are corrected for anharmonicity.

Table 4 .
Calculated excited state polyatomic harmonic frequencies compared with experiment (in cm -1 ).

Table 5 .
Calculated excited state polyatomic scaled harmonic frequencies compared with experiment (in cm -1 )

Table 6 .
Calculated excited state diatomic harmonic frequencies compared with harmonic experimental frequencies (in cm -1 ).