Analytical Gradients for Core-Excited States in the Algebraic Diagrammatic Construction (ADC) Framework

Here we present a derivation of the analytical expressions required to determine nuclear gradients for core-excited states at the core-valence separated algebraic diagrammatic construction (CVS-ADC) theory level. Analytical gradients up to and including the extended CVS-ADC(2)-x order have been derived and implemented into a Python module, adc_gradient. The gradients were used to determine core-excited state optimized geometries and relaxed potential energy surfaces for the water, formic acid, and benzne molecules. 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 determin-ing 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 conﬁguration is preferred. For formic acid, we ﬁnd 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” neighbouring 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. The supplementary material contains comparisons between analytical gradients and numerical gradients computed using the central ﬁnite diﬀerence algorithm (section I). We list examples for the C, N, and O K-edge, respectively. The molecular geometries are from the supplementary material of Ref. [1]. The ﬁrst C1s/N1s/O1s core-excited state of each molecule is given, where 10 singlet states were included in the ADC matrix diagonalization. The basis set used was def2-SVP. 2 The convergence tolerance in the ADC calculation was 10 − 6 and the displacement step in the ﬁnite diﬀerence algorithm was 10 − 4 bohr. Additionally, we show several bending mode potential energy curves for the dissociative O 1s → π ∗ excited state of the water molecule, computed at diﬀerent combinations of OH bond lengths (section II).


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 behaviour 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 4 th 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 pumpprobe setup provides an unprecedented amount of insight into out-of-equilibrium processes on timescales as short as the 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 coreelectron spectroscopies. 14 In this context, the hierarchy of ADC methods provides a sizeconsistent ab initio alternative, attractive due to the facts 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 doubles (EOM-CCSD) method in terms of accuracy and precision. 15 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. 16 A Fock-matrix driven implementation of this approach without the CVS approximation has been presented in the Gator program 17 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 1,000 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 spectra calculations, 18,19 and simulating core-excited state dynamics by using, for example, the surface hopping technique [20][21][22] (potentially on a highly reliable PES with interpolation). [23][24][25] Our derivation is based on the Lagrangian method which is well-established for analytical gradients, but has not yet been applied to the CVS approximation.
In the next section, 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 Section III, we describe the code implementation, followed by its demonstrative applications to the core-excited states of water, formic acid, and benzene in Section IV. Finally, a summary and outlook is provided in the final section.

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) partition-ing of the Hamiltonian. The ADC equations are derived either by expanding the matrix representation of the polarization propagator in a time-dependent perturbation series, [26][27][28] or via the intermediate state representation (ISR) by applying the set of excitation opera- ...} to the MP ground state. 29,30 Here,â † andâ are respectively the creation and the annihilation operators, 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. 1a. 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. 31,32 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. 1b, 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 employed to determine X-ray absorption spectra of a variety of molecules, [33][34][35][36][37][38] showing a particularly high accuracy and precision at the CVS-ADC(2)-x level. 15,38 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,39 Assuming that the energy functional has the form E = E(C, T, ξ), 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,40 Here, C is the molecular orbital (MO) coefficient matrix relating the atomic orbitals {φ µ } 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: 3,41-45 where Λ andT are the sets of undetermined Lagrange multipliers, and f c (C) = 0, f t (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. 41 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 X n is the ADC eigenvector corresponding to excited state |n and E 0 is the ground state energy with contributions from Hartree-Fock (E HF ) and the MP correction (E MP ). If we use p, q, ... for indexing the entire orbital space, the associated Lagrangian is constructed as follows: is the overlap matrix, and f z (t z ) = 0 is the constraint for the non-variational parameters T. 3,40 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 formt z =t ijab 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 .
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 Eq. (7) and Eq. (4), where the latter equation is written explicitly for the desired ADC order. The density matrices corresponding to amplitude response are derived from Eq. (8), after having found the amplitude response Lagrange multipliers. 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 I of the same Appendix.

D. Amplitude response Lagrange multipliers
Expressions for the amplitude response Lagrange multipliers are determined from the fact that L n is stationary with respect to the amplitudes tīj ab : ph,ph X ∂tīj ab + r,s,c,dtkl cd ∂f z (tkl cd ) ∂tīj ab = 0 .
Here, we are momentarily usingī andj 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 E MP2 = −1/4 īj ab tīj ab īj ||ab 32 and M (2) ph,ph is the second-order contribution to the ADC matrix (Appendix A). The term X † M (2) ph,ph X can be expressed in terms of the t-amplitudes and the zeroth-order one-particle density matrices ab , 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:

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, δ t o 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 γ pq = γ qp , 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 coreoccupied 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 in the 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) 46  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. [49] for a thorough review), exhibiting characteristically different behaviours: e.g., the first coreexcited state is dissociative while the second one is a bound state. 11,50 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 states, its vibrational RIXS is also an interesting target to compare, with which the 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**). Also, its first core-excited state geometry has already been reported at the multi-configurational self-consistent field (MCSCF) level of theory. 51 The molecular geometries were first optimized at the MP2 theory level with the cc-pVTZ basis set. 52 These ground state optimized structures were then used to compute the X-ray absorption spectra (XAS) at CVS-ADC(2) and CVS-ADC(2)-x in combination with the 6-311++G** basis set. 53 The XAS calculations were performed in Q-Chem 5.2 54 to be able to perform excited-state wave function analysis, 55,56 which currently is 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 experi- mental 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) 57 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**. 53

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.
The O K-edge spectrum of water (Fig. 3a) 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 4a 1 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,50,61 The second XAS peak is assigned to a transition to the 2b 2 MO and it corresponds to a bound core-excited state. 11,61 The higher photon energy features of the spectrum are related to transitions to several close lying orbitals of Rydberg characters. 49 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). 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. 3d, 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. 62 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 4a shows the unrelaxed PES, When it comes to the O1s → 4a 1 core-excited state, CVS-ADC(2)-x captures the dis-sociative 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.

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 respectively presented in Figs. 6a and 6b. Compared to the ground state geometry shown in Fig. 6e, 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. 6c 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 neighbour 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 C 60 . 63 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. 64 In the present oxygen core excitation, the approximation by HCFOH indeed displays rather large differences (Fig. 6d).
The O C 1s → π * transition we are considering involves the 1s electron from the carbonyl oxygen (O C , Fig. 2b). 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 symmetry, the first core-excited state of benzene is practically sixfold degenerate. Figure 7a shows the MOs and their diagrams connected to the first XAS peak. In the D 6h point group, the C 1s MOs split into four types, two non-degenerate and two doubly-degenerate, as marked by σ 1 to σ * 6 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 π * 4 and π * 5 orbitals (e 2u ) with anti-bonding characters, which constitute the lowest unoccupied molecular orbitals (LUMOs). At the GS equilibrium geometry, the vertical transitions which give rise to the first XAS peak involve the six core orbitals and the doubly degenerate LUMO.
However, due to the orbital symmetries, only one transition out of these six has a non-zero oscillator strength. The transition is from (σ * 4 ,σ * 5 ) to (π * 5 ,π * 4 ) (or from e 2g to e 2u ). 62 In the presence of a core hole, Jahn-Teller symmetry breaking takes place and the geometry relaxes to one potential energy well: C 1s→ π * 5 (x = 1.0) or C 1s→ π * 4 (x = −1.0). The associated bond length alternation is depicted with a color bar. The norm of the gradient is denoted by ∇ and given next to the corresponding molecular structure.
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. 3d).
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. 7b, 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.
In the delocalized MO picture, at the high symmetry configuration, each of the two groups of core-excited states has a six-fold degeneracy. Note also that these nearly degenerate transitions are from core orbitals to the degenerate LUMO, consisting of π * 4 and π * 5 . 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 (σ * 4 ,σ * 5 ) remains small enough even at large distortions (x = 1 in the figure), so these couples still remain nearly degenerate. (ii) The degeneracy of the LUMO 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 π * 5 is lowest in energy, and another where a transition to π * 4 is lower. 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 • .
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. 51. Compared to the GS equilibrium structure, the C-H bond in the core-excited site contracts, while the C-C bond at the same site elongates. The other C-C bonds contract and elongate in an alternating manner, as illustrated in Fig. 7b.

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) 46 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 O C 1s→ π * core-excited state relaxes into a non-planar geometry, in contrast to the equivalent core approximation with fluorine which 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 D 6h (ground state) to C 2v (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. Towards such utilizations, our gradient module will become embedded into future releases of the Gator program, 17

SUPPLEMENTARY MATERIAL
The supplementary material contains comparisons between analytical and numerical gradients at different orders of CVS-ADC for the C, N, and respectively O K-edges of selected molecules. Additionally, we provide several bending mode potential energy curves for the O1s→4a 1 dissociative core-excited state of the water molecule.

Appendix B: Integral derivatives
As usual, the following partial derivatives are required for the gradient: 3 Here, h pq 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: (B7)

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 X † n MX n term using the CVS-ADC(2)-x matrix elements of Eqs. (A1), (A2), and (A4). At the zeroth order: x * Jb a δ ab δ IJ x Ia − I,a,J,b 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: ph,ph X = I,a,J,b where t ijab 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 I.
[a] A factor of 1/ √ 2 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. 32 Appendix D: Derivatives required for amplitude and orbital responses 1. Amplitude response µ C µu ∂ pq||rs ∂C µt = uq||rs δ pt + pu||rs δ qt + pq||us δ rt + pq||ru δ st , In the equations above, indices {ī,j, ...} refer to general occupied orbitals (both core and valence), while indices {i, j, ...} refer to valence occupied orbitals only.

Appendix E: Equations for the λ orbital response Lagrange multipliers
For CVS-ADC, the only non-zero blocks of λ are the mixed blocks: occupied-core, corevirtual 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. The supplementary material contains comparisons between analytical gradients and numerical gradients computed using the central finite difference algorithm (section I).
We list examples for the C, N, and O K-edge, respectively. The molecular geometries are from the supplementary material of Ref. [1]. The first C1s/N1s/O1s core-excited state of each molecule is given, where 10 singlet states were included in the ADC matrix diagonalization. The basis set used was def2-SVP. 2 The convergence tolerance in the ADC calculation was 10 −6 and the displacement step in the finite difference algorithm was 10 −4 bohr.
Additionally, we show several bending mode potential energy curves for the dissociative O 1s→ π * excited state of the water molecule, computed at different combinations of OH bond lengths (section II).