Expressions for analytical molecular gradients of core-excited states have been derived and implemented for the hierarchy of algebraic diagrammatic construction (ADC) methods up to extended second-order within the core–valence separation (CVS) approximation. We illustrate the use of CVS-ADC gradients by determining relaxed core-excited state potential energy surfaces and optimized geometries for water, formic acid, and benzene. For water, our results show that in the dissociative lowest core-excited state, a linear configuration is preferred. For formic acid, we find that the O K-edge lowest core-excited state is non-planar, a fact that is not captured by the equivalent core approximation where the core-excited atom with its hole is replaced by the “Z + 1” neighboring atom in the periodic table. For benzene, the core-excited state gradients are presented along the Jahn–Teller distorted geometry of the 1s → π* excited state. Our development may pave a new path to studying the dynamics of molecules in their core-excited states.
I. INTRODUCTION
The concept of the potential energy surface (PES) represents a cornerstone in theoretical chemistry, constituting an effective way to visualize and understand the interplay between the electronic and nuclear degrees of freedom that, in turn, underlie the behavior and interactions between molecules. In recent years, an in-depth exploration of the potential energy landscape of valence excited states has become possible thanks to the implementation of analytical gradients at various levels of theory ranging from time-dependent density functional theory (TD-DFT)1,2 to accurate wave function based methods, such as coupled cluster (CC)3,4 and algebraic diagrammatic construction (ADC).4,5 Coupled with the tremendous steps forward made in the field of experimental time-resolved spectroscopy, these theoretical advances have enabled the study of (ultrafast) photo-induced processes and photochemical reactions, such as excited state proton transfer,6 excited state electron transfer in photovoltaic materials,7 photo-dissociation, and photochemical isomerization.8 In particular, the fourth generation radiation sources (x-ray free electron lasers, XFELs) have made time-resolved x-ray spectroscopies on ultra-short timescales available and widespread. The combination between the element- and site-specificity of x-ray techniques with the time-resolved pump–probe setup provides an unprecedented amount of insight into out-of-equilibrium processes on timescales as short as femtoseconds,9,10 characteristic for bond formation and breaking and many other chemical phenomena.
Among the several x-ray spectroscopy techniques, resonant inelastic x-ray scattering (RIXS) stands out as to enable the observation of core-excited state nuclear dynamics of targeted atoms, notably those involved in ultrafast dissociation.11 With the use of RIXS, it is possible to identify the signatures of weak chemical bonds (e.g., hydrogen bonds) and quantify their strengths.12,13 However, the abundance of information encoded into RIXS and time-resolved x-ray spectra is difficult, if not impossible, to decode in the absence of adequate computational tools. It is therefore not surprising that tremendous efforts have been devoted to the development of a plethora of theoretical methods to address core-electron spectroscopies.14 In this context, the hierarchy of ADC methods provides a size-consistent ab initio alternative, attractive due to the fact that it is based on a Hermitian Jacobian matrix and displays high accuracy yet (comparatively) favorable computational scalings. In conjunction with the core–valence separation (CVS) approximation, the second-order extended variant of ADC [or CVS-ADC(2)-x] is on par with the equation-of-motion coupled cluster singles and doubles (EOM-CCSD) method in terms of accuracy and precision.15
It is important to take note that ADC, similar to coupled cluster (CC), is essentially a single-reference method, and its application to systems whose ground state (GS) has multi-reference character is not adequate.16,17 Additionally, ADC is based on a Møller–Plesset (MP) reference state and, therefore, its performance becomes worse for cases where MP breaks down, for example, when the molecular geometry is markedly distorted from the ground state equilibrium, as is the case for dissociation. For an extensive discussion with examples, we refer the reader to Ref. 18. However, we also note that the spin–flip variant of ADC is capable of describing a certain category of systems with multi-reference character and conical intersections,19 and a multi-reference formulation of ADC has been recently presented.16
Compared to CVS-ADC(2)-x, the strict second-order CVS-ADC(2) method provides a reduced accuracy for core-excited states but has the benefit of allowing for the solution of the associated secular equation in the space of single-electron excitations.20 A Fock-matrix driven implementation of this approach without the CVS approximation has been presented in the Gator program21 and is suited for execution in massively parallel high-performance computing (HPC) cluster environments. It is expected that complex and dense biochemical systems with a medium-sized quantum-mechanical (QM) region (described by some 1000 basis functions) will represent routine calculations with this approach. Thus, in the context of core-excited state dynamics, the exploration of core-excited state potential energy surfaces at the CVS-ADC level of theory becomes a realistic and interesting application. This aspect has motivated us to derive and implement molecular gradients up to the CVS-ADC(2)-x level of theory, opening the door for performing core-excited state structure optimizations (as illustrated here), including vibrational effects in x-ray spectrum calculations,22,23 and simulating core-excited state dynamics by using, for example, the surface hopping technique24–26 (potentially on a highly reliable PES with interpolation).27–29 Our derivation is based on the Lagrangian method that is well-established for analytical gradients but has not yet been applied to the CVS approximation.
In Sec. II, we present the derivation of the equations required for CVS-ADC molecular gradients up to the order of CVS-ADC(2)-x with special focus on the amplitude and the orbital responses. Next, in Sec. III, we describe the code implementation, followed by its demonstrative applications to the core-excited states of water, formic acid, and benzene in Sec. IV. Finally, a summary and outlook is provided in Sec. V.
II. THEORY
A. ADC and CVS-ADC
For completeness, let us first briefly describe the ADC method and its CVS approximation. The ADC scheme for the polarization propagator is an electronic structure theory method for describing electronic excited states based on the Møller–Plesset (MP) partitioning of the Hamiltonian. The ADC equations are derived either by expanding the matrix representation of the polarization propagator in a time-dependent perturbation series30–32 or via the intermediate state representation (ISR) by applying the set of excitation operators to the MP ground state.33,34 Here, and are the creation and the annihilation operators, respectively, and the indices i, j, … and a, b, … refer to occupied and unoccupied orbitals in the Hartree–Fock state. The ADC hierarchy is obtained by truncating the perturbation series at the desired order in the fluctuation potential and by applying specific types of excitations (singles, doubles, etc.) to the selected MP ground state. In practice, the excitation energies, the excitation vectors, and the related properties are determined by diagonalizing the ADC matrix, where the matrix size equals the size of the excitation space. It is noted that core excitations are embedded into the larger valence-excitation space as depicted in Fig. 1(a).
(a) Schematic representation of the structure of the full ADC(2) matrix. The numbers in each block represent the highest order of perturbation contributing to the matrix elements. The structure of the ADC(2)-x matrix is the same, with the exception that the elements of the doubles block are treated up to the first order perturbation instead of the zeroth order. (b) Reduced ADC(2) matrix obtained by invoking the CVS approximation. Note that the core, the valence occupied, and the unoccupied orbitals are denoted by {I, J, …}, {i, j, …}, and {a, b, …}, respectively.
(a) Schematic representation of the structure of the full ADC(2) matrix. The numbers in each block represent the highest order of perturbation contributing to the matrix elements. The structure of the ADC(2)-x matrix is the same, with the exception that the elements of the doubles block are treated up to the first order perturbation instead of the zeroth order. (b) Reduced ADC(2) matrix obtained by invoking the CVS approximation. Note that the core, the valence occupied, and the unoccupied orbitals are denoted by {I, J, …}, {i, j, …}, and {a, b, …}, respectively.
To reduce the size of the ADC matrix, one possibility is to pick up only the blocks relevant to the core-excitation process. This is the idea behind the CVS approximation, which is based on the fact that the core and the valence orbitals are well separated spatially and energetically so that the core-excitation and valence-excitation spaces may be considered decoupled.35,36 Several Coulomb integrals involving combinations of core and valence orbitals are therefore neglected, resulting in zeros in the full ADC matrix. This means that a reduced ADC matrix, illustrated in Fig. 1(b), may be used to determine the core-excitation transition energies and amplitudes, as well as core-excited state properties. With this approximation, ADC has been successfully applied to describe the x-ray absorption spectroscopy (XAS) of molecules,37–42 showing a particularly high accuracy and precision at the CVS-ADC(2)-x level.15,42 As they are of direct relevance in the present work, the explicit expressions for the CVS-ADC(2)-x matrix elements are provided in Appendix A.
B. Lagrange formalism for analytical gradients
In the case of non-variational electronic structure methods such as CC and ADC, the derivation of the analytical expressions for energy derivatives can efficiently be performed using the Lagrange formalism.3,43 Assuming that the energy functional has the form , where C and T are non-variational parameters, the total energy derivative with respect to the variable ξ (in our case a nuclear coordinate) is obtained via the chain rule,3,44
Here, C is the molecular orbital (MO) coefficient matrix relating the atomic orbitals (AOs) {ϕμ} to the MOs {χp} via χp = ∑μCμpϕμ, and T contains the remaining non-variational contribution, which, in the case of ADC, is related to the MP perturbation treatment. To avoid computing the derivatives of the non-variational parameters, a Lagrangian L is constructed such that ∂L/∂C (orbital response) and ∂L/∂T (amplitude response) are both zero,2,3,45–48
where Λ and are the sets of undetermined Lagrange multipliers, and fc(C) = 0 and ft(T) = 0 define the constraints for the non-variational parameters C and T. Equations for the undetermined Lagrange multipliers are derived from the stationary conditions, i.e., ∂L/∂C = 0 and ∂L/∂T = 0.45 Once the Lagrange multipliers have been determined, the total derivative of the energy functional with respect to a generic variable ξ is
C. The Lagrange formalism applied to the ADC energy functional
We are interested in computing the derivative of the total energy of an excited state |n⟩ described at the CVS-ADC(2)-x level of theory, with respect to ξ. To do so, we follow closely the derivations for EOM-CC and ADC analytical gradients from Refs. 3 and 5. The energy of the excited state is
where Xn is the ADC eigenvector corresponding to excited state |n⟩ and E0 is the ground state energy with contributions from Hartree–Fock (EHF) and the MP correction (EMP). If we use p, q, … for indexing the entire orbital space, the associated Lagrangian is constructed as follows:
with Λ = {λpq}, Ω = {ωpq}, and as the Lagrange multipliers. Here, F = {Fpq} represents the Fock matrix with the corresponding canonical orbital energies {ϵp}, S = {Spq} is the overlap matrix, and fz(tz) = 0 is the constraint for the non-variational parameters T.3,44 The conditions for the Lagrange multipliers {λpq} and {ωpq} can be obtained from the fact that the Fock matrix is diagonal and, respectively, the fact that the overlap matrix is unity in the case of the underlying Hartree–Fock reference state.3 The t-amplitudes are specific to any given MP order, and in our case, the amplitude response multipliers take the form with the constraint
If we write the total energy and the amplitude contribution in terms of one- and two-particle density matrices (1PDMs γpq and 2PDMs Γpqrs),
we can obtain the gradient only in terms of partial derivatives with respect to ξ,
with γ′ = γ + γA and Γ′ = Γ + ΓA. Here, γA and ΓA denote the amplitude response components of the 1PDM and 2PDM, respectively, and are constructed based on the t-amplitudes and their corresponding Lagrange multipliers.
Expressions for the derivatives required in Eq. (9) are given in Appendix B. The remaining ingredients that we need to determine the nuclear gradient are (i) the 1PDMs and 2PDMs and (ii) the orbital and amplitude response Lagrange multipliers. The density matrices corresponding to the total energy of the excited state are identified by comparing Eqs. (7) and (4), where the latter equation is written explicitly for the desired ADC order. The density matrices corresponding to the amplitude response (γA and ΓA) are derived from Eq. (8), after having found the amplitude response Lagrange multipliers and by replacing the generic condition fz(tz) with the condition specific to the desired MP order, e.g., Eq. (6) for the second order. The derivation of density matrices is illustrated in Appendix C for the singles block of the CVS-ADC(2)/CVS-ADC(2)-x matrix. All the required one- and two-particle density matrices including amplitude response are listed in Table II.
D. Amplitude response Lagrange multipliers
Expressions for the amplitude response Lagrange multipliers are determined from the fact that Ln is stationary with respect to the amplitudes ,
Here, we are momentarily using and to refer to general occupied orbitals spanning over the core {I, J, …} and the valence {i, j, …}.
The MP2 energy is of course given as 36 and is the second-order contribution to the ADC matrix ( Appendix A). The term can be expressed in terms of the t-amplitudes and the zeroth-order one-particle density matrices , as shown in Appendix C.
The explicit forms of the three derivatives required in Eq. (10) are given in Appendix D. Using these derivatives, the final expressions for the amplitude response Lagrange multipliers are obtained as follows:
E. Orbital response Lagrange multipliers
The equations for the orbital response Lagrange multipliers are obtained in a similar manner by formally imposing that the partial derivatives of the Lagrangian with respect to the elements of the orbital transformation matrix are zero,3,5
Practically, working with the half MO-transformed quantities can be avoided by using the following expression:3,5
Quite naturally, calculating the partial derivative of the Lagrangian with respect to Cμt requires the derivatives of the Fock matrix elements, the anti-symmetrized two-electron integrals, and the overlap matrix, and these derivatives are given in Appendix D.
By expressing the Lagrangian in terms of the density matrices [see Eqs. (5) and (9)] and using the derivatives from Appendix D, the partial derivative of the Lagrangian can be written as
resulting in a coupled equation for both the {λ} and the {ω} multipliers,
Here, is non-zero only if t refers to an occupied orbital. In reaching the equation, we have used the symmetry properties in combination with the assumption that we are adopting only real valued orbitals, namely, , ⟨pu‖qt⟩ = ⟨qt‖pu⟩, and Γpqrs = Γqpsr = Γrspq. We have also taken the symmetric representation by imposing λpq = λqp and ωpq = ωqp.
To obtain the equations for the λ orbital response Lagrange multipliers, λ and ω can be decoupled by taking the difference ∑μCμu∂L/∂Cμt − ∑μCμt∂L/∂Cμu.
The system of equations for {λ} is then obtained by choosing u and t from different orbital spaces in the resulting equation,
where the desired ADC variant determines the structure and composition of the density matrices γ′ and Γ′.
While λia are determined from this equation, without frozen-core or frozen-virtual approximations, λij and λab can be safely set to zero in general.3,5 In the case of CVS-ADC, however, additional equations must be solved for the mixed core-occupied block λIj, by setting only the pure blocks λIJ, λij, and λab to zero. As a rule of thumb, any mixed block of λ will be generally non-zero and must be determined using Eq. (18). For example, if a set of virtual orbitals would be frozen, then the mixed block of virtual and frozen-virtual orbitals of λ would also need to be determined. Coming back to our case, the core–virtual (λIa) and the occupied–virtual (λia) blocks are coupled and must be determined together using an iterative technique such as the conjugate gradient algorithm. However, the core–occupied block (λiJ) can be obtained independently. The specific equations used in the case of CVS-ADC(2)-x are listed in Appendix E.
Once λ is computed, the ω Lagrange multipliers can be determined using Eq. (17). Note that all blocks of ω, including the pure blocks such as core–core, occupied–occupied, and virtual–virtual, are non-zero and must be computed. The equations required for the ω multipliers at the CVS-ADC(2)-x level of theory are listed in Appendix F.
As a final note, we comment that the equations presented above are general and can be applied to higher orders of CVS-ADC. The difference between various CVS-ADC orders enters only through the 1PDM and 2PDM (including amplitude response). For x-ray spectroscopy applications, the CVS-ADC(2)-x order represents the current state-of-the art, and we have therefore stopped at this level for the moment. All orders lower than CVS-ADC(2)-x, including the MP2 reference, are immediately available and have been implemented. They are obtained by simply using fewer blocks from the density matrices described above, by ignoring the blocks that are zero at the lower order.
III. IMPLEMENTATION AND COMPUTATIONAL DETAILS
The steps described above have been implemented in Python using the ADC-connect (adcc)49 module, which performs ADC calculations using external electronic structure programs such as PySCF50 and VeloxChem51 for the self-consistent field (SCF) reference state. The resulting Python module named adc_gradient is made available for download under GPLv3 license at https://gitlab.com/iubr/cvs-adc-grad. The gradients have been tested in comparison to numerical gradients determined by the central finite differences algorithm. Comparisons between analytical and numerical gradients at different levels of theory are listed in the supplementary material.
To illustrate the utility of the CVS-ADC gradients, we first determined the optimized geometries of three small molecules, namely, water, formic acid, and benzene. Water was chosen because it has multiple well-separated core-excited states and because its restricted active space (RASPT2) reference potential energy surfaces are available.11 Furthermore, the core-excited states of water have been extensively studied in the literature (see Ref. 52 for a thorough review), exhibiting characteristically different behaviors: e.g., the first core-excited state is dissociative, while the second one is a bound state.11,53 Formic acid was chosen because it provides two different K-edges within the same molecule, namely, O and C edges. In addition to seeing how CVS-ADC treats these two edges, formic acid is an interesting system to study because the vibrational RIXS involving excitations from the hydroxyl O site could be used to detect hydrogen bonding as in the case of acetic acid.13 Finally, benzene was chosen as it constitutes the size limit for the current implementation toward the core-excited state geometry optimization using CVS-ADC(2)-x in combination with a reasonably large basis set (6-311G**). In addition, its first core-excited state geometry has already been reported at the multi-configurational self-consistent field (MCSCF) level of theory.54
The molecular geometries were first optimized at the MP2 theory level with the cc-pVTZ basis set.55 These ground state optimized structures were then used to compute the x-ray absorption spectra at CVS-ADC(2) and CVS-ADC(2)-x in combination with the 6-311++G** basis set.56 The XAS calculations were performed in Q-Chem 5.257 to be able to perform excited-state wave function analysis,58,59 which is currently not available in adcc. The calculated oscillator strengths were broadened using a Gaussian broadening of 0.5 eV full width at half maximum (FWHM), and the peaks were shifted to align with the first experimental peak. The well-separated low energy XAS peaks were selected to represent the states where geometry optimizations and further analyses were to be attempted. In addition to full optimizations, we also performed constrained optimizations to construct contour representations of the relaxed molecular PES. In the case of highly symmetric benzene, we employed effective-core potentials (ECPs) to be able to follow one specific core excitation throughout the optimization procedures. For this purpose, ECPs from the Stuttgart/Cologne group (ECP2MWB)60 were used to replace the C 1s electrons of all C atoms except the one from which the core electron was promoted. The basis set adopted in this case was 6-311G**.56
IV. RESULTS AND DISCUSSION
A. X-ray absorption spectra
Let us first consider the spectral features of the three molecules at the ground state optimized geometries presented in Fig. 2. The CVS-ADC(2) and CVS-ADC(2)-x calculated XAS spectra are shown in Fig. 3 and are compared with the experimental gas-phase data, while the calculated CVS-ADC(2)-x excitation energies and oscillator strengths for the states of interest are listed in Table I. The O K-edge spectrum of water [Fig. 3(a)] is well known and has been extensively studied in the literature. The first XAS peak is assigned to a transition from the O 1s core orbital to the 4a1 anti-bonding molecular orbital (MO). It is also well known that this is a dissociative state with a strong signature in vibrational RIXS.11,12,53,64 The second XAS peak is assigned to a transition to the 2b2 MO and it corresponds to a bound core-excited state.11,64 The higher photon energy features of the spectrum are related to transitions to several close lying orbitals of Rydberg characters.52 Because the first two core-excited states are energetically well separated and have very different characters, these two were selected for constructing relaxed potential energy surfaces by using the CVS-ADC gradients (presented in Sec. IV B).
Ground state MP2 relaxed geometries of (a) water, (b) formic acid, and (c) benzene.
Ground state MP2 relaxed geometries of (a) water, (b) formic acid, and (c) benzene.
XAS of (a) O K-edge of water, (b) O K-edge of formic acid, (c) C K-edge of formic acid, and (d) C K-edge of benzene. CVS-ADC(2)-x, CVS-ADC(2), and the experimental results are shown from bottom to top. The calculated spectra have been broadened with Gaussian functions of 0.5 eV FWHM and horizontally shifted such that the lowest energy peak aligns with the experimental one, with the shift amounts specified on each panel. In the case of benzene, to reduce the computational burden, we limited the calculations to the 40 lowest excited states, covering up to ∼288 eV. With CVS-ADC(2)-x, a spectrum calculated by using ECPs for the five non-core-excited C atoms is also shown (yellow dotted line) by summing up the spectral contributions by all six C atoms. In this case, due to the decreased excitation space size, the calculation could easily cover a much wider spectral range than with the all-electron approach. The experimental spectra are from Ref. 61 (water), Ref. 62 (formic acid), and Ref. 63 (benzene).
XAS of (a) O K-edge of water, (b) O K-edge of formic acid, (c) C K-edge of formic acid, and (d) C K-edge of benzene. CVS-ADC(2)-x, CVS-ADC(2), and the experimental results are shown from bottom to top. The calculated spectra have been broadened with Gaussian functions of 0.5 eV FWHM and horizontally shifted such that the lowest energy peak aligns with the experimental one, with the shift amounts specified on each panel. In the case of benzene, to reduce the computational burden, we limited the calculations to the 40 lowest excited states, covering up to ∼288 eV. With CVS-ADC(2)-x, a spectrum calculated by using ECPs for the five non-core-excited C atoms is also shown (yellow dotted line) by summing up the spectral contributions by all six C atoms. In this case, due to the decreased excitation space size, the calculation could easily cover a much wider spectral range than with the all-electron approach. The experimental spectra are from Ref. 61 (water), Ref. 62 (formic acid), and Ref. 63 (benzene).
CVS-ADC2-(x) calculated excitation energies and oscillator strengths for selected XAS peaks. The excitation energies have been shifted as marked in Fig. 3 to match the experimental first peaks.
. | . | Excitation . | Oscillator . | . |
---|---|---|---|---|
Molecule . | Edge . | energy (eV) . | strength . | Assignment . |
Water | O K-edge | 534.26 | 0.0097 | O 1s → 4a1 |
535.11 | 0.0188 | O 1s → 2b2 | ||
Formic acid | O K-edge | 532.21 | 0.0340 | OC 1s → π* |
535.28 | 0.0020 | OC 1s → σ* | ||
535.79 | 0.0136 | OH 1s → π* | ||
535.98 | 0.0122 | OH 1s → σ* | ||
Formic acid | C K-edge | 288.28 | 0.0673 | C 1s → π* |
Benzene | C K-edge | 285.26 | 0.2365 | C 1s → π* |
. | . | Excitation . | Oscillator . | . |
---|---|---|---|---|
Molecule . | Edge . | energy (eV) . | strength . | Assignment . |
Water | O K-edge | 534.26 | 0.0097 | O 1s → 4a1 |
535.11 | 0.0188 | O 1s → 2b2 | ||
Formic acid | O K-edge | 532.21 | 0.0340 | OC 1s → π* |
535.28 | 0.0020 | OC 1s → σ* | ||
535.79 | 0.0136 | OH 1s → π* | ||
535.98 | 0.0122 | OH 1s → σ* | ||
Formic acid | C K-edge | 288.28 | 0.0673 | C 1s → π* |
Benzene | C K-edge | 285.26 | 0.2365 | C 1s → π* |
Figures 3(b) and 3(c) show the O and the C K-edge spectra of gas-phase formic acid, respectively. The first peak in the O K-edge spectrum corresponds to a transition from the 1s orbital of the carbonyl oxygen [denoted OC in Fig. 2(b)] to the π*-orbital, resulting from the hybridization between the out-of-plane C and O 2p orbitals. The second O K-edge peak has several contributions, and the most important ones are by the transitions from OH 1s to π* and σ*. The C K-edge spectrum has a strong and energetically well separated first peak corresponding to a transition from C 1s to π*, followed by a series of weaker transitions at higher energies. Given the unambiguous assignment and separation in energy, the first core-excited states arising in the O and the C K-edge XAS spectra have been selected for geometry optimizations. The CVS-ADC(2)-x relaxed geometries will be compared to those obtained with an equivalent core approximation in Sec. IV C.
In the case of the C K-edge spectrum of benzene shown in Fig. 3(d), we can observe a strong π*-peak, corresponding to transitions from the nearly degenerate core orbitals to the doubly degenerate π* MOs. This is followed by low intensity Rydberg transitions at higher photon energies.65 Due to symmetry, the vertical transitions that contribute to the first XAS peak correspond to six nearly degenerate core excitations, out of which only one has non-zero oscillator strength. Although the ground state (GS) PES leads to an equilibrium molecular configuration with a high symmetry, the core-excited PES has six potential wells with slight distortions corresponding to Jahn–Teller minima. The Jahn–Teller effect will be discussed in relation to the core-excited state geometry optimization in Sec. IV D.
B. Relaxed core-excited state potential energy surfaces
The two-dimensional (2D) potential energy surfaces of the ground state (GS) and the two core-excited states of water are shown in Fig. 4. Figure 4(a) shows the unrelaxed PES, determined at a fixed bond angle θ = 103.7° (the MP2 equilibrium value), obtained by scanning the two OH bond lengths, R1 and R2. The GS PES is obtained at the MP2 level of theory, while the core-excited state PES is calculated using CVS-ADC(2)-x. The relaxed potential energy surfaces were obtained by optimizing the bond angle for each combination of (R1 and R2) and are shown in Fig. 4(b). For comparison, RASPT2 data from Ref. 11 are also shown in Fig. 4(c), which were determined at a fixed bond angle of θ = 104.2° (equilibrium value from Ref. 11). Comparing the MP2 and the CVS-ADC(2)-x results to the reference RASPT2 data reveals that the natures of both the bound and the dissociative states are well captured. In the case of the ground state, the RASPT2 and the MP2 PES in the vicinity of the equilibrium point are, for all intents and purposes, identical. As the bonds become elongated, the MP2 energy increases slightly more compared to RASPT2. Comparing the relaxed vs the unrelaxed PES, changes in energy by the geometry relaxation are rather small.
Potential energy surfaces of the ground state (GS), O1s → 4a1, and O1s → 2b2 core-excited states of water. (a) Unrelaxed PES determined at a fixed bond angle (103.7°), (b) relaxed PES after constrained optimization of the bond angle, and (c) reference RASPT2 unrelaxed PES from Ref. 11. For GS and O1s → 2b2, the zero energy references are the energy minima in each panel, marked with E0 and Em. For O1s → 4a1, they are the saddle points, marked with Es.
Potential energy surfaces of the ground state (GS), O1s → 4a1, and O1s → 2b2 core-excited states of water. (a) Unrelaxed PES determined at a fixed bond angle (103.7°), (b) relaxed PES after constrained optimization of the bond angle, and (c) reference RASPT2 unrelaxed PES from Ref. 11. For GS and O1s → 2b2, the zero energy references are the energy minima in each panel, marked with E0 and Em. For O1s → 4a1, they are the saddle points, marked with Es.
When it comes to the O1s → 4a1 core-excited state, CVS-ADC(2)-x captures the dissociative character of the state relatively well. From the saddle point along the symmetric stretching direction, the energy increases somewhat more steeply than with RASPT2, while along the asymmetric stretching direction, the energy decreases less steeply. The energy on the relaxed PES changes slightly more abruptly than on the unrelaxed PES, but the differences are, in general, again quite small.
The O1s → 2b2 core-excited state is a bound state with an energy minimum with slightly elongated OH bonds compared to the GS one (ROH =∼1.08 Å). Also in this case, from the minimum along symmetric stretching, the CVS-ADC(2)-x energy increases more steeply compared to RASPT2. On the other hand, the situation is reversed along the asymmetric stretching coordinate, where RASPT2 increases slightly faster. Comparing the relaxed vs the unrelaxed CVS-ADC(2)-x PES, on the relaxed PES, the energy increase along the symmetric stretching coordinate is less steep, while along the asymmetric stretching coordinate, the relaxed and unrelaxed PES remain virtually identical. Although no dramatic changes are observed when comparing the relaxed and unrelaxed 2D PES contours, the changes in the optimized bond angles can be quite large, especially in the case of the core-excited states with potentially flat PES morphology. Figure 5(a) shows 2D representations of the optimized θ as a function of R1 and R2 for the three states of water. Indeed, dramatic variations in the bond angle are observed for the dissociative O1s → 4a1 core-excited state. For many combinations of (R1 and R2) in a wide region, the molecule prefers to become linear. In contrast, when both bonds are short, the relaxed bond angle lies between 100° and 120°. Interestingly, when both bonds are elongated, the angle changes its preference toward 120°–165°. We should note that some of the bending mode potential energy curves for this core-excited state are quite flat (Fig. S1 in the supplementary material) and sometimes geometric convergence is difficult to achieve.
(a) Two-dimensional representation of the optimized bending angle θ as a function of the OH bond lengths (R1 and R2). (b) Bending potential energy curves (PECs) determined at fixed OH bond lengths (equilibrium value, 0.96 Å), in comparison to RASPT2 data from Ref. 11. Closely matching harmonic PECs are also drawn for visual comparisons.
(a) Two-dimensional representation of the optimized bending angle θ as a function of the OH bond lengths (R1 and R2). (b) Bending potential energy curves (PECs) determined at fixed OH bond lengths (equilibrium value, 0.96 Å), in comparison to RASPT2 data from Ref. 11. Closely matching harmonic PECs are also drawn for visual comparisons.
In the case of the O1s → 2b2 core-excited state, the relaxed θ value at the energy minimum is very close to the GS equilibrium value, but the energy minimum position is shifted toward slight bond elongations, as mentioned before. Along the symmetric stretching coordinate, the value of the relaxed bond angle decreases until it reaches a minimum of ∼88° at around RH = 1.5–1.6 Å, after which it increases again. Along the asymmetric stretching coordinate, the θ isovalue lines are rather vertical and parallel to each other. This is different than in the GS, where the relaxed θ decreases both along symmetric and asymmetric stretching when moving from the equilibrium toward bond elongation.
Finally, Fig. 5(b) depicts the unrelaxed potential energy curves computed using MP2/CVS-ADC(2)-x in comparison to RASPT2 data, represented as functions of θ − θ0 with the ground state equilibrium value θ0. The MP2/CVS-ADC(2)-x results resemble the RASPT2 ones rather well, with slight differences observed with the core-excited states. For these, the CVS-ADC(2)-x curves are shifted downward by ∼1.2 eV and have shallower minima compared to RASPT2.
C. Core-excited state geometry optimization
For the lowest C and O core-excited states of formic acid, we performed full geometry optimizations at the CVS-ADC(2)-x level of theory. The results for the two states are presented in Figs. 6(a) and 6(b), respectively. Compared to the ground state geometry shown in Fig. 6(e), the C1s → π* core-excited state one has more elongated CO bonds (by 0.05 Å) and a shortened CH bond. The O–C–O bond angle is a few degrees smaller, but the dihedral angles do no change and the molecule remains planar. If we compare it with the optimized geometry bearing an equivalent core as in Fig. 6(c) at the MP2 level, the geometry is indeed almost identical. This represents the equivalent core approximation, where the core-excited atom and its localized hole are replaced by the atom’s “Z + 1” nearest neighbor in the periodic table. In the case of C, which gets replaced by N, this approximation is known to work quite well and has provided good agreement to experimental data even for large molecules, such as C60.66 In the case of O, in contrast, a replacement by F with the equivalent core approximation is known to fail at times due to the high electronegativity of fluorine that sometimes induces unexpected charge transfer.67 In the present oxygen core excitation, the approximation by HCFOH indeed displays rather large differences [Fig. 6(d)].
CVS-ADC(2)-x optimized geometries of formic acid (HCOOH) for the (a) C 1s → π* and (b) the OC 1s → π* states in comparison with the geometries based on the equivalent core approximation with (c) HNOOH and (d) HCFOH. (e) GS geometry of HCOOH and (f) CVS-ADC(2)-x attachment density for the OC 1s → π* core excitation at the relaxed core-excited state geometry.
CVS-ADC(2)-x optimized geometries of formic acid (HCOOH) for the (a) C 1s → π* and (b) the OC 1s → π* states in comparison with the geometries based on the equivalent core approximation with (c) HNOOH and (d) HCFOH. (e) GS geometry of HCOOH and (f) CVS-ADC(2)-x attachment density for the OC 1s → π* core excitation at the relaxed core-excited state geometry.
The OC 1s → π* transition we are considering involves the 1s electron from the carbonyl oxygen [OC, Fig. 2(b)]. In the CVS-ADC(2)-x relaxed geometry, the C=O bond is elongated to 1.38 Å, whereas the C–O bond length decreases by a small amount. In the equivalent core approximation, the C=O bond also elongates but to a lesser degree, while the C–O bond remains almost the same. The O–C–O bond angle decreases in both cases but to a larger extent in CVS-ADC(2)-x. The most important difference, however, is the fact that the CVS-ADC(2)-x optimized structure is non-planar, with an O–C–O–H dihedral angle of ∼11° [see the side-view with the attachment density in Fig. 6(f)]. Note that the attachment density in Fig. 6(f), which represents the charge rearrangement during core excitation, is non-symmetric in the z-direction (i.e., perpendicular to the molecular plane). This aspect is, in fact, correlated with the non-planarity of the molecule in the OC 1s → π* core-excited state. In contrast, within the equivalent core approximation, the molecule remains basically planar with an almost negligible out-of-plane distortion of ∼0.1°.
D. The first core-excited state of benzene and the Jahn–Teller effect
Finally, we discuss the first core-excited state of benzene in connection to the Jahn–Teller effect. As already mentioned in Sec. IV A, due to the high symmetry, the first core-excited state of benzene is practically sixfold degenerate. Figure 7(a) shows the MOs and their diagrams connected to the first XAS peak. In the D6h point group, the C 1s MOs split into four types, two non-degenerate and two doubly degenerate, as marked by σ1 to in the figure. Of course, the energy differences between these are very small. A similarly patterned splitting takes place for the valence orbitals, but the energy differences are much larger (up to a few eVs). Important here are the doubly degenerate and orbitals (e2u) with anti-bonding characters, which constitute the lowest unoccupied molecular orbitals (LUMOs).
(a) Molecular orbitals and their diagrams involved in the C 1s → π* transition of benzene with D6h symmetry. (b) Potential energy curves of the first and the second group of core-excited states, computed using ECPs (yellow) or the full set of C 1s MOs (blue, gray) as a function of a geometric distortion x, where x represents the interpolation weight of the relaxed core-excited geometry. The inset shows the XAS spectrum of benzene determined at the ground state equilibrium geometry with CVS-ADC(2)-x, shifted by 0.33 eV to lower photon energies as in Fig. 3. The peak denoted 1 in the inset corresponds to the first group of excited states, while the dotted line denoted 2 marks the energy position of the second group of excited states that are silent in XAS due to symmetry (see text). In the presence of a core hole, Jahn–Teller symmetry breaking takes place and the geometry relaxes to one potential energy well: C 1s (x = 1.0) or C 1s (x = −1.0). The associated bond length alternation is depicted with a color bar. The norm of the gradient is denoted by ∇ and given (in atomic units) next to the corresponding molecular structure.
(a) Molecular orbitals and their diagrams involved in the C 1s → π* transition of benzene with D6h symmetry. (b) Potential energy curves of the first and the second group of core-excited states, computed using ECPs (yellow) or the full set of C 1s MOs (blue, gray) as a function of a geometric distortion x, where x represents the interpolation weight of the relaxed core-excited geometry. The inset shows the XAS spectrum of benzene determined at the ground state equilibrium geometry with CVS-ADC(2)-x, shifted by 0.33 eV to lower photon energies as in Fig. 3. The peak denoted 1 in the inset corresponds to the first group of excited states, while the dotted line denoted 2 marks the energy position of the second group of excited states that are silent in XAS due to symmetry (see text). In the presence of a core hole, Jahn–Teller symmetry breaking takes place and the geometry relaxes to one potential energy well: C 1s (x = 1.0) or C 1s (x = −1.0). The associated bond length alternation is depicted with a color bar. The norm of the gradient is denoted by ∇ and given (in atomic units) next to the corresponding molecular structure.
At the GS equilibrium geometry, transitions from the six core orbitals to the doubly degenerate LUMO form two groups of sixfold degenerate core-excited states. The group denoted 1 in Fig. 7(b) corresponds to the first peak in the XAS spectrum [peak 1 in the inset of Fig. 7(b)], while the group denoted 2, ∼0.8 eV higher in energy, is silent in the XAS spectrum. Due to the orbital symmetries, only one transition out of these has a non-zero oscillator strength. This transition, responsible for the first XAS peak, is from () to () (or from e2g to e2u) and resides in group one.65
Due to the degeneracy of these transitions, benzene represents a difficult case for the core-excited state geometry optimization. As soon as the molecular structure is distorted, the degeneracy is lifted and the energy ordering of the core-excited states changes with small atom displacements. Consequently, following one particular state throughout the optimization procedure becomes very tricky. To circumvent this issue, we used ECPs to represent the five non-core-excited 1s electrons, meaning basically that we enforced the core–hole to localize on a particular atom. Physically, this is equivalent to an actual x-ray absorption situation, where the x-ray photon creates a localized core hole on a particular site. In this picture, the XAS intensity is a sum of equal contributions from the six equivalent atoms. At the CVS-ADC(2)-x level of theory, this ECP based approach produces essentially identical XAS spectra as the all-electron approach [Fig. 3(d)].
Using ECPs, we optimized the geometry of the first core-excited state of benzene at the CVS-ADC(2)-x level of theory. By interpolating between the core-excited geometry and the ground state geometry, we constructed the potential energy curves corresponding to the first (1s → π* peak in XAS) and second (not present in XAS due to zero oscillator strength) group of core excitations. These curves are shown in Fig. 7(b), alongside the ground state geometry, the core-excited state geometry, and several intermediate points. Energies from both the full-electron and the ECP calculations are shown for comparison. The potential energy curves are represented as a function of the interpolation weight x, where the intermediate geometry corresponding to x (Gx) has been constructed as a linear combination between the ground state (GGS) and, respectively, the ECP core-excited state (G*) relaxed geometries,
As mentioned already, in the delocalized MO picture at the high symmetry configuration, each of the two groups of core-excited states has a sixfold degeneracy. Note also that these nearly degenerate transitions are from core orbitals to the degenerate LUMO, consisting of and . More specifically, the only transition with non-zero oscillator strength in the first group of excited states bears and characters. In the second group, the transition characters switch, i.e., and , resulting in zero oscillator strength due to symmetry. This is the reason why no peak corresponding to the second group of excitations is present in the XAS spectrum at the photon energy position marked by the dotted line in the inset of Fig. 7(b).
Once the geometry distorts there are two effects that take place. (i) The energy differences between the nearly degenerate core MOs increase, so the near degeneracies are lifted. However, the energy difference between the two doubly degenerate couples (σ2, σ3) and () remains small enough even at large distortions (x = 1 in the figure), so these couples still remain nearly degenerate. (ii) The degeneracy of the LUMOs is lifted and one of the two π* orbitals is stabilized over the other, depending on the direction of the distortion. To exemplify, two of the six potential energy wells are illustrated in Fig. 7, one where a transition to is lowest in energy and the other where a transition to is lower. Notice that the excited states arising from fully non-degenerate core orbitals are marked with full lines in Fig. 7, while those arising from the two pairs of nearly degenerate core orbitals are marked with dotted lines. In the ECP approach, only one of the six curves is captured, but the others can be obtained by rotating the molecular frame by multiples of 60°. While the ECP curve matches the lowest energy curve in the first group of excited states in the delocalized full-electron MO picture, in the second group the ECP curve follows one of the higher energy states. Actually, this still is the lowest energy state if we consider only the states arising from non-degenerate core orbitals. The other two states that are lower in energy in the full-electron picture arise from the degenerate pairs, which cannot be captured in the ECP picture.
Finally, we note that the way that the CVS-ADC(2)-x core-excited state optimized geometry changes from the GS geometry is in agreement with the results obtained using the MCSCF method, as well as the equivalent core approximation in Ref. 54. Compared to the GS equilibrium structure, the C–H bond at the core-excited site contracts, while the C–C bonds at the same site elongate. The other C–C bonds contract and elongate in an alternating manner, as illustrated in Fig. 7(b).
V. SUMMARY AND OUTLOOK
In summary, we have presented a derivation and implementation of analytical nuclear gradients up to the level of CVS-ADC(2)-x in the hierarchy of the ADC methods. The resulting Python module named adc_gradient is integrated with ADC connect (adcc)49 and is made available for download under GPLv3 license at https://gitlab.com/iubr/cvs-adc-grad.
Sample illustrations of this development were provided for the lowest core-excited states of water, formic acid, and benzene, where we have determined optimized core-excited state geometries and relaxed potential energy surfaces. In the case of water, we observed that in the dissociative lowest core-excited state, the molecular geometry relaxes to a linear configuration. In the case of formic acid, the OC 1s → π* core-excited state relaxes into a non-planar geometry in contrast to the equivalent core approximation with fluorine that does not capture this effect. Finally, in the case of benzene, the core-excited state gradients are determined along the Jahn–Teller distorted geometry of the C 1s → π* excited state. At the relaxed geometry in this state, the C–H bond at the excitation site contracts, while the C–C bonds elongate and contract in an alternating manner. This reduces the symmetry from D6h (ground state) to C2v (1s → π*).
The core-excited optimized geometries and relaxed potential energy surfaces are important for the inclusion of vibrational effects in computing x-ray spectra as well as for the simulation of vibrational RIXS. Furthermore, the molecular gradients can be used to perform core-excited state dynamics using, e.g., surface hopping. Toward such utilizations, our gradient module will become embedded into future releases of the Gator program21 in combination with the OpenMP/MPI-parallel Fock-matrix driven implementation of ADC(2). This will open up new doors for the calculation of analytical molecular gradients of core-excited states for complex molecular systems comprising medium-sized cores treated at the level of quantum mechanics (employing up to about 1000 basis functions) in combination with a large-scale environment treated at the level of a classical polarizable embedding.
SUPPLEMENTARY MATERIAL
See the supplementary material for comparisons between analytical and numerical gradients at different orders of CVS-ADC for the C, N, and O K-edges of selected molecules. Additionally, we provide several bending mode potential energy curves for the O1s → 4a1 dissociative core-excited state of the water molecule.
ACKNOWLEDGMENTS
Financial support from the Swedish Research Council (Grant Nos. 2017-06419 and 2018-4343) and from the National Research Foundation (NRF) of Korea (Grant No. 2020R1A5A1019141) funded by the Ministry of Science and ICT (MSIT) is acknowledged. Computational resources were provided by the Swedish National Infrastructure for Computing (SNIC). We also express our gratitude to Maximilian Scheurer for providing access and guidance to adcc, Vinicius Vaz Da Cruz for providing RASPT2 data alongside insightful explanations, Viktoriia Savchenko for helpful discussions on the topic of RIXS, and Manuel Hodecker for help in finding a few errors in the equations.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX A: CVS-ADC(2)-x MATRIX ELEMENTS
The explicit forms of the CVS-ADC(2)-x matrix elements are given below. The upper case indices refer to core orbitals (I, J, …) and the lower case ones refer to valence occupied (i, j, …) or virtual (a, b, …) orbitals. The orbital eigenvalues are denoted by ϵ, while δ is the Kronecker delta. The elements that are set to zero with the CVS approximation are underlined (A3) or crossed out (all other equations).
The (p-h; p-h) block elements are36
with
The (p-h; 2p-2h) and the (2p-2h, p-h) matrix elements are36
Finally, the (2p-2h; 2p-2h) block elements are36
APPENDIX B: INTEGRAL DERIVATIVES
As usual, the following partial derivatives are required for the gradient:3
Here, hpq represents a core Hamiltonian matrix element and the superscript indicates a partial derivative with respect to variable ξ. These derivatives are composed of AO-based derivative integrals,
APPENDIX C: DENSITY MATRICES
We illustrate below how to find γ and Γ for the singles block of the CVS-ADC(2)-x matrix. We first need to re-write the energies on the right hand side of Eq. (4) in terms of the Fock matrix elements and the two-electron integrals. This is done by writing explicit expressions for the term using the CVS-ADC(2)-x matrix elements of Eqs. (A1), (A2), and (A4). At the zeroth order,
resulting in the following 1PDM terms:
The first order ADC matrix gives rise to 2PDM components related to the core orbitals,
leading to
Finally, the second-order ADC matrix results in 2PDM components involving only valence indices,
where tijab are MP(2) t-amplitudes [Eq. (6)]. From this equation, we can identify the 2PDM elements as
All the other density matrices, including those corresponding to the ground state energy and the amplitude response, are obtained in a similar way. The final expressions for all density matrices required for gradients up to CVS-ADC(2)-x order are listed in Table II.
Final 1PDM and 2PDM expressions required at the CVS-ADC(2)-x level of theory. Only the non-zero matrix elements are listed. The Einstein summation convention is used for brevity.
Density matrix . | Reference state terms . | ADC matrix terms . | t-response terms . |
---|---|---|---|
δIJ | |||
δij | |||
−2δIKδJL | |||
−δikδJL | |||
−2δikδjl | |||
a | |||
a | |||
Density matrix . | Reference state terms . | ADC matrix terms . | t-response terms . |
---|---|---|---|
δIJ | |||
δij | |||
−2δIKδJL | |||
−δikδJL | |||
−2δikδjl | |||
a | |||
a | |||
A factor of has been added to the expressions for ΓkJIa and Γiabc to account for a scaling factor used in the coupling blocks in the actual implementation of CVS-ADC(2)-x.36
APPENDIX D: DERIVATIVES REQUIRED FOR AMPLITUDE AND ORBITAL RESPONSES
1. Amplitude response
The derivatives required for amplitude response are
2. Orbital response
The derivatives required for orbital response are
APPENDIX E: EQUATIONS FOR THE λ ORBITAL RESPONSE LAGRANGE MULTIPLIERS
For CVS-ADC, the only non-zero blocks of λ are the mixed blocks: occupied–core, core–virtual, and occupied–virtual. The occupied–core block can be computed using Eq. (18) by setting u = i and t = J,
Because the first sum vanishes owing to the symmetries of γ′ and λ, we get
Once λiJ are obtained, the coupled blocks λia and λIa are determined together through an iterative technique, such as the conjugate gradient algorithm. The equations required are
or equivalently,
for the occupied–virtual block.
Similarly,
for the core–virtual block.
APPENDIX F: EQUATIONS FOR THE ω ORBITAL RESPONSE LAGRANGE MULTIPLIERS
All blocks of the ω orbital response Lagrange multipliers are generally non-zero and must be computed. For the CVS approximation, the blocks are core–core, occupied–core, occupied–occupied, core–virtual, occupied–virtual, and virtual–virtual. From Eq. (17), the equations required for each block are as follows: