Identifying and tracing potential energy surfaces of electronic excitations with specific character via their transition origins: application to oxirane

Identifying and tracing potential energy surfaces of electronic excitations with specific character via their transition origins: application to oxirane

Jian-Hao Li    T. J. Zuehlsdorff    M. C. Payne    N. D. M. Hine TCM Group, Cavendish Laboratory, J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom
July 14, 2019

We show that the transition origins of electronic excitations identified by quantified natural transition orbital (QNTO) analysis can be employed to connect potential energy surfaces (PESs) according to their character across a wide range of molecular geometries. This is achieved by locating the switching of transition origins of adiabatic potential surfaces as the geometry changes. The transition vectors for analysing transition origins are provided by linear response time-dependent density functional theory (TDDFT) calculations under the Tamm-Dancoff approximation. We study the photochemical CO ring opening of oxirane as an example and show that the results corroborate the traditional Gomer-Noyes mechanism derived experimentally. The knowledge of specific states for the reaction also agrees well with that given by previous theoretical work using TDDFT surface-hopping dynamics that was validated by high-quality quantum Monte Carlo calculations. We also show that QNTO can be useful for considerably larger and more complex systems: by projecting the excitations to those of a reference oxirane molecule, the approach is able to identify and analyse specific excitations of a trans-2,3-diphenyloxirane molecule.

I Introduction

Kohn-Sham (KS) density functional theory (DFT)HK64 (); KS65 (); PY94 (); GV05 () can accurately reproduce many properties of molecular and solid state systemsKBP96 (); B97 (); CH08 () as compared to experimental methods and higher-level theory, e.g. molecular structures, vibrational frequencies, atomisation energies and ionisation potentials. Meanwhile, its time-dependent extension, TDDFTRG84 (); MUN06 (); C95 (); DH05 (), particularly in the linear response (LR) formalismC95 (), can deliver properties of excited states and optical absorption spectra of molecular systems. DFT and TDDFT calculations are increasingly able to treat large, complex biological systems, such as proteins, enzymes, and nucleic acids, due to the increasing availability of large computational facilities in recent yearsLC14 (); CC13 (); FP13 (); TB13 (). However, there remain serious obstacles to such uses: the intrinsic cubic scaling of traditional KS-DFT with respect to the number of electrons; the explosion of the size of the available configurational phase space associated with large, flexible molecules; and the challenges of accurate electrostatics of solvated systemsLCHHP (). When excited state properties are considered, further challenges come into play such as the proliferation of spurious charge transfer excitations, and the difficulty of identifying connected excited states on a PES.

Recent advances in linear-scaling (LS)-DFTWH94 (); WC06 (); SSF96 (); HG95 (); FG06 (); VK05 (); SA02 (); SHMP05 () have made the study of very large molecular systems much more feasible. LS-DFT relies on the exponential decay of the density matrix, , with respect to in a system with a finite band-gapBH97 (); IA99 (); HV01 (). In such systems, the density matrix, expressed using local orbitals, can be truncated by value or by range to render it local. In this work we use the ONETEP codeSHMP05 (), where a minimal number of local optimisable functions, called nonorthogonal generalised Wannier functions (NGWFs), are optimised in situ to best represent the density matrix of the system. Recently, Linear-Response TDDFT Tamm-Dancoff approximation (TDA) calculations have also become available in ONETEP, using a similar density matrix formalismZH13 (). Such methods can provide valuable knowledge about the electronic excitations of very large molecular systems such as proteins and DNA segments. However, the low symmetry and complex geometry of biological systems exacerbates the aforementioned problem of tracing the PES of specific excited states. This work therefore focuses on quantified natural transition orbital (QNTO) analysis, recently proposed by Li et al.LCGH11 (), which offers progress towards this goal.

QNTO analysis rewrites the transition vectors of an electronic excitation identified by a theoretical calculation, e.g. TDDFT, in terms of the natural transition orbitals (NTOs)M03 (); DH05 (); PWD14 (). The projection of hole (NTO-H) and electron (NTO-E) orbitals of the th NTO pair onto a standard orbital basis set such as natural bond orbitals (NBOs)RCW88 (); GLW12 (); LC12 () then helps quantify the transition origins of the NTOs of the excitation. Alongside the absorption energy and oscillator strength, knowledge of these transition origins helps provide a physical interpretation of the role of a given excitation dominated by single excitation character. Similar analysis may also be developed for single de-excitations, which is not as important in generic biological systems where de-excitations are usually negligibleDH05 (), signified by the good approximation of TDDFT/TDA to TDDFT. In addition, even in large and complex systems, electronic excitations may involve numerous electronic transitions between valence and conduction orbitals; using NTOs as the orbital basis, the number of dominant electronic transitions is often reduced to one or two, helping clarify what transition components contribute to the electronic excitations. QNTO analysis can help to assess the variation of transition origins determined by the different theoretical methods available in an electronic excitation calculation. For example, this approach has been used to assess the performance of different exchange-correlation functionals in predicting transition origins by comparing to a high-level theoretical calculationLCGH11 (). The variations of transition origins caused by the presence of different molecular environments and/or different molecular geometries can also be investigated. For example, QNTO has been used to study the role of the DNA backbone in mediating the transition origins of electronic excitations of B-DNALCGH12 ().

In the present paper, we introduce theoretical background to Linear-Response calculations in ONETEP (Sec. II), and show that QNTO analysis can be implemented with a non-orthogonal basis set such as the NGWFs used in ONETEP (Sec. III). In Sec. IV we then demonstrate two applications of LR-TDDFT/TDA/QNTO: (1) studying the photo-driven ring opening reaction of oxirane, and (2) identifying and comparing similar excitations between trans-2,3-diphenyloxirane and oxirane. We demonstrate that QNTO analysis can be readily used to locate state crossings that are otherwise unclear from plotting potential energy curves. The results are consistent with previous experimental and theoretical studies, showing that the present method is a feasible and reliable approach for investigating excited state related reactions.


Linear-scaling DFT as implemented in ONETEPSMHDP02 (); SHMP05 () relies on a reformulation of the single-electron density matrix in terms of local orbitals, as follows


where are single-particle eigenstates with occupation , while are spatially-localised nonorthogonal orbitals. is a generalised matrix representation of the occupancies, known as the density kernel. Note that implicit summation over repeated Greek indices is used throughout this work. The density kernel is defined by , with being the KS single-particle density operator. The nonorthogonal orbitals and their duals form a biorthonormal basis set such that .

ONETEP employs a minimal number of NGWFs, each expressed in terms of periodic bandwidth limited delta functions (psincs) centred on the grid points of the simulation cell, with a spacing determined by a kinetic energy cutoff akin to that in a plane-wave code.SHMP5757 () NGWFs are localised to atom-centred spherical regions defined by a cutoff . A total energy calculation proceeds by minimising the total energy simultaneously with respect to the matrix elements of the density kernel and the expansion coefficients of the NGWFs in terms of the psinc functions, via nested loops of conjugate gradients optimisationsMHSP03 (). NGWFs are initialised to pseudoatomic orbitals which solve the pseudopotential for a free atom.RHS12 () We confine our attention to closed-shell systems (even number of electrons) with non-fractional occupancy, and the total electronic energy of the real system is written as


where denotes ‘valence’, and with . The factor accounts for the spin degeneracy and is the standard “double-counting” correction. An FFT box techniqueSMHDP02 () is employed to allow one to construct matrix elements and density components in computational effort independent of system size.

As described in detail in Ref. HSMP08, , the minimisation is carried under the constraint of conserved total electron number,


where , and also the constraint of idempotent density matrix,


which is equivalent to requiring the orthonormal orbitals have occupancy of or . An important advantage of optimising the NGWFs in situ is that basis set superposition error (BSSE) is eliminated.HSMP06 () Sophisticated algorithms have been implemented so that a high performance can be achieved running on parallel computersHHMSP09 ().

Because the NGWFs are optimised in situ to describe the occupied part of the density matrix, they do not in general provide a good description of unoccupied states. Therefore, we instead generate and optimise in situ a second set of NGWFs, denoted as , obtained by minimising the bandstructure energy associated with a chosen range of low-energy conduction band states, having first projected out the valence states. The bandstructure energy to be minimised is written as


where denotes ‘conduction’, and . A full description of this procedure can be found in Ref. RHH11, .

Recently Zuehlsdorff et al. proposed a formulation for linear-scaling LR-TDDFT/TDA (hereafter simply denoted as TDDFT) calculations, which uses the conduction and valence NGWFs as a basis for describing the response densities of electronic excitationsZH13 (). This approach adapts the Casida formulation for finding transition vectors, , and excitation energies by using the TDDFT equation, , to a nonorthogonal local orbital basis. The quantity sought is the sum of the energies of the lowest excitations


A conjugate gradients minimisation of is performed subject to the constraint of keeping the eigenvectors orthonormal, expressing each transition using a transition density kernel

The transition density kernel, with , defines the transition density according to


and the Hartree-exchange-correlation self-consistent field matrix is defined by:

with being the Hartree and exchange-correlation response to the transition density as defined in Eq. 9 of Ref. ZH13, .

The minimisation of Eq. 6 proceeds as described in Ref. ZH13, . To avoid unwanted valence-valence transitions, we project out these unwanted components by constructing the kernel as


Additionally, we define a combined NGWF set by amalgamating valence and conduction NGWF sets to create a joint set, . Using this representation, and the corresponding overlap matrix , in place of and , means that we no longer need the explicit representation of the conduction space density operator. Instead, the conduction space projector can be replaced by subtracting the valence density operator from the identity operator within the subspace of the full joint NGWF representation. This defines through , giving


where and is the inverse of .

We also note that in the TDDFT the KS-system is linked to the real-system by


where and are the ‘exact’ many-body wavefunctions that are obtained from diagonalising the exact Hamiltonian in the configuration space. The linear dependence among the products becomes likely as the orbital basis approaches completion H86 () and the simple assignment may not be justified. Nevertheless, this assignment would become reasonableC95 () if the number of orbitals involved in is not too large, e.g. for low-lying excitations, and therefore the study of transition origins becomes possible via directly analysing or from the TDDFT calculation.

Iii Quantified natural transition orbital (QNTO) analysis

In this section, we first introduce the basic formulation of QNTO analysis in terms of an orthonormal orbital basis such as KS-orbitals. We then detail how it is implemented in ONETEP within the NGWF representation. The standard notation ‘’ is used throughout to denote conjugate transpose.

iii.1 Natural transition orbitals

The NTO basis


is an orbital basis that makes diagonal, with non-negative real numbers on the diagonal. In this representation the transition density operator is simplified to require summation over just one index:


denotes the number of non-zero values , which are ranked by magnitude. are called the singular values from the singular value decomposition (SVD) of . We also see that the and respectively diagonalise and ,


with the same eigenvalue set , where and are density matrices for the electron and hole orbitals, respectively.

In TDDFT calculations in ONETEP, the local orbitals (NGWFs) are not an orthonormal basis. If the NTOs are expressed in terms of NGWFs as


is then rewritten as and the target for decomposition becomes


If we hope to directly employ an SVD algorithm for orthonormal basis, we would have to first construct , where and . Once the singular vector matrices and are obtained, we would then need to calculate and in order to solve for and .

However, since the square root and inverse of a sparse matrix are in general much less sparse, their calculation would be very costly when dealing with large matrices. As an alternative approach, we can instead choose to diagonalise the and as shown in Eq. 13. Since and , the eigenvalue equations can be derived as


Note that in the SVD of , the relative phase between each pair of left (NTO-E) and right (NTO-H) singular vectors are fixed: if the left vector is multiplied by , the right vector must be multiplied by to remain a solution. Hence, the relative phase between each pair of eigenvectors needs to be determined after solving Eq. 16. The following relation from Eq. 15 can be used:


We define two matrices and as


From Eqs. 17 and 18, we see there is the relation:


We can thus start from NTO and proceed iteratively through higher NTOs, checking which of the associated two eigenvectors obtained in Eq. 16, denoted as and , produces the smaller value in Eq. 19. If Eq. 19 is not obeyed, we can multiply the whole vector (or ) by to obtain a correct pair of and . As this is a separate operation from the calculation of (Eq. 16), we can first solve for a number of eigenvalues , where , and subsequently define a number , where ) of these where the associated singular vectors are of interest. This can reduce the computational effort compared to directly running an SVD algorithm given that many excitations will be dominated by just NTO or a small number of NTO pairs.

iii.2 Interpretation of the transition origins of NTO electron-hole pairs

After obtaining the NTO-based transition vector element to any NTO, , next we want to interpret the transition origins of each NTO electron-hole pair from the projection to a standard orbital set. In the present work, we will in some places employ natural bond orbitals (NBOs)LC12 () as the standard orbital set. By such a projection the composition of NTO-H and NTO-E can be quantitatively defined. We generate NBOs by applying the NBO GB11 () program as a postprocessing step after the ONETEP calculation. A relatively large atomic orbital (AO) basis set is used in ONETEP for generating NBOs, so that its quality is comparable to that of standard in-situ optimised NGWFs. The electron spillage is thus reduced to a low level, typically in the systems studied here, while the HOMO energy is typically within eV of its value for NGWFs.

Note that the transition vector, , can be based on different sets of reference orbitals in different calculations. Rewriting the transition vector using the unique NTO basis removes this reference dependence. Projection of NTOs onto a common orbital set such as NBOs then provides a unification of interpretative orbital set for these NTOs. For example, excitations in two systems of different environments and/or geometries can be compared and the precise contribution of an orbital of interest to an excitation can be computed. We may be interested in the hole orbital involved in a excitation of a molecule, and can compute its quantitative involvement in an excitation of the same molecule embedded in a chemical/physical environment with or without a covalent bond formed.

Two NTOs can be compared using the root mean square deviation, , of the projection coefficients with respect to a common standard orbital setLCGH12 () such as the NBOs used in this work. The similarity of two NTOs can also be quantified by the magnitude of projection from one to the other. The projection sign can be arbitrary, as switching the overall phases of both NTO-H and NTO-E orbitals at the same time will result in them remaining a singular vector pair in Eq. 15. We will use the latter method for comparing two NTO-H(E) orbitals in the present work, while the projection to NBOs is only employed to help interpret the transition origins of NTOs.

The transition origins of each excitation can be seen as its fingerprint, helping to identify excitations of similar character. The adiabatic potential energy surfaces of electronic states can then be ‘reconnected’ on the basis of their character. If an electronic state of specific character changes its adiabatic energetic ordering from one geometry to the next, one would expect that at least one conical intersection (CX), key to many photochemical and photophysical processes BRS94 (); Y01 (); DY12 (); MBL00 (); LM07 (); CB03 (); LCM08 (); RC98 (); PD06 (); KT02 (); CO95 (), of electronic states has been passed between the two geometries. A by-product of analysing transition origins of a series of geometries along some coordinate(s) is, therefore, that one observes the influence of a CX (hyper-)point and locates the point where its projection onto the chosen path lies. Note that the current method is not by itself intended for locating the precise geometry of CX points, as to do that would require searching through the full internal coordinate space. In many well-defined cases, an appropriate series of geometries for analysing transition origins may be obtained by (constrained) geometry optimisation within the ground state or a chosen excited stateFA02 (), by molecular dynamics for the ground state or excited statesTTR07 (), via transition state search GPFKA03 () etc. Overall, the method can be highly illustrative as to how a photo-absorption driven process/reaction might proceed.

We also note that due to the approximate nature of practical TDDFT functionals, the behaviour of TDDFT excitation energies around CXs may be incorrectLKQM06 (), e.g. they may have overly rapid energy variation. TDDFT can also give the wrong dimension of the branching planeLKQM06 (). In other words, the PESs along some branching plane coordinates can be qualitatively wrong as the degeneracy is not lifted. While these issues of TDDFT on the study of CXs are still highly debated, TDDFT has, however, been able to offer a reasonable description of oxirane photochemistryTT08 (); CD07 () which will be our focus in the present work. It has also been shown that TDA as implemented in the current work has some advantage over full TDDFT in the description of regions around intersections HW14 (); CD07 (). It is known that the use of the TDA in TDDFT avoids the problems associated with singlet and triplet instabilitiesCG00 (). The former can occur for systems which have strong static correlation, signified by the vanishing HOMO-LUMO gap, whereas the latter can occur for systems which have an unstable closed-shell ground state, reflected by a negative triplet excitation energy within the TDDFT/TDA.

Overall, since locating CXs is not the point of the present method, the ability to describe the immediate vicinity of a CX is not critical. As long as the influence of a CX on the excitation character of geometries reasonably distant from the CX itself in configuration space is qualitatively correct, then the method outlined will be useful. Note also that the current QNTO method is not limited to use with TDDFT; any electronic structure methods that provide a decomposition of excitations into single transitions of electrons can be combined with this approach, and TDDFT is particularly useful for the study of very large molecular systems.

As an example, we can examine the nd excitation of an oxirane molecule using TDDFT/QNTO. The geometry of the oxirane molecule, with CCO angle=, is taken from a scan along the CCO angle which will be discussed in greater detail in Sec. IV. There are valence (occupied) KS-orbitals in the calculation, and therefore the HOMO(LUMO) is indexed () etc. The transition density operator expressed in the KS-orbital basis results from summation over several pairs of orbitals:


By contrast, in the NTO basis it is mostly accounted for by a single NTO pair:


The transition density operator is dominated by the NTO transition, which accounts for about of the total transition density, whereas in the KS-orbital basis the largest component, , only accounts for about of the density. In Fig. 1 we plot the NTO-H and NTO-E orbitals as well as the most dominant NBOs for each.

Figure 1: The NTO-H and NTO-E orbitals for the nd excitation of an oxirane molecule with CCO=. The st dominant NBO for each NTO-H(E) and the projection coefficient is shown below. Isovalue= e/bohr. The NBO labels specify bonding (BD) and anti-bonding (BD) and denote which pairs of atoms (eg C1 and O3) are involved in a given bond.

Iv Results and Discussion

In this section, we discuss two applications of QNTO. In subsection A, we focus on tracing the PESs associated with the ring opening process of oxirane. In subsection B, we use QNTO to locate and compare similar excitations between trans-2,3-diphenyloxirane and those of oxirane. As discussed in Sec. III B, it is interesting to see how environmental factors can affect the transition origins of excitations of a molecule, which in turn may modify its reaction pathways. The study of trans-2,3-diphenyloxirane thus serves as an excellent example of the use of QNTO to identify specific excitations of the oxirane molecule embedded in a different environment.

Figure 2: Gomer-Noyes mechanism. The underlying snapshots are from the CCO angle scans discussed in the text.
Figure 3: Constrained optimisation scan for different CCO angles followed by TDLDA/TDA calculation. All internal coordinates during the scan are optimised except for the CCO bond angle that is fixed in each step.

In all calculations a cut-off Coulomb approachHDHS11 () has been employed to avoid the influence of periodic images from neighbouring simulation cells. A minimal set of NGWF for H atom and NGWFs for C and O atoms is used. All calculations are well-converged with respect to the variational parameters of the localisation radii: suitable values are bohr for valence NGWFs and bohr for conduction NGWFs. The kinetic energy cutoff for the grid spacing is well-converged at eV. Norm-conserving pseudopotentials are used to replace core electrons. For DFT and TDDFT calculations we make use of the (A)LDA functional in the Perdew-Zunger parameterisationPZ81 () throughout. While checking the accuracy of results using other higher level method may be desirable, we note that DFT/TDDFT calculations with the standard (semi-)local functionals LDA and PBEPBE96 () are capable of describing the photochemically relevant PESs of ring-opening reactions in qualitative agreement with high quality quantum Monte Carlo calculationsTT08 (); CD07 (). In addition, even though the use of Hartree-Fock exchange hybridised functionals such as PBE0AB99 () or asymptotically corrected functionals such as LB94LB94 () may improve the DFT/TDDFT results, the results can deteriorate to be worse than that of pure PBE for structures near bond breaking region at large CCO angle where correlation effects become highly significantTT08 (). Although further examinations of functionals such as range-separated hybrid ones that have been used for similar molecular systems CH08 (); WPS09 (); SG11 (); LCGH11 () may still be desirable, the (A)LDA is considered to be sufficient for the purpose of the present work.

Figure 4: Projection results for different pairs of NTO-H(E) orbitals. Each row/column represents one of the first three excitations at each geometry at intervals in a scan in the range CCO=. Each pair of CCO angles form a block in the map, which is divided into cells for different pairing of excitations. The first and second number in a cell is the projection result (%) for NTO-H and NTO-E, respectively. A cell with yellow background indicates that both numbers are . The blocks marked by blue lines are where some character changes of excitations are observed one step off the diagonal blocks due to, for example, passing a CX. The small figure shows an example of adiabatic state behaviour near a CXAR97 (). C: (Diabatic) configuration coefficient for the two adiabatic states.

iv.1 Photoinduced ring opening of oxirane

CCO Excitation number
angle 1 2 3 4 5 6 7 8
60 0.9997 0.9999 0.9997 0.9981 0.9997 0.9985 0.9992 0.9961
65 0.9997 0.9997 0.9999 0.9980 0.9994 0.9986 0.9976 0.9980
70 0.9995 0.9997 0.9999 0.9976 0.9992 0.9916 0.9900 0.9983
75 0.9997 0.9998 0.9999 0.9964 0.9989 0.8737 0.8832 0.9947
80 0.9998 0.9999 0.9998 0.9921 0.9420 0.9874 0.9874 0.7348
85 0.9999 0.9998 0.9992 0.8237 0.7862 0.9972 0.9978 0.9560
90 0.9999 0.8798 0.8477 0.9971 0.9977 0.9926 0.9988 0.9789
95 0.9999 0.9800 0.9967 0.9963 0.9997 0.9954 0.9809 0.9996
100 0.9999 0.9843 0.9994 0.9996 0.9997 0.9859 0.9733 0.9877
105 0.9999 0.9868 0.9996 0.9997 0.9997 0.7387 0.7487 0.9853
110 0.9989 0.9503 0.9810 0.9067 0.9952 0.8161 0.7883 0.8233
115 0.9997 0.9959 0.9992 0.9993 0.9983 0.9909 0.9965 0.9877
Table 1: The transition vector component for NTO () of the first eight excitations along the oxirane CCO scan (Fig. 3).
Figure 5: Projection results for different pairs of NTO-H(E) orbitals. Each orbital comes from the first three excitations of a CCO= scan. An abrupt switch of transition origins betweeen the second and third excited state can be seen to occur between and .
Figure 6: The corresponding potential energy curves for the map shown in Fig. 5. The NTO-H(E) orbitals at and have also been plotted, with isovalue= e/bohr. While the NTO-H remains the same for the nd and rd excitation, the NTO-E orbitals can be seen to have switched their features between the nd and rd excitation at and , namely, the NTO-E character of the nd(rd) excitation at corresponds to that of the rd(nd) excitation at .
Figure 7: Projection results for different pairs of NTO-H(E) orbitals. Each orbital comes from the first three excitations of a CCO= scan. The green block indicates that, using the NTO-H and NTO-E of the first two excitations at CCO= as the reference set of orbitals, a switch of the dominating electron orbitals is observed between CCO=. For more details see the text.
Figure 8: The corresponding potential energy curves for the map shown in Fig. 7. Unlike the case in Fig. 5, the switch of dominant NTO character cannot be seen from the plot of NTO-H and NTO-E, reflected also by the map (Fig. 7) where the projection results among , , and cannot be seen to reveal a switch due to a smoother transition of dominant NTO character.
Figure 9: The st, nd, rd, and th singlet excitation potential energy curve (ranked at CCO=) that is re-connected based on the similarity of transition origins.
Figure 10: NBOs corresponding to Table 2. The lone pair on the oxygen (LP/O) accounts for a large proportion of the NTO-H for all three excitations. The st excitation is represented by NBOs denoted BD/C H, BD/C H, BD/C H, and BD/C H, whereas the nd excitation is represented by the same set of NBOs each with a reduced contribution, plus the more important BD/C H. The rd excitation is represented by BD/C O and BD/C O.

The Gomer-Noyes mechanism, as depicted in Fig. 2, was postulated and experimentally confirmedGN50 (); KIIT73 () as a way to explain the photoinduced ring opening process of the oxirane molecule. More discussions about this textbook molecule of stereochemistry can be found, for example, in Ref. MB90, . In particular, the specific electronic states involved in the Gomer-Noyes mechanism have been examined in some detail through TDDFT excited state molecular dynamics (MD) simulations validated by high-quality quantum Monte Carlo calculationsTT08 (); CD07 (), and several states have been found to be involved in mediating the ring-opening reaction.

Here, we perform a constrained ground state optimisation scan along the CCO bond angle of oxirane. We relax all internal coordinates except for a fixed CCO bond angle, which is set to values in the range to . For each oxirane geometry in the scan, a TDDFT calculation (vertical excitations) is performed followed by QNTO analysis to obtain the transition origins of the first few low-lying excitations.

The potential energy curves of the low-lying singlet and triplet states are plotted in Fig. 3. For reference we also show the negative HOMO energy and the value of the HOMO-LUMO gap at each angle. A caveat is that the calculation quality would deteriorate for those excitations of energies higher than the negative HOMO energy. As this tends to be underestimated by the TDDFT/LDACJCS98 (), one may obtain fewer bound states than are present in the real system. In addition, as mentioned in Sec. III. B, a vanishing HOMO-LUMO gap and triplet excitation energy signifies a (nearly) singlet and triplet instability, respectively, whereas employment of TDA helps circumvent these problems and may be recommendable when using generic approximate functionals in TDDFT calculations.

We note that from the excitation energies alone there is no clear indication of where crossings, indicating nearby CX points, have been passed while progressing along these curves. This is where the utility of QNTO analysis can be seen, as discussed in Sec. III B. If the transition density operator for these systems is dominated by the NTO component (Eq. 21), the character of each excitation can be represented by the transition origins of just NTO. The transition vector component for NTO, , along the scan are listed in Table 1. As can be seen, in most cases is higher than . For the first three excitations, the smallest , , is from the rd excitation of the oxirane geometry at CCO=, which still captures about of the density contribution. Thus, all the excitations on the oxirane geometries from the scan are dominated by the NTO transition, and when a switch of the NTO transition origins occurs, it shows that the excitation changes its character.

To put this in a form suitable for identifying electronic states of similar character section by section, we can construct a map (Fig. 4) showing the projections of the NTO pairs for each of the first three excitations (, , ) at each geometry onto their equivalents at different geometries, forming a block for each pair of geometries. We refer to the geometry listed down the row of the map as the ‘System’ (Sys) and compare it to a ‘Reference’ (Ref) geometry along the columns. To obtain a meaningful projection while accounting for the geometry change, the NTO-H(E) orbitals are translated from the Sys-geometry to the Ref-geometry and renormalised before being used for projection. The optimised atom-centred NGWFs from the calculation in the Sys-geometry are translated to the corresponding atom centres in the Ref-geometry to produce NGWFs . The NTO-H orbital (Eq. 14) is then re-expressed in terms of these using the original coefficients


The translated orbital is then renormalised:


and each block of projections is calculated as

for and .

N NBO type Proj NBO type Proj
1 LP O3 0.894 BD* C2 H7 -0.168
BD C2 H7 -0.215 BD* C1 H4 -0.167
BD C2 H5 0.213 BD* C1 H6 -0.165
BD C1 H6 0.199 BD* C2 H5 -0.164
BD C1 H4 -0.198 BD* C1 O3 0.055
BD* C2 H5 0.080 BD* C2 O3 -0.048
Total 0.976 Total 0.116
2 LP O3 -0.894 BD* C1 C2 -0.210
BD C2 H7 0.215 BD* C1 H4 0.124
BD C2 H5 -0.213 BD* C2 H5 -0.123
BD C1 H6 -0.200 BD* C1 H6 0.121
BD C1 H4 0.198 BD* C2 H7 -0.119
BD* C2 H5 -0.080 BD* C1 O3 0.019
Total 0.976 Total 0.104
3 LP O3 0.890 BD* C1 O3 -0.306
BD C2 H7 -0.219 BD* C2 O3 0.270
BD C2 H5 0.218 BD* C2 H7 -0.077
BD C1 H6 0.203 BD* C1 H4 -0.076
BD C1 H4 -0.202 BD* C1 H6 -0.067
BD* C2 H5 0.080 BD* C2 H5 -0.065
Total 0.977 Total 0.187
Table 2: Decomposition of NTO-H(E) into NBOs for the first three excitations of oxirane at CCO=. indicates the excitation number. Proj denotes the projection result. NBO type gives information on the NBO feature. For each excitation, the first ranked NBOs of lone pair (LB), bonding (BD), and anti-bonding (BD) types are listed. The highly basis-dependent Rydberg-type NBOs (RY) that accounts for most of these excitations are not listed. Note that the projection magnitude of the th-ranked NBO has gone down to below for all cases. Total means the total density contribution (square sum) of all the NBOs listed. These NBOs are plotted in Fig. 10.

