PSIXAS: A Psi4 plugin for efficient simulations of X‐ray absorption spectra based on the transition‐potential and Δ‐Kohn–Sham method

Near edge X‐ray absorption fine structure (NEXAFS) spectra and their pump‐probe extension (PP‐NEXAFS) offer insights into valence‐ and core‐excited states. We present PSIXAS, a recent implementation for simulating NEXAFS and PP‐NEXAFS spectra by means of the transition‐potential and the Δ‐Kohn–Sham method. The approach is implemented in form of a software plugin for the Psi4 code, which provides access to a wide selection of basis sets as well as density functionals. We briefly outline the theoretical foundation and the key aspects of the plugin. Then, we use the plugin to simulate PP‐NEXAFS spectra of thymine, a system already investigated by others and us. It is found that larger, extended basis sets are needed to obtain more accurate absolute resonance positions. We further demonstrate that, in contrast to ordinary NEXAFS simulations, where the choice of the density functional plays a minor role for the shape of the spectrum, for PP‐NEXAFS simulations the choice of the density functional is important. Especially hybrid functionals (which could not be used straightforwardly before to simulate PP‐NEXAFS spectra) and their amount of “Hartree‐Fock like” exact exchange affects relative resonance positions in the spectrum.


| INTRODUCTION
X-ray absorption spectroscopy (XAS) is a powerful analytical tool to investigate molecules in the gas phase or adsorbed on surfaces, where an X-ray photon is used to promote the system into a core-excited state. In a one-electron picture, these states are characterized by a core hole, that is, an empty core orbital close to a nucleus that results from the excitation of a core electron into an unoccupied valence orbital. A Near Edge X-ray Absorption Fine Structure spectrum (NEXAFS) is obtained if the absorption cross section is measured as a function of the photon energy. [1] If the core excitation occurs after a valence excitation, so-called pump-probe-or transient NEXAFS (PP-NEXAFS) spectra can be measured. [2][3][4][5] For organic molecules containing N, C or O atoms, one is usually interested in the K-edge, meaning that a 1s electron is excited by an X-ray photon. Because core electron binding energies of different elements in a molecule are energetically well separated, it is possible to probe only a specific element in the molecule. In addition, if atoms of the same element have a different chemical environment, it is even possible to probe specific sites in the molecule.
Therefore, NEXAFS and PP-NEXAFS can be element-and sitespecific.
Therefore, only small to medium-sized molecules can be investigated. Alternatively, density functional theory can be used to simulate NEXAFS and PP-NEXAFS spectra. One way is to use techniques based on timedependent density functional theory, [29][30][31][32] another one uses the so-called transition-potential approach [33,34] in combination with the Δ-Kohn-Sham (Δ-KS) method. This approach has been successfully used by us [35,36] and other authors [34,37,38] to simulate NEXAFS absorptions of various systems.
Recently, we proposed and extension of this method to calculate also PP-NEXAFS spectra. [39] In this paper, we present a software plugin for the Psi4 quantum chemistry program [40,41] to calculate primarily NEXAFS and PP-NEXAFS phenomena using the transition-potential and Δ-KS techniques. Besides the core functionality to calculate X-ray absorption phenomena, the plugin offers a user-friendly input, a multicore parallelization and access to a large number of basis sets, generalized gradient, hybrid, meta and long-range corrected functionals, and an easy to use interface to visualize the resulting orbitals and spectra. Because of these features, we hope that this software plugin may be of use for other theoreticians or even for experimental researchers working in this field. We may note that other programs like CP2K, [42] StoBe [43] and Erkale [44] are also able to simulate ordinary NEXAFS spectra within the TP-DFT/Δ-KS approach. However, to the best of our knowledge, it is either not possible or rather cumbersome to calculate PP-NEXAFS using these programs. In addition, using hybrid, meta and long-range corrected functionals to simulate PP-NEXAFS spectra is a unique feature of PSIXAS.
In the following, we will briefly outline the theoretical foundations to calculate NEXAFS and PP-NEXAFS spectra, followed by details about the implementation. Then, we revisit the thymine molecule and analyze the importance of larger basis sets with higher angular momentum basis functions and the impact of hybrid, meta and long-range corrected density functionals on the shape of the PP-NEXAFS spectrum. Finally, we conclude our work and give an outlook.

| METHODS
As we will outline below, the simulation of NEXAFS and PP-NEXAFS spectra within the TP-DFT/Δ-KS approach requires self-consistent Kohn-Sham solutions with orbital occupation patterns different from the ground state. The most important occupation patterns for this work are shown in Figure 1. In the following two sections, we will briefly summarize the approach. For more detailed descriptions, we refer the reader to the original publications. [33,39]

| Ordinary NEXAFS simulations
Within the TP-DFT method, the NEXAFS cross section is calculated by using the Golden Rule expression Here, σ i gs ω ð Þ is the contribution of an excitation center i, C is a constant, ℏω the excitation energy, ϵ is the polarization of the incoming light, ψ i represents the initial state for excitation center i (see The final spectrum is usually further corrected by shifting the whole spectrum to the first core excitation energy Here, E gs is the ground state energy and E i ce is the first coreexcited state, that is, where one electron is fully removed from the core orbital and placed into the lowest, former empty, valence orbital. The occupation patterns that are used to find self-consistent solutions are shown in Figure 1 (denoted as E gs and E i ce ).
F I G U R E 1 Orbital occupation schemes used to simulate NEXAFS and PP-NEXAFS spectra within the TP-DFT and Δ-KS approach. The left three panels (framed in blue) correspond to the calculation of ground state NEXAFS spectra, while the three occupation patterns on the right (framed in red) are used to calculate PP-NEXAFS spectra. The index i labels the excitation center, where the core hole is created [Color figure can be viewed at wileyonlinelibrary.com] Finally, if multiple excitation centers are present, the overall spectrum is obtained by summing over the NEXAFS absorptions of all excitation centers as shown in Equation (1).

| Pump-probe NEXAFS simulations
In contrast to an ordinary NEXAFS experiment, where all molecules are excited from the electronic ground state, in a pump-probe experiment a certain fraction of molecules is pre-excited (pump) to a valence excited state before the X-ray absorption (probe) takes place. Further, a delay between pump and probe pulse allows the system to relax via internal conversion or intersystem crossing to other states which may be initially optically dark or spin forbidden. These states can then be probed by X-ray.
The PP-NEXAFS spectrum is therefore given by a weighted sum where p represents the aforementioned fraction of pre-excited molecules. σ i ve ω ð Þ is the cross section for the X-ray absorption that takes place at excitation center i for a valence excited molecule. The occupation pattern used to calculate the transition-potential state is shown in Figure 1 (σ i ve ω ð Þ ) and the spectrum is calculated according to Equation (1). In addition, a second energy correction is used to shift the spectrum to the lowest-energy transition where E ve is the energy of the valence excited state and its orbital occupation pattern is depicted in Figure 1 (denoted E ve and E ce ).

| IMPLEMENTATION DETAILS
Before we discuss the details of the implementation, it is worth to list the capabilities a program should have to successfully calculate NEXAFS and PP-NEXAFS spectra. A program needs to: 1 Perform ordinary unrestricted ground state Kohn-Sham calculations.
3 Perform a nonground state unrestricted Kohn-Sham calculation for an arbitrary occupation pattern.
4 Be able to maintain the occupation pattern during the selfconsistent field (SCF) procedure.
5 Offer flexibility to control the SCF convergence.
6 Calculate transition dipole matrix elements to obtain the spectrum.
7 Printout the orbitals for postprocessing.
As mentioned before, we chose to write a Python plugin for the Psi4 program. [40,41] The main advantage is that the plugin can make use of the comprehensive software infrastructure of Psi4 and its Python bindings. For example, typical tasks like user input validation, basis set assignment, the construction of the exchange-correlation functional, the calculation of two electron integrals, the calculation of one-electron matrices and many more are done by using Psi4 capabilities. This allows us to focus efficiently on the algorithms without spending too much time on technical issues. Another example project that makes use of this approach is the Psi4NumPy project, [45] which also inspired the development of PSIXAS.
It is therefore straightforward to implement an ordinary unrestricted Kohn-Sham algorithm using the Psi4 capabilities. The program can use all pre-and user-defined basis sets. All standard LDA, GGA, hybrid, meta, range-separated, and user-defined exchange-correlation functionals (including Hartree-Fock) can be used. Double-hybrid functionals are not supported. To control the convergence behavior, we use density mixing and a direct inversion of the iterative subspace (DIIS) algorithm. [46,47] To increase performance and decrease hardware requirements, we chose as default the density-fitting (RIJK approximation) algorithm, [48] but all other Psi4 options can also be used in the usual way. Once the calculation has converged, the orbitals are written for user inspection in the Molden format [49] and for restart/localization purposes in a NumPy format. [50] Then, the core orbitals can be localized (in case they are degenerate), using the Pipek-Mezey algorithm. [51] To calculate transition-potential, core-excited and valenceexcited states, we implemented a nonground state Kohn-Sham algorithm that is able to maintain a user-defined occupation pattern during the self-consistent field procedure. If for example, ψ γ n f g and f f γ n g with γ = α or β are the Kohn-Sham orbitals and occupation numbers, the electron density is given by where K is the number of basis functions/molecular orbitals. The orbital indices, orbital spins and orbital occupation numbers are defined with respect to the ground state solution, where f γ 1 , f γ 2 , …, f γ Nγ are equal to one, and f γ Nγ + 1 , …, f γ K are zero (N γ is the number of alpha/ beta electrons). Then, for the nonground state calculation the f f γ n g can take any values between zero and one. In addition, an overlap criterion can be applied to find the orbitals during the self-consistent field procedure in case the energetic order changes. For example, if in the beginning of the calculation the orbital ψ a,0 is chosen to have a different occupation number, its overlap during the SCF with the current set of orbitals {ψ j,n } can be calculated by where c a μ,0 and C νj,n are the coefficients for the basis set expansion (ψ a,0 = P μ c a μ,0 ϕ μ ) and S μν = 〈ϕ μ jϕ ν 〉 is the overlap matrix for the basis functions {ϕ μ }. The second index (0 or n) refers to the 0-th (initial) or n-th SCF step. Then the occupation number of that orbital having the largest overlap with the initial one is modified. This approach has some similarities with the so-called maximum-overlap-method (MOM) by Gilbert et al. [52] and its alternative, the initial-maximum-overlapmethod (IMOM) by Barca et al. [53] For both methods all occupied orbitals in SCF step n are those with the maximum overlap with the orbitals of the preceding n-1 step (MOM) or the orbitals of the initial guess (IMOM). Our implementation is different, as we use the overlap criterion for a subset of user-defined orbitals, while the rest is found by the aufbau principle.
Finally, one can choose orbitals to be frozen during the SCF. If the orbital with index j should be frozen, its off-diagonal elements in the Kohn-Sham matrix (F KS ) are set to zero in each step of the SCF Here, F KS ij = 〈ψ i jf KS j ψ j 〉 is the Kohn-Sham matrix in the basis of the molecular orbitals, withf KS being the Kohn-Sham operator. In addition to damping and DIIS techniques, we have found that energetically shifting the virtual orbitals is helpful in some cases. [54] Similar to the ground state Kohn-Sham procedure, orbitals are written for visualization and restart purposes. Once the transition-potential states are obtained in a self-consistent way, our program calculates and exports the transition dipole moments which can be used to calculate the absorption cross section. In addition to using the self-consistent transition state orbitals, it is also possible to use a the double basis technique as done by Ågren et al., [55,56] in case a better description of Rydberg and continuum states is required. Here, the basis set at the center where the core excitation takes place is augmented with a set of diffuse s, p, and d basis functions as used in StoBe [43] and Erkale. [44] The Kohn-Sham matrix is then once diagonalized in this new basis and the resulting molecular orbitals are used to calculate the transition dipole moments.
We introduced a set of keywords to provide sufficient flexibility and to control the flow of the program. All keywords and their options are listed in

| A PP-NEXAFS example
In the following, we present the results of a PP-NEXAFS simulation applied to the O-K edge of thymine, which was recently investigated by us. [39] For this system experimental spectra and simulations based on the correlated coupled cluster method exist. [14] The experimental peak positions will serve as reference for the simulations. The molecule and the labels for the two oxygens are shown in Figure 2. where the O2 1s !n resonance is significantly more intense (the O1 1s !n resonance will be neglected in the following).
Instead of discussing the detailed spectrum, we will focus primarily on the new features offered by PSIXAS. This includes the elaboration of larger basis sets with higher angular momentum functions and the use of more advanced density functionals for the simulation of the spectra. In the previous approach, [39] the choice of exchangecorrelation functionals was restricted to a limited set of local density and generalized gradient approximations. In addition, basis functions with angular momentum quantum numbers higher than two were neglected.
All spectra are calculated for a geometry, which was optimized using Psi4 and the same functional and basis set as for the PP-NEXAFS calculations. See Appendix and the supporting information for a set of example input files to calculate the PP-NEXAFS-spectrum for B3LYP/def2-SVP (see below). We assume a fraction of p = 13% of pre-excited molecules in a randomly oriented ensemble and use Gaussians with a full width at half maximum (FWHM) of 0.8 eV for the spectra shown in Figure 3 as done in Ehlert et al. [39] For the resonance positions given in Tables 2 and 3, we use the overall spectrum and extracted the corresponding local maxima of the peaks. We may note that relativistic effects are completely neglected in the current approach. As shown in a study by Takahashi and Pettersson [57] these effects are supposed to uniformly shift the O-K spectrum by approximately 0.3 eV toward higher energies.

| Basis set dependence
In Table 2, the basis set dependence of the aforementioned resonance positions is shown. All spectra simulations were done using the B3LYP [58,59] density functional and the Karlsruhe basis sets. [60] In addition to the absolute peak positions, the table also provides relative resonance positions.
From the numbers given in Table 2, one can derive the following conclusions: 1 All resonance positions are red shifted if larger basis sets are used.  Note: In the lower part of the table, the basis set has been decontracted on both oxygen atoms (labeled with "O-decon"). Geometries and spectra are obtained using the B3LYP functional. δ 1 and δ 2 represent the relative peak positions. Experimental values are taken from Wolf et al. [14] All resonance positions are given in eV. The labels for the oxygen atoms correspond to those in Figure 2. Note, that δ 1 and δ 2 are calculated in full accuracy before rounding.

T A B L E 3
Relative and absolute resonance positions in the (nπ*)-NEXAFS spectrum of thymine are calculated within the TP-DFT/Δ-KS approach O2 1s ! n O21s ! π * O1 1s ! π * δ 1 δ 2 Note: Different exchange-correlation functionals as implemented in Psi4 (with different amounts of "Hartree-Fock like" exact exchange, HF-X) have been used. For long-range corrected functionals two numbers are given, which correspond to the short and long range parts of the functional, respectively. For the IP-fitted version of the ωB97X functional an ω of 0.2442 is used (see text below). All calculations have been done with the def2-QZVPD basis set. Note, that δ 1 and δ 2 are calculated in full accuracy before rounding.

| Functional dependence
One additional feature of PSIXAS is the capability to make use of several hybrid, meta and long-range corrected density functionals within the TP-DFT/Δ-KS approach. Especially for the case of PP-NEXAFS spectra hybrid functionals can play an important role, because of their known influence on the energy gap between the highest occupied (HOMO) and lowest unoccupied (LUMO) molecular orbital. [61] We use all functionals as implemented through the LibXC library [62,63] in the Psi4 program. In Table 3, we list the absolute and relative resonance positions as a function of the exchange-correlation functional. The def2-QZVPD basis set has been used for spectra simulations and geometry optimizations, which have been done for each functional.
Note, that for HF and BHANDHLYP the energetic order of the O 1s MOs is exchanged compared to all other functionals, which show the same ordering as given in the example inputs (see Appendix).
The first observation is that the LDA functionals (SVWN, When going from TPSSh with an exact exchange of 10% to BHANDHLYP with an exact exchange of 50%, the difference between the O2 1s ! n and O2 1s ! π * resonances increases by 0.6 eV. However, the difference between the resonances arising from the O2 1s ! π * and O1 1s ! π * is also affected and increases by 0.4 eV, which slightly reduces the agreement with the experiment. Therefore, the amount of exact exchange is a key ingredient to model PP-NEXAFS spectra accurately within the TP-DFT/Δ-KS approach.
These findings are illustrated in Figure 3 for selected global hybrid functionals with different amounts of exact exchange, that is, BLYP, B3LYP, PBE0, B97-K, and BHANDHLYP. In Figure 3a the simulated PP-NEXAFS spectra are given and vertical blue lines indicate the experimental resonance positions. All spectra seem to be rather similar and scatter around the experimental values with acceptable error margins. However, in Figure 3b the same spectra are shown, but this time shifted in a way that the O2 1s ! π * resonance is centered around 0 eV. In this graph one can clearly see that δ 2 , that is, difference between the O2 1s ! n and O2 1s ! π * resonances, indeed depends on the amount of exact exchange, as stated above. In addition the increase of δ 1 , that is, the splitting of the two O 1s ! π * resonances, as a function of the exchange can be observed. For this selection of functionals, we find that B97-K and BHANDHLYP perform best with respect to relative resonance positions. The latter also provides the best absolute resonance positions compared to the experiment in this group of functionals.
One can also examine the dependence of simulated resonance positions as a function of the amount of exact Hartree-Fock like exchange by constructing user defined functionals. This is shown in Note: A set of user defined exchange-correlation functionals based on the BLYP functional has been used. All calculations have been done with the def2-QZVPD basis set and at the optimized BLYP geometry. Note, that δ 1 and δ 2 are calculated in full accuracy before rounding.
distances, which is controlled by a parameter ω. Within the so-called IP-fitting procedure, an ω is found, such that the negative energy eigenvalue of the HOMO is equal to the ionization potential obtained by the energy difference of the cation and neutral molecule, which a property of the exact functional. [64] We therefore used the IP-fitting procedures implemented in Psi4 to find an improved ω for the ωB97X functional and used it thereafter to compute the PP-NEXAFS spectrum (ω was found to be 0.2442). The results are also shown in Table 3. Although ωB97X provides already resonance positions in very good agreement with the experiment, ωB97X-IPFIT still slightly improves the results. All three resonances are blue shifted between 0.1 and 0.3 eV after the IP-fitting. Therefore, IP-fitting might be a useful system dependent tool to obtain even more accurate spectra.
However, more elaboration on this aspect is required, but should be left out for future explorations.

| SUMMARY AND OUTLOOK
In summary, we present a plugin for the Psi4 software package which allows the simulation of NEXAFS and PP-NEXAFS spectra within the TP-DFT/Δ-KS approach. The core capability of the plugin is to find self-consistent nonground state Kohn-Sham solutions for a given occupation pattern. This is accomplished by overlap and freezing techniques. In addition, the plugin provides access to a wide selection of basis sets and LDA, GGA, META, LRC, and hybrid density functionals.
Applied to the case of thymine, we find that larger basis sets are needed to provide a better description of electron relaxation effects.
As a result, better absolute resonance positions can be obtained, while relative peak positions are not affected. On the other side, the use of functionals that incorporate exact exchange has a direct influence on the relative resonance positions in the PP-NEXAFS spectra. For thymine, best results with respect to relative peak positions have been obtained with the (IP-fitted) ωB97X and BHANDHLYP functionals, the latter having a rather high amount of global exact exchange. It remains to future investigations to find out if these functionals are in general more accurate for PP-NEXAFS simulations, or if the amount exact exchange needs to be adjusted for each individual system.
In the near future, we will focus on expanding the capabilities and functionalities of PSIXAS. On the DFT side, a recent paper by Hait et al. [65] demonstrates good accuracy by using the restricted openshell Kohn-Sham (ROKS) method in combination with the square gradient minimization. [66] On the side of wave function based methods, we will focus on efficient techniques based on configuration interaction [17] and extend them for PP-NEXAFS spectra.

ACKNOWLEDGMENT
The authors thank M. Gühr and P. Saalfrank for fruitful and stimulating discussions.