In Fig. 4 one can immediately see which cells involve pairs of excitations with similar transition origins. We highlight those that have both projections greater than , indicating that the density contribution to one from the other is higher than . One then sees where there are switches of transition origins of these adiabatic excited state curves, or changes of transition character of all excited states, from one block to another. The former case can indicate that there has been one or more CXs passed between the two excitations. The latter would be caused by a character change of the ground state itself. This can be due to the fact that a CX between S and S has been passed, e.g. and , or that the molecular structures have become drastically different, e.g. and .

One can further subdivide the geometry scan to more accurately locate the position along the CCO angle coordinate where the character of the excitation changes. We demonstrate two such examples of switches of excitation which show subtly different behaviour. Fig. 5 shows an example where the switch of excitation character previously shown from the coarser map to occur between and is narrowed down to occurring abruptly between and . The corresponding potential energy curves for the two states involved, and the NTO-H and NTO-E orbitals before and after the switch, are plotted in Fig. 6, showing that the bands approach very closely around this point.

In other cases, however, when a more closely spaced geometry map is considered, the change of excitation character may be so smooth that no abrupt character switch occurring between two minimally-differentiated geometries can be found. An example is shown in Fig. 7. The two adiabatic excited state curves as plotted in Fig. 8 also do not come closer than eV to each other. However, as the transition characters have switched for the two adiabatic excited states curves between CCO=, there must be a CX located nearby, i.e. on other (internal) coordinates, even if it does not occur exactly along the scanned coordinate. Thus, without analysing the electronic character, only from the plot of the energy, one would not necessarily know that molecular motions in this region of the scanned coordinate would be near a CX, and that an electronic state change from the upper surface to the lower surface could occur. However, note that from considering the overall shape of the energies in Fig. 3, it is not difficult to envisage that a crossing could have occurred between S and S. For further analysis, the NTO-H(E) orbitals of excitations from a geometry somewhat further away from those being considered can be used as a reference set, similar to a diabatic configurations set used for locating CXsAR97 (). For the map shown in Fig. 7, the NTO-H and NTO-E orbitals of the first two excitations at CCO= can provide an appropriate reference set. We thus project the NTO1-H and NTO1-E orbitals at varying geometries onto these reference orbitals, giving:


This shows that the two reference electron orbitals from CCO= switch domination between CCO= and , which could not have been determined from Fig. 8 alone. Thus, a better strategy without overlooking a character switch between two curves would be to employ a wide coordinate window within which all the cross-projection data between each pair of geometries are checked.

We note that the exact point where the domination switches occur can be somewhat affected by the choice of reference orbitals used, although this does not in general prevent us from tracing down the CX to a small region in the scanned coordinate. In the language of diabatic configurationsAR97 (), the situation can be analysed by considering the matrix elements of a Hamiltonian acting on a reduced subspace consisting of just two such configurations, denoted as and . The switch of domination of two diabatic configurations as a given geometry path is traversed indicates that the surface has been crossed, though the exact shape of the surface of will be affected by the reference diabatic configurations used.

Evidence of a CX between S and S, on the other hand, would be reflected in the map by the fact that a whole column of excitations would change transition character at the same time, because the ground state itself on which excitations are based has changed character. We are able to locate two such points in Fig. 4: between and , and between and , both related to the change of the ground state character expected from the Gomer-Noyes mechanism.

After locating the regions where two adiabatic excited state curves switch their characters, we can plot the adiabatic singlet states again, with the curves for the various excitations ‘re-connected’ based on their transition origins, as shown in Fig. 9. State labels of (n,), (n,) and (n,) are based on the symmetry at the CCO= geometry and propagated to other geometries based on the connection of the states.

It then becomes clear that population of the S state at CCO= can directly lead to a crossing between S and S. The energy decreases monotonically along with the widening CCO angle, which will result in a rapid ring opening reaction. This is consistent with the previous excited states MD studyTT08 () which found that the population of the (n,) state would facilitate the occurrence of ring opening process. If the (n,s) state is populated instead, it has to overcome an energy barrier to a widening of CCO angle before meeting a crossing with the (n,) state which can then lead to the ring opening of oxirane. Thus, this chemical process is more difficult and would produce a lower quantum yield. We identify that the ground state is re-populated from the (n,) state after oxirane passing the critical point between and .

Meanwhile, passing the critical point between and leads to another change of ground state character, which corresponds to a proton transfer process producing acetaldehyde, the product in the rd stage of the Gomer-Noyes mechanism (Fig. 2). This is also consistent with the MD studyTT08 () finding that the proton shift occurs after the oxirane has been de-excited to the ground state. Under thermally elevated conditions this proton shift will likely occur at lower CCO angles: here we address only the quasistatic constrained geometries. A follow-up reaction breaking the C-C bond, the final stage in Fig. 2, would then be possible while remaining in the ground state, given an appropriate combination of internal motionTT08 ().

To demonstrate the chemical composition of these transitions, a decomposition of NTO-H and NTO-E into natural bond orbitals for the first three excitations at CCO= is shown in Table 2. The AO basis used for generating NBOs is fairly large, producing NBOs, but these are mostly of Rydberg (RY) type. In practice, the lone-pair (LP) and bonding (BD) NBOs describe the majority of the NTO1-H excitation, which in all three cases is well-described by just five main NBOs. For the NTO-E orbitals of the three excitations, however, a substantial contribution comes from RY orbitals, which are highly dependent on the AO basis set used in generating NBOs. The LP and BD NBOs for NTO1-H, and the most relevant BD NBOs for NTO1-E are shown in Fig. 10.

iv.2 Trans-2,3-Diphenyloxirane

Figure 11: A plot of the system molecule, trans-2,3-diphenyloxirane, and the reference oxirane molecule that is extracted from the trans-2,3-diphenyloxirane. The five atoms labelled are the common (core) atoms.
Figure 12: Plots of the NTO-H(E) orbitals for the first eight Ref-excitations as well as several Sys-excitations that have a substantial overlap to some Ref-excitations. Isovalue= e/bohr for NTO-H; for showing more details of the NTO-E the isovalue is set to e/bohr. See Table 3, 4 and the text for more discussions.

In this subsection we focus on a significantly larger system, the trans-2,3-diphenyloxirane molecule, which contains the same CCO subsystem as oxirane as its central motif. We show that we can use the QNTO method to compare the excitations with those of an oxirane molecule. As would be expected, molecules with different physical/chemical environments can have distinct properties. For example, aryl substitution of oxirane can lead to a different ring opening mechanismFF09 (). The fact that the trans-2,3-diphenyloxirane, referred to as the system (Sys) molecule, is formed by two aryl substitutions on oxirane makes it a good example for examining the influence of chemical environment on the excitations of the core oxirane moiety using TDDFT/QNTO.

The structure of the Sys-molecule is first fully optimised without constraints. A equivalent geometry for the oxirane molecule, referred to as the reference (Ref) molecule, is then extracted from the optimised Sys-molecule by replacing the two aryl substitutents with hydrogen atoms and optimising the C-H bond lengths, at fixed bond angles (Fig. 11). The widespread problem that the energy of charge-transfer excitations is underestimated by local density functionals such as the LDADH05 () means in this case that we expect to find large numbers of charge-transfer excitations in the energy window of interest. In practice we find we need to calculate up to excitations of the Sys-molecule so that the highest (th) excitation energy calculated, eV, exceeds the energy of the th excitation of the Ref-molecule ( eV), which was the highest Ref-excitation studied in the current work.

The NTO-H and NTO-E orbitals for the excitations of Sys-molecule (Sys-excitations) and Ref-molecule (Ref-excitations) are then compared to each other to see if similar excitations can be identified between the first Ref-excitations and the first Sys-excitations. This is achieved by constructing, for each NTO-H(E), a renormalised orbital based on just the component of each on the local orbitals of the atoms common to both geometries (the oxirane core). We denote the NTO-H(E) orbital of a Sys-excitation as , and the corresponding core orbital cNTO-H(E) in the form . Similarly, we define and for the Ref-molecule. The orbtitals are the NTO-H(E) orbitals of the original reference oxirane geometry, whereas , , and are constructed via a similar procedure to that used in constructing the projections of the previous section.

In particular, constructing and requires a translation of the NTO-H(E) from the Sys-molecule geometry to the Ref-molecule geometry. Using NTO-H as an example, the set below are constructed:


followed by renormalisation as in Eq. 23 to obtain , , and , respectively. The “sr” in Eq. 26 denotes that the orbital is translated from the Sys-molecule to the Ref-molecule. and are defined by setting respectively the coefficients and in Eq. 14 to zero for those NGWFs not centred on the core atoms; the subscript ‘s’ for Sys-excitations and ‘r’ for Ref-excitations indicate that the original NGWF coefficient matrices have been used.

By evaluating the overlap between an NTO-H or NTO-E orbital of the Sys-molecule and the corresponding NTO of the Ref-molecule, one can then tell whether there is a strong overlap of the character in the core region. In addition, the overlaps and allow one to determine how large a proportion of the total excitation is accounted for by the component within the core region.

In the general cases where excitations are dominated by the NTO, one can then say that the two excitations from different molecules are similar if the following three conditions hold: (1) a Sys-excitation must have a substantial for both NTO-H and NTO-E, (2) a Ref-excitation must have a substantial for both NTO-H and NTO-E, and (3) the magnitude of must be large.

We find that of all the excitations of the larger system, there are whose cNTO-H and cNTO-E respectively match with the cNTO-H and cNTO-E of one or more Ref-excitations. We choose a threshold of for both cNTO-H and cNTO-E. Overall there are matching combinations, as some excitations overlap significantly with several excitations of the reference system. The st, nd, th, and th Ref-excitation are involved in these matching pairs.

If the threshold is lowered to , the number of involved Sys-excitations becomes and the number of matchings becomes . The rd and th Ref-excitation also now match with some Sys-excitations. All the excitation matchings are listed in Table 3.

Ref 1 2 3 4 5 6 7 8
(96;76) (96;72) (96;99) (96;78) (96;98) (96;98) (96;51) (00;75)
Sys 15 34 26 6 15 16
(80;38) (78;81) (79;32) (81;15) (80;38) (32;39)
(-00;73) (-00;-87) (-00;65) (-99;61) (-00;-59) (-90;74)
Sys 24 43 39 14 24 21
(48;39) (81;34) (51;37) (49;19) (48;39) (04;39)
(-00;74) (00;84) (-00;68) (97;65) (-00;-60) (80;73)
Sys 37 44 59 34 37 31
(81;65) (65;87) (39;43) (78;81) (81;65) (50;24)
(00;-75) (-00;-86) (79;60) (00;83) (00;63) (-92;-65)
Sys 52 58 (81) 43 52 41
(47;66) (13;39) (51;62) (81;34) (47;66) (32;63)
(-00;-75) (-72;81) (00;-58) (-00;-76) (-00;63) (90;-78)
Sys 57 44 57 46
(61;38) (65;87) (61;38) (04;65)
(99;72) (00;83) (99;-60) (75;-78)
Sys 65 58 79 64
(28;64) (13;39) (42;53) (29;65)
(-66;68) (72;-72) (99;-61) (-87;62)
Sys 79 (71)
(42;53) (73;37)
(00;78) (99;74)
Sys 74
Table 3: Matchings between Sys-excitations and Ref-excitations having for both the NTO-H and NTO-E. The first and second number in a “(.)” is the projection result () for NTO-H and NTO-E, respectively. denotes . Numbers shown in boldface are the excitations plotted in Fig. 12. Given the freedom of orbital shape outside the core region, it should not be surprising to see that the nd and th Ref-excitation match with multiple Sys-excitations.

In Table 3, we see that for both NTO-H and NTO-E of the st-th and th Ref-excitation. As for , of all the involved Sys-excitations discussed above only the th Sys-excitation (hereafter denoted as Sys-) satisfies for both NTO-H and NTO-E. In other words, for the remaining Sys-excitations there are density contributions of NTO-H and/or NTO-E which come from outside of the core region. This would be unsurprising as there remain atoms outside of the core region compared to atoms in the core, and the de-localisation of the density contribution for many Sys-excitations is anticipated. If we lower the threshold to for both NTO-H and NTO-E, the Sys- and also fall into the set. Some excitation properties of the Sys-, , and along with two more excitations, Sys- and , respectively for comparing with the Ref- and Ref-, are listed in Table 4 and the NTO-H(E) orbitals plotted in Fig. 12. It can be seen in Table 4 that the NTO transition vector component, , has generally decreased for the larger-sized Sys-molecule, although the NTO still accounts for at least density contribution.

In Fig. 12 we see that the NTO-H orbitals are essentially the same for the first Ref-excitations. Moreover, the core parts of the Sys-excitations, Sys-, , , and , fully coincide with the core parts of this NTO-H orbital shared by the first Ref-excitations, namely, . The differences between Sys-excitations and Ref-excitations mainly come from the more diffuse NTO-E orbitals.

Strictly speaking, only the Sys- has a significant orbital component of NTO-E coming from the core region, namely . In Table 4, it can be seen that the Sys- has an excitation energy that is close to those of the Ref- and Ref-. The oscillator strength of Sys- also falls in between those of the Ref- and Ref-. As for the orbital shape, the NTO-E of the Sys- in Fig. 12 appears to be a mix between those of the Ref- and Ref-, also reflected by a substantial for both the Ref- and Ref-.

The other Sys-excitations have a high proportion of NTO-H and/or NTO-E on the region outside the core, so comparing these Sys-excitations with the Ref-excitations would be misleading. In fact, the excitation energy can be close, as between the Ref- and Sys-, or significantly different as between the Ref- and Sys-. Nevertheless, the energies of similar-core excitations between the Sys- and Ref-molecule still fall reasonably close to each other, given the wide ranges of Sys- and Ref-excitation energies, i.e. - eV for the Sys-molecule and - eV for the Ref-molecule. Hence, the present method has been able to locate similar excitations of the two molecules among numerous other excitations, by appropriately tuning the thresholds of , , and .

It would be interesting to perform a similar study, and its ring opening reaction, using self-interaction error corrected or reduced functionals such as range-separated ones CH08 (); LCGH11 () in the near futureDHS13 (). For reducing charge-transfer contamination to local excitations, using some technique to confine electronic transitions to a small spatial range to counteract the de-localisation error from the LDA may also be a feasible approach.

Ref-01 Ref-07 Sys-37
E (eV) 6.04 7.46 6.55
f 0.0253 0.0114 0.0070
NTO1 0.9997 0.9998 0.9876
Ref-02 Ref-04 Sys-34 Sys-44
E (eV) 6.52 6.62 6.44 6.75
f 0.0001 0.0199 0.0062 0.0489
NTO1 0.9998 0.9982 0.9130 0.7740
Ref-03 Sys-81
E (eV) 6.57 7.36
f 0.0039 0.0332
NTO1 0.9998 0.8154
Ref-08 Sys-71
E (eV) 7.61 7.20
f 0.0028 0.0054
NTO1 0.9982 0.8851
Table 4: Excitation energy (E), oscillator strength (f), and NTO transition vector component of several excitations listed in Table 3. Sys-excitations and Ref-excitations that have similar core-orbitals are put in the same row. We note that the Ref- and Ref- are substantially similar in their cNTO-E orbitals to have .

V Conclusions and future work

In this work we have developed QNTO analysis based on a non-orthogonal basis, as implemented in the linear-scaling DFT/TDDFT code ONETEP. We use the QNTO method to analyse the transition vectors of electronic excitations of various molecular systems offered by TDDFT calculations.

We first studied the photo-induced ring-opening process of the oxirane molecule. Through a ground state scan along the reactive CCO angle of an oxirane molecule followed by TDDFT calculations (vertical excitations), a profile of several low-lying adiabatic excited state curves on which the excited state driven reaction may be carried can be constructed. However, these adiabatic excited state curves ranked by energies do not give us very much knowledge on how the reaction will proceed. We stress that the reaction process can be made much clearer if the excited state potential energy curves are labelled according to their excitation character instead of their energies. We demonstrated that the transition vectors of these low-lying excitations are dominated by the first NTO electron-hole pair (NTO), and therefore the transition origins of NTO can be used to represent the character of these excitations. Based on this approach, we have successfully located several special geometrical points along the scanned CCO coordinate where the transition origins of two adiabatic excited state curves switch. By reconnecting the adiabatic excited state energy points (calculated by TDDFT) based on transition origins instead of energetic ordering, we are then able to depict how the excited state-driven ring-opening process can proceed. This picture is consistent with the previous excited state molecular dynamics study that found that the population of the (n,) state facilitates the ring-opening process, while the population of (n,) state does not. We have demonstrated that the (n,) state exhibits a monotonically-decreasing energy curve along the widening CCO angle, leading to a ring-opening reaction. For more complicated reactions with several coordinates involved and geometry strongly distorted, more geometries may be used, e.g. those from an excited state scan and molecular dynamics simulation, to help identify the photo-driven reaction pathways of interest. A promising future direction could be to apply this methodology to larger molecules with much more complex reaction pathways.

QNTO analysis has also been employed to study an aryl-substituted oxirane, the trans-2,3-diphenyloxirane molecule. The results demonstrate that the approach can be used to identify specific excitations among a large number of other excitations of a large/complex molecular system. By comparing the excitations of this larger molecule to those of oxirane, we have shown that a small number of mutually similar excitations can be identified and compared between the two molecules. The matching conditions that define similar excitations can be tuned by the magnitudes of: (1) the projection result of the core orbitals, cNTO-H(E), between Sys-excitations and Ref-excitations; (2) the proportion of the cNTO-H(E) accounting for the whole NTO-H(E) in the Sys-excitation; (3) the proportion of the cNTO-H(E) accounting for the whole NTO-H(E) in the Ref-excitation. For example, a Sys-excitation can have a similar cNTO-E to that of a Ref-excitation, but the cNTO-E of the Sys-excitation only accounts for a tiny proportion to the whole NTO-H(E) of Sys-excitation. In such a case, the orbital-section of the NTO-E of Sys-excitation coming from outside of the core region cannot be ignored and the two excitations can not be assigned as similar excitations. The same approach for finding similar excitations through the transition origins would be particularly powerful when there are numerous excitations in a large molecular system. For example, those excitations that are located in a specific geometric region can be readily identified, helping to analyse the chemical properties of that region in a specific energy range of light absorption.

Vi Acknowledgement

The authors would like to thank Louis Lee for very helpful discussions. All calculations in this work were carried out using the Cambridge HPC Service under EPSRC Grant EP/J017639/1. JHL acknowledges the support of Taiwan Cambridge Scholarship. TJZ acknowledges the support of EPSRC Grant EP/J017639/1 and the ARCHER eCSE programme. NDMH acknowledges the support of the Winton Programme for the Physics of Sustainability.


  • (1) P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.
  • (2) W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
  • (3) R. Parr and W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University Press, Oxford, 1994.
  • (4) G. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid, Cambridge University Press, Cambridge, 2005.
  • (5) W. Kohn, A. D. Becke and R. G. Parr, J. Phys. Chem., 1996, 100, 12974.
  • (6) A. D. Becke, J. Chem. Phys., 1997, 107, 8554.
  • (7) J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2007, 128, 084106.
  • (8) E. Runge and E. K. U. Gross, Phys. Rev. Lett., 1984, 52, 997.
  • (9) M. A. L. Marques, C. A. Ullrich, F. Nogueira, A. Rubio, K. Burke and E. K. U. Gross, Time-Dependent Density Functional Theory, Springer-Verlag, Berlin-Heidelberg, 2006.
  • (10) M. E. Casida, in Recent Advances in Density Functional Methods, ed. D. P. Chong, World Scientific, Singapore, 1995, vol. 1.
  • (11) A. Dreuw and M. Head-Gordon, Chem. Rev., 2005, 105, 4009.
  • (12) G. Lever, D. J. Cole, R. Lonsdale, K. E. Ranaghan, D. J. Wales, A. J. Mulholland, C.-K. Skylaris and M. C. Payne, J. Phys. Chem. Lett., 2014, 5, 3614.
  • (13) D. J. Cole, A. W. Chin, N. D. M. Hine, P. D. Haynes and M. C. Payne, J. Phys. Chem. Lett., 2013, 4, 4206.
  • (14) S. J. Fox, C. Pittock, C. S. Tautermann, T. Fox, C. Christ, N. O. J. Malcolm, J. W. Essex and C.-K. Skylaris, J. Phys. Chem. B, 2013, 117, 9478.
  • (15) M. Todorovic, D. R. Bowler, M. J. Gillan and T. Miyazaki, J. R. Soc. Interface, 2013, 10, 20130547.
  • (16) G. Lever, D. J. Cole, N. D. M. Hine, P. D. Haynes and M. C. Payne, J. Phys.: Condens. Matter, 2013, 25, 152101.
  • (17) C. A. White and M. Head-Gordon, J. Chem. Phys., 1994, 101, 6593.
  • (18) V. Weber and M. Challacombe, J. Chem. Phys., 2006, 125, 104110.
  • (19) M. C. Strain, G. E. Scuseria and M. J. G. Frisch, Science, 1996, 271, 51.
  • (20) E. Hernandez and M. J. Gillan, Phys. Rev. B, 1995, 51, 10157.
  • (21) J. L. Fattebert and F. Gygi, Phys. Rev. B, 2006, 73, 115124.
  • (22) J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103.
  • (23) J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón and D. Sánchez-Portal, J. Phys.: Condens. Matter, 2002, 14, 2745.
  • (24) C.-K. Skylaris, P. D. Haynes, A. A. Mostofi and M. C. Payne, J. Chem. Phys., 2005, 122, 084119.
  • (25) R. Baer and M. Head-Gordon, Phys. Rev. Lett., 1997, 79, 3962.
  • (26) S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett., 1999, 82, 2127.
  • (27) L. He and D. Vanderbilt, Phys. Rev. Lett., 2001, 86, 5341.
  • (28) T. J. Zuehlsdorff, N. D. M. Hine, J. S. Spencer, N. M. Harrison, D. J. Riley and P. D. Haynes, J. Chem. Phys., 2013, 139, 064104.
  • (29) J.-H. Li, J.-D. Chai, G. Y. Guo and M. Hayashi, Chem. Phys. Lett., 2011, 514, 362.
  • (30) R. L. Martin, J. Chem. Phys., 2003, 118, 4775.
  • (31) F. Plasser, M. Wormit and A. Dreuw, J. Chem. Phys., 2014, 141, 024106.
  • (32) A. E. Reed, L. A. Curtiss and F. Weinhold, Chem. Rev., 1988, 88, 899.
  • (33) E. D. Glendening, C. R. Landis and F. Weinhold, WIREs Comput. Mol. Sci., 2012, 2, 1.
  • (34) L. P. Lee, D. J. Cole, M. C. Payne and C.-K. Skylaris, J. Comp. Chem., 2013, 34, 429.
  • (35) J.-H. Li, J.-D. Chai, G. Y. Guo and M. Hayashi, Phys. Chem. Chem. Phys., 2012, 14, 9092.
  • (36) C.-K. Skylaris, A. A. Mostofi, P. D. Haynes, O. Diéguez and M. C. Payne, Phys. Rev. B, 2002, 66, 035119.
  • (37) C.-K. Skylaris, P. D. Haynes, A. A. Mostofi and M. C. Payne, J. Phys.: Condens. Matter, 2005, 17, 5757.
  • (38) A. A. Mostofi, P. D. Haynes, C.-K. Skylaris and M. C. Payne, J. Chem. Phys., 2003, 119, 8842.
  • (39) A. Ruiz-Serrano, N. D. M. Hine and C.-K. Skylaris, J. Chem. Phys., 2012, 136, 234101.
  • (40) P. D. Haynes, C.-K. Skylaris, A. A. Mostofi and M. C. Payne, J. Phys.: Condens. Matter, 2008, 20, 294207.
  • (41) P. D. Haynes, C.-K. Skylaris, A. A. Mostofi and M. C. Payne, Chem. Phys. Lett., 2006, 422, 345.
  • (42) N. D. M. Hine, P. D. Haynes, A. A. Mostofi, C.-K. Skylaris and M. C. Payne, Comput. Phys. Comm., 2009, 180, 1041.
  • (43) L. E. Ratcliff, N. D. M. Hine and P. D. Haynes, Phys. Rev. B, 2011, 84, 165131.
  • (44) J. E. Harriman, Phys. Rev. A, 1986, 34, 29.
  • (45) E. D. Glendening, J. K. Badenhoop, A. E. Reed, J. E. Carpenter, J. A. Bohmann, C. M. Morales and F. Weinhold, NBO 5.9 ( the NBO 5.9 Manual, Theoretical Chemistry Institute, University of Wisconsin, Madison, WI, 2011.
  • (46) M. J. Bearpark, M. A. Robb and H. B. Schlegel, Chem. Phys. Lett., 1994, 223, 269.
  • (47) D. R. Yarkony, J. Phys. Chem. A, 2001, 105, 6277.
  • (48) W. Domcke and D. R. Yarkony, Annu. Rev. Phys. Chem., 2012, 63, 325.
  • (49) A. M. Mebel, M. Baer and S. H. Lin, J. Chem. Phys., 2000, 112, 10703.
  • (50) B. G. Levine and T. J. Martinez, Annu. Rev. Phys. Chem., 2007, 58, 613.
  • (51) A. Cembran, F. Bernardi, M. Garavelli, L. Gagliardi and G. Orlandi, J. Am. Chem. Soc., 2004, 126, 3234.
  • (52) B. G. Levine, J. D. Coe and T. J. Martinez, J. Phys. Chem. B, 2008, 112, 405.
  • (53) A. Raab, G. A. Worth, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1999, 110, 936.
  • (54) S. Perun, A. L. Sobolewski and W. Domcke, J. Phys. Chem. A, 2006, 110, 13238.
  • (55) B. K. Kendrick, C. A. Mead and D. G. Truhlar, Chem. Phys., 2002, 277, 31.
  • (56) P. Celani, M. A. Robb, M. Garavelli, F. Bernardi and M. Olivucci, Chem. Phys. Lett., 1995, 243, 1.
  • (57) F. Furche and R. Ahlrichs, J. Chem. Phys., 2002, 117, 7433.
  • (58) E. Tapavicza, I. Tavernelli and U. Rothlisberger, Phys. Rev. Lett., 2007, 98, 023001.
  • (59) N. Govind, M. Petersen, G. Fitzgerald, D. King-Smith and J. Andzelm Comput. Mater. Sci., 2003, 28, 250.
  • (60) B. G. Levine, C. Ko, J. Quenneville and T. J. Martinez, Mol. Phys., 2006, 104, 1039.
  • (61) E. Tapavicza, I. Tavernelli, U. Rothlisberger, C. Filippi and M. E. Casida, J. Chem. Phys., 2008, 129, 124108.
  • (62) F. Cordova, L. J. Doriol, A. Ipatov, M. E. Casida, C. Filippi and A. Vela, J. Chem. Phys., 2007, 127, 164111.
  • (63) C. Hu, O. Sugino and K. Watanabe, J. Chem. Phys., 2014, 140, 054106.
  • (64) M. E. Casida, F. Gutierrez, J. Guan, F.-X. Gadea, D. Salahub and J.-P. Daudey, J. Chem. Phys., 2000, 113, 7062.
  • (65) N. D. M. Hine, J. Dziedzic, P. D. Haynes and C-K. Skylaris, J. Chem. Phys., 2011, 135, 204103.
  • (66) J. P. Perdew and A. Zunger, Phys. Rev. B, 1981, 23, 5048.
  • (67) J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
  • (68) C. Adamo and V. Barone, J. Chem. Phys., 1999,