Efficient electronic structure calculation for molecular ionization dynamics at high x-ray intensity.

We present the implementation of an electronic-structure approach dedicated to ionization dynamics of molecules interacting with x-ray free-electron laser (XFEL) pulses. In our scheme, molecular orbitals for molecular core-hole states are represented by linear combination of numerical atomic orbitals that are solutions of corresponding atomic core-hole states. We demonstrate that our scheme efficiently calculates all possible multiple-hole configurations of molecules formed during XFEL pulses. The present method is suitable to investigate x-ray multiphoton multiple ionization dynamics and accompanying nuclear dynamics, providing essential information on the chemical dynamics relevant for high-intensity x-ray imaging.

INTRODUCTION

The advent of x-ray free-electron lasers (XFELs) opens up a new era in science and technology. Unprecedentedly ultraintense and ultrafast hard x-ray pulses generated from XFELs enable us to measure molecular structures on the atomic scale and to explore the structural dynamics on the femtosecond scale. One of the most prominent XFEL applications is imaging of biological macromolecules. X-ray crystallography with XFELs, after demonstration of the proof-of-principle, has started to reveal previously unknown protein structure, promising a breakthrough in structural biology (see reviews in Refs. 8–11). Recent advances in time-resolved serial femtosecond crystallography enable us to take a step forward towards molecular movies. To investigate molecular structure and structural dynamics with XFELs, one needs to understand radiation damage dynamics—x-ray multiphoton ionization dynamics and accompanying nuclear dynamics. Furthermore, the phase problem is the bottleneck to reconstruct molecular structures in serial femtosecond crystallography as much as in conventional x-ray crystallography. To overcome the phase problem for x-ray crystallography with XFELs, one uses conventional phasing technique at intermediate x-ray intensity, or one takes an advantage of the large degree of ionization at high x-ray intensity. The latter brings in high-intensity phasing (HIP) methods, including high-intensity multiwavelength anomalous diffraction and high-intensity radiation damage induced phasing. The HIP techniques require detailed description of deep-inner-shell decay dynamics of heavy atoms influenced by the molecular environment. Therefore, understanding of radiation damage dynamics is the key for successful molecular imaging.

Modeling of biological macromolecules exposed to XFEL radiation was pioneered by the seminal work of Neutze et al. Since then, there have been several computational tools to address molecular imaging problems, for example, classical molecular dynamics model, particle-in-cell approach, transport model, and Coulomb complex model. Some of these methods have been recently applied to ion fragment spectra and electron spectra from C60 molecules interacting with intense x-ray pulses. So far, most computational methods have been based on the independent-atom model or the plasma model. The description of the molecular environment is omitted by assumption or incorporated in an ad hoc manner. When a molecule absorbs x-ray photons, inner-shell multiple ionization induces fragmentation dynamics. Chemical bonds are weakened, and electrons and holes rearrange before the molecule breaks apart. Detailed electronic structure calculations for chemical bonding and charge rearrangement are thus crucial to describe molecular effects in modeling of the XFEL–matter interaction.

The electronic response of atoms and molecules to XFEL pulses is in essence characterized by multiphoton multiple ionization dynamics. The xatom toolkit has been developed to simulate the XFEL–atom dynamical interaction, and the ionization dynamics model has been tested with a series of experiments. The unprecedentedly large number of x-ray photons delivered by an XFEL pulse induces sequential x-ray absorptions, creating a variety of different q-hole configurations for each charge state . To simulate ionization dynamics, one needs to calculate photoionization cross section, Auger rate, and fluorescence rate for each configuration and solve a set of coupled rate equations for the time-dependent populations of the configurations. The q-hole electronic configurations are energetically highly excited, and theoretical treatment of such highly-excited states is challenging. Another challenge is the complexity of the ionization dynamics. Even for the atomic case, one must solve more than 2 × 107 coupled rate equations for Xe L-shell-initiated ionization dynamics. To address this formidable problem, a Monte-Carlo approach has been proposed for solving the rate equations and sampling the most probable configurations. In this scheme, the electronic structure is calculated for every single configuration selected on the fly as part of the Monte Carlo sampling. For the molecular case, the complexity increases even further because of the degrees of freedom associated with atomic motions, so the Monte Carlo approach seems to be the only way to simulate molecular response to an intense XFEL pulse. However, the most important question still remains: How to describe the electronic structure of molecules on the fly for the Monte Carlo method?

Here, we present a new efficient method for electronic structure calculation of polyatomic molecules and implement a dedicated toolkit, xmolecule. The proposed method is well suited for calculations of molecular multiple-hole configurations that are formed during x-ray multiphoton ionization dynamics. To efficiently describe molecular orbitals of core-hole configurations, the method employs atomic orbitals as basis functions that are numerical solutions of atomic core-hole states, calculated by xatom. For any given molecular electronic configuration and any given molecular geometry, xmolecule calculates molecular orbitals and orbital energies, which are essential components for dynamical simulations of x-ray multiphoton multiple ionization. We demonstrate that xmolecule is capable to calculate the whole spectrum of multiple-hole configurations at a given molecular geometry and potential energy surfaces for given multiple-hole configurations of molecules. Also, performance scalability with the system size is discussed. In this paper, we focus on the implementation of a molecular electronic-structure approach. Calculating cross sections and rates and solving coupled rate equations to simulate ionization dynamics will be described elsewhere. Having achieved these results, xmolecule aims to play a key role in molecular imaging at high x-ray intensity.

The paper is organized as follows. Section II formulates our scheme to calculate molecular multiple-hole configurations. It includes theoretical and computational schemes for basis function generation with numerical atomic orbitals, multicenter integration on a molecular grid, and direct Coulomb integral evaluation. In Sec. III, we show benchmark calculations for xmolecule, and then numerical results for the potential energy curves of various electronic configurations of carbon monoxide, and single- and double-core ionization potentials of several polyatomic molecules. We discuss the scalability of our scheme to a molecular size of hundreds of atoms. This is followed by the conclusion in Sec. IV.

COMPUTATIONAL METHODS

The Hartree–Fock–Slater method

We consider a molecular system composed of atoms with electrons. The Ath nuclear charge and coordinates are denoted by Z and , respectively. The molecular charge state is given by . We use the Hartree-Fock-Slater (HFS) method in which molecular orbitals (MO), , and orbital energies, , are obtained by solving the effective single-electron Schrödinger equation (atomic units are used unless specified otherwise), where is the external potential due to the nuclei, and the Hartree potential represents the classical Coulomb interaction among the electrons, and the last term represents the exchange interaction, which is approximated by the Slater exchange potential, The electronic density is obtained by the sum of squared MO's weighted by the occupation numbers as where . In contrast to conventional ground-state electronic structure calculations, in which the spin-orbitals with the lowest energies are filled, we consider all possible subject to , in order to take account of electronic excited states representing q-hole configurations.

The total energy within the HFS method is given by the sum of the nucleus–nucleus repulsion energy and the electronic energy,

Linear combination of numerical atomic orbitals

For atomic systems, the orbital is represented with spherical harmonics as where n, l, and m are the principal quantum number, the orbital angular momentum quantum number, and the associated projection quantum number, respectively. The radial wavefunction can be solved by a numerical grid-based method. The xatom toolkit has been developed to solve the atomic HFS equation. By employing the generalized pseudospectral (GPS) method and imposing a spherically symmetric potential, xatom accurately calculates for a given (n, l)-subshell, and accordingly for a given . This numerical atomic orbital has been used to successfully calculate multiple-hole configurations formed during x-ray multiphoton ionization dynamics in the atomic case.

For molecular systems, we employ the linear combination of atomic orbitals (LCAO) scheme to construct molecular orbitals, where is the μth atomic orbital (AO) and is the coefficient of the μth AO for the ith MO. Using Eq. (8) transforms the self consistent field (SCF) Eq. (1) into the corresponding Roothaan-Hall equation, where E is a diagonal matrix of MO energies and C is the MO coefficient matrix. The elements of the Hamiltonian matrix H and the overlap matrix S are given as where the effective potential . Equation (9) is solved in a self-consistent manner. To accelerate convergency, we employ the direct inversion in the iterative subspace (DIIS) method. When we encounter convergence problems at large bond distances, where the energy gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) is very small, we apply level shifts in the SCF iterations.

Here, our choice of basis set for the LCAO scheme is the numerical atomic orbitals (NAOs) obtained by xatom described above. In Fig. 1, we plot the squared radial function for the 1s, 2s, and 2p orbitals of the ground state of the neutral nitrogen (N) atom, the single-core-hole (SCH) state of , and the double-core-hole (DCH) state of , respectively. Comparison among different core-hole states shows significant deformation of valence orbitals in states with core holes. To cover these effects efficiently in the molecular calculation, we use NAOs that are numerical solutions of the corresponding atomic core-hole states. For instance, with one core hole at each atomic site (a DCH state) is calculated with basis functions optimized for on both N atoms, whereas with a single-site DCH state is calculated with basis functions optimized for on which the core hole is located and basis functions optimized for neutral N on the other side. In this way, we expect core-hole MOs are well described by core-hole-adapted NAOs.

FIG. 1.

Numerical atomic orbitals for different core-hole states of the nitrogen atom.

To achieve utmost efficiency towards complex ionization dynamics, we employ the minimal basis set. Each AO with (n, l, m) in Eq. (7) corresponds to a single basis function. Fully or partially occupied (n, l)-subshells contribute to a set of basis functions and each l gives basis functions (). For example, the N atom has 1s, 2s, and 2p (partially) occupied subshells, which constitute 5 basis functions (, and ) in total. This basis set is denoted as [2s1p]. According to the minimal-basis-set scheme, the chemical elements from B to Ne have the same number of basis functions (). In Sec. III A, we will discuss limitations and extensions of the minimal-basis-set scheme.

Molecular grid and multicenter integration

Equations (10) and (11) require evaluation of the corresponding integrals in three dimensions. In our case, the and are represented with a radial grid and spherical harmonics. To perform 3D integrals involving many atomic centers, we employ the multicenter integration proposed by Becke. Molecular grid points are constructed as a combination of sets of atomic grid points. Each set of atomic grid points, centered at one of the nuclei, consists of N radial grid points and angular grid points. The radial grid points are exactly the same as those used for NAO calculations with the GPS method. The angular grid points are obtained by the Lebedev grid scheme with an angular momentum cutoff at . The number of angular grid points is approximately given by . A detailed description of constructing multicenter molecular grid points is found in Refs. 60 and 61. We use an atomic radial grid size () large enough (∼10 Å) so that the atomic grids of many neighboring atoms overlap with each other. In principle, different atomic grid parameters can be used for individual atoms in a molecule. For convenience, however, we use the same grid parameters for all atoms. Then, the total number of molecular grid points is given by .

Becke's multicenter integration scheme introduces a set of smooth nuclear weight functions , subject to the constraint . The nuclear weight functions are generated by the third-order polynomial cutoff profile in the fuzzy cell scheme. Then, any integral of a given function f can be evaluated by the sum of individual atomic integrals, where . Each atomic integral can be readily performed using the spherical coordinate system of , centered at the Ath atom, where k is the index of the grid points of the Ath atom and w is defined as a product of the radial Legendre-Gauss-Lobatto quadrature weights and the angular Lebedev quadrature weights.

Implementation of direct Coulomb integrals

In electronic structure calculations, one of the most time-consuming parts is the evaluation of electron repulsion integrals. In order to achieve fast calculation within a desired accuracy, we have developed a multipole expansion scheme with an adaptive cut off. First, the integral involved in the Hartree potential in Eq. (3) can be decomposed into individual atomic integrals, where . Each single-center density can then be regarded as the atomic contribution to the total electronic density. To implement the integral, we expand the single-center density with real spherical harmonics as where is the (l, m)-component of the spherical expansion, With this single-center decomposition and spherical harmonic expansion of the electronic density, , the Hartree potential in Eq. (3) is obtained as where is given by where and . This radial integral is numerically evaluated in combination with various truncation methods (see the Appendix).

Molecular electronic configuration

Keeping the energetically lowest orbitals doubly occupied, the SCF procedure obtains the HFS solution for the electronic ground state. In order to obtain a solution for an excited electronic state of a q-hole configuration, each molecular orbital has to be assigned a specific occupation number. This can be done, as in the ground state calculation, by identifying the orbitals by their HFS energy eigenvalue. However, during the SCF iterations, the energetic order of MOs may change. Thus, identifying the orbitals by ordering them according to their HFS energy eigenvalue may lead to failure of the above SCF procedure or yield a solution for a different electronic state than required. This is called variational collapse.

To prevent this situation, we employ a variant of the maximum overlap method. In the maximum overlap method, the desired excited electronic state is specified by a set of initial guess orbitals in combination with a set of occupation numbers . In each SCF iteration, the occupation number n of the calculated orbital ψ is chosen according to its projection onto the subspace spanned by the guess orbitals with respective occupation number. Specifically, we calculate the overlap of the ith current MO with the jth guess MO, where . Note that the basis set for the initial guess orbitals is not necessarily the same as the one used for the expansion of the actual molecular orbitals, because different NAOs can be used for different q-hole configurations. Therefore, can be different from the overlap matrix defined in Eq. (11). Then, the projections of the ith orbital onto the span of the guess orbitals for the unoccupied (n = 0), singly occupied (n = 1), and doubly occupied (n = 2) cases are given by where j runs over all initial guess orbitals whose occupation number n equals n. To preserve the character of the required electronic configuration during the SCF procedure, we choose the set of the occupation numbers of the current orbitals, , such that is maximized, while the total number of doubly and singly occupied orbitals is maintained.

This procedure to determine the orbital occupation critically depends on the initial guess MOs. Thus, it is essential that the provided guess MOs together with the provided occupation numbers describe a wavefunction that is close to the required solution. For the calculations performed here, we choose initial guess MOs obtained from a previous calculation for a lower ionized electronic state or for the same electronic state with an altered molecule geometry. For the single-core-hole state in N2, we obtain a localized core hole on a specific nucleus by performing a Boys-orbital-localization procedure of the two guess core orbitals. Having obtained a converged solution, we verify that the obtained set of MOs is indeed close to the initial guess, by inspecting the individual overlap O.

RESULTS AND DISCUSSION

Benchmark calculations

We first estimate the accuracy of our calculations using the numerical multicenter integration in comparison with conventional calculations using the analytic Gaussian integration by gamess. Here, we employ the 6–31G Gaussian basis set to calculate the SCF-level ground-state energy of a water molecule. The internuclear distance of Å and the bond angle of are used. Only in this test, we employ the restricted Hartree-Fock (RHF) method instead of the HFS method, in order to directly compare with the gamess results. Figure 2 shows that our numerical calculations converge to the gamess results as the number of radial grid points per atom (N) and the number of angular grid points per atom (determined by ) are increased. The total number of molecular grid points for N = 50 and is . The maximum radius a.u., and the GPS mapping parameter L = 1 a.u. is used. Note that all grid parameters utilized here provide a numerical accuracy eV. If chemical accuracy is required (typically 1 kcal/mol 0.04 eV), our study for the water molecule shows that it is achievable with and , keeping the same L and . As to be shown in Sec. III B, the energy scale of x-ray-induced dynamics of highly-charged molecules will extend into the keV regime. Therefore, the worst grid parameters (for example, N = 30 and ) shown in Fig. 2 would be sufficient to describe the molecular ionization dynamics at high x-ray intensity.

FIG. 2.

Convergence of the total HF energy with respect to the number of grid points. N is the number of radial grid points, and controls the number of angular grid points per atom. The ground-state energy calculation of H2O with RHF/6-31G is performed using the numerical multicenter integration, and is the difference from the result obtained using the analytic Gaussian integrals.

We next examine the performance of our NAO basis set scheme. In Fig. 3, we show the calculated HFS energy of (a) the ground state of neutral N2 molecule with NAOs optimized for neutral N atom and (b) the quadruple-core-hole (QCH) state of ion with NAOs optimized for the DCH state of . The internuclear distance R = 1.096 Å is fixed. N = 200, L = 1 a.u., a.u., and are used. The results are shown together with those obtained by the equivalent calculations using conventional Gaussian-type-orbital (GTO) basis sets of different sizes (STO-3G and a series of Dunning's correlation-consistent basis sets; All GTO basis sets are obtained from the EMSL Basis Set Library). is the energy difference from the total energy calculated with the uncontracted version of cc-pV6Z, [16s10p5d4f3g2h1i] with 161 basis functions, which is considered here the complete basis set limit. Thus, indicates the numerical error due to lack of basis functions.

FIG. 3.

Comparison of convergency in total energy with respect to the number of basis functions, using the GTO scheme and the NAO scheme: (a) neutral N2 and (b) QCH . is defined by the total energy difference from the complete basis set limit (see the text).

In both Figs. 3(a) and 3(b), one can see that the minimal NAO basis set is superior to the conventional minimal basis set of STO-3G, illustrating that fully optimized NAOs are a practical choice for the basis set in the LCAO scheme. Also, Fig. 3 shows convergency of GTOs with respect to the number of basis functions. Interestingly, the conventional GTOs for QCH perform almost one order of magnitude less accurate than GTOs for neutral N2. The reason is that GTOs are optimized to be used for neutral ground-state calculations. In contrast, NAOs optimized for corresponding atomic q-hole configuration provide similar accuracy for both neutral N2 and QCH . Thus, NAO functions provide an ideal basis set for our minimal-basis-set HFS scheme.

To improve accuracy, we try to increase the number of NAOs in a systematic manner by including unoccupied atomic orbitals with higher (n, l) such as 3s, 3p, and so on. As shown in Fig. 3, the NAOs are somewhat inefficient to achieve higher accuracy by simply extending to higher (n, l), as previously reported in Ref. 73. This is attributed to the fact that additional series of higher (n, l)-orbitals, whose mean square radius is far from the atomic center, are inefficient for representing bonding molecular orbitals. Instead, we propose a scheme for adding compact p-type and d-type functions to the minimal NAO basis set in order to improve the description of chemical bonding. Additional functions are constructed by use of radial wavefunctions of occupied subshells multiplied by r, where r is the radial coordinate in the atomic system. For the chemical elements from B to Ne, the p-type functions are , where , and the d-type functions are , where . By adding these functions, as denoted by extended NAO (NAO[e]) and as marked with the black rectangle in Fig. 3, the accuracy is much improved; the total energy of neutral N2 is close to the cc-pVDZ level, and the total energy of QCH is close to the cc-pVQZ level. The number of basis functions for NAO[e] is only 13 per atom, whereas cc-pVQZ has 55 basis functions. There have been several approaches for extension of the minimal NAO basis set, where additional basis functions are constructed in a systematic way.

Potential energy curves for various hole configurations

Figure 4 shows the HFS total energies in Eq. (6) using core-hole-adapted NAO basis functions for all possible q-hole configurations that can be accessed by x-ray multiphoton ionization of the neutral carbon monoxide molecule. The internuclear distance R = 1.128 Å is fixed, and the grid parameters of N = 50, L = 1 a.u., a.u., and are used. For convenience, the figure shows these configurations grouped into charge states. The lowest horizontal line for each charge state indicates the ground-state energy for a given charge . This figure then illustrates how much energetically excited the q-hole configurations are. For example, the energy of DCH () is about 1 keV higher than the ground-state energy of . Ionization dynamics induced by intense x-ray pulses may occur step by step, visiting lots of these electronic states. Therefore, it is crucial to efficiently calculate this set of electronic states of q-hole configurations.

FIG. 4.

Spectrum of total energies for various electronic states of CO, which are accessible by x-ray multiphoton ionization. The colors indicate different core hole configurations. Potential energy curves of the DCH states inside the box will be shown in Fig. 5.

We further investigate the behavior of the potential energy curves (PEC) obtained using the NAO basis set. In Fig. 5, we show the calculated HFS total energies for three different types of DCH states: (a) , (b) , and (c) . The dashed red line indicates PECs calculated with core-hole-adapted NAO[e]. The solid red line indicates PECs with core-hole-adapted NAO without additional functions. Both results are compared with the solid blue line calculated with the conventional cc-pVTZ basis set. Previous theoretical studies of core-hole states suggested that calculations of the cc-pVTZ level are reasonably converged. Our NAO[e] scheme reproduces well the cc-pVTZ results, even though the size of NAO[e] () is much smaller than that of cc-pVTZ (). For comparison, we also plot PECs with NAOs optimized for neutral ground-state atoms, denoted by NAO[n], which shows a trend similar to what a conventional STO-3G minimal basis set would be. The NAO[n] results represent a poor estimate of PECs due to missing the core-hole effect on orbitals. On the other hand, PECs from NAOs, which are optimized for atomic core-hole states, show dramatic improvement over NAO[n], even though NAO and NAO[n] have the same number of basis functions ().

FIG. 5.

Potential energy curves of double-core-hole states as a function of the internuclear distance R. Energy is given relative to the ground-state energy of neutral CO.

Single- and double-core ionization potentials of molecules

To further test the accuracy of our calculation scheme, we compare core ionization potentials for a series of small molecules obtained from the HFS calculation. The molecular geometries are taken from Ref. 77 and the grid parameters are the same as those used in Sec. III B. We derive the single-core ionization potential from the HFS orbital energy of a neutral ground-state calculation using NAOs with and without additional basis functions. The double-core ionization potential is calculated as the sum of the first and the second core ionization potential, where the second core ionization potential is taken from the orbital energy of the SCH state calculation. For the DCH states with core holes on different nuclear sites, thus, two values are obtained for the two different ionization sequences.

Table I lists the ionization potentials compared with the values obtained from complete-active-space SCF (CASSCF) calculations and experimental values. For our calculations of two-site DCH states, a mean value and a deviation are listed for two values from the different ionization sequences. For the CASSCF results of two-site DCH states, only triplet spin states are listed and the difference between singlet and triplet states is smaller than 0.7 eV. Note that N2O is a linear molecule N—N—O, where N indicates the terminal N atom and N means N at the center. As can be seen, the CASSCF results show agreement within less than 4 eV with the available experimental values. The single ionization potentials we extract from the much simpler HFS calculation using the minimal NAO basis set show for all molecules a similar agreement within 5.1 eV, except in LiF (28 eV). For the DCH states, where the core holes are located on different nuclear sites, with the minimal NAO basis set, we also see a similar agreement within 7 eV to the CASSCF values and, where available, the experimental values. Again, LiF is an exception showing a much larger discrepancy of eV. For the DCH state with core holes on the same nucleus, we find a systematically larger disagreement of about 20–30 eV (for in LiF 78 eV).

TABLE I.

Single core hole and double core hole ionization potentials in eV. The molecular geometries are taken from Ref. 77.

Molecule	Configuration	Present (NAO)	Present (NAO[e])	CASSCF75	Exp.
CO	O1s^−1	537.43	533.80	542.82	542.5a
	C1s^−1	295.81	289.84	296.36	296.5b
	O1s^−2	1139.12	1136.43	1176.56
	C1s^−2	647.50	636.89	664.42	667.9b
	C1s^−1O1s^−1	850.70 ± 2.28	840.34 ± 0.52	855.20	855.3b
LiF	F1s^−1	663.77	670.23	688.04	691.8c
	Li1s^−1	59.34	58.56	65.33	61.9d
	F1s^−2	1403.81	1420.99	1481.50
	Li1s^−2	154.84	153.19	172.60
	Li1s^−1F1s^−1	735.13 ± 2.72	739.48 ± 1.06	763.28
N2	N1σg−1	409.57	403.30	411.03	409.9e
	N1σu−1	409.54	403.26	410.93
	N1s^−2	878.29	868.83	901.16	903.2f
	N1sA−1N1sB−1	836.96 ± 0.02	823.87 ± 0.02	836.44
N2O	O1s^−1	537.81	534.72	542.54	541.4g
	N_t1s^−1	408.70	403.66	408.61	409.0f
	N_c1s^−1	413.74	407.89	412.52	412.5e
	O1s^−2	1138.96	1136.23	1173.25
	N_t1s^−2	874.42	866.83	893.93
	N_c1s^−2	883.76	875.54	902.31
	O1s^−1N_t1s^−1	961.29 ± 0.25	951.47 ± 0.25	963.27
	O1s^−1N_c1s^−1	964.30 ± 0.40	954.53 ± 0.37	965.62
	N_t1s^−1N_c1s^−1	836.55 ± 0.01	825.25 ± 0.11	833.22	834.2f

Reference 78.

Reference 79.

Reference 80.

Reference 81.

Reference 82.

Reference 83.

Reference 84.

The inclusion of the p-type and d-type functions in the basis set leads in most cases to a larger deviation to the literature values than the results obtained with the minimal basis set. For these calculations, we get ionization potentials that tend to be lower than the literature values (from 3.4 eV for to 60.5 eV for in LiF). Clearly, the extended NAO basis set should improve the quality of the electronic structure model, as the electronic wave function has more flexibility. Thus, we conclude that the good agreement with the minimal basis set might be an artifact due to cancellation of errors.

For the results obtained with the larger basis set, we attribute the remaining deviations to the literature values mainly to relaxation energy contributions associated with the core hole electron removal. The applied scheme of taking orbital energies as ionization potentials cannot account for these effects. For core holes on the same nuclear site, where the core hole relaxation contributions are particularly strong, we see the strongest deviations (18.2–60.5 eV). Also, the extreme deviations for LiF may be explained from these contributions: The core hole on the F atom in LiF shows a particular large core hole relaxation effect, whereas for the core hole on Li it is very small.

Performance scaling

Our implementation of xmolecule aims for large-scale molecular calculations, especially for a large number of repeated calculations where time and resources available for each calculation are severely limited. At the same time, it requires the capability of calculating moderate-size systems in order to describe molecular-environment effects. Here, we demonstrate the performance scalability of our scheme toward molecular calculations with a few hundred atoms. Our grid-based method has the potential to achieve linear scaling in the number of atoms.

In the HFS method, the two-body interaction is divided into the exchange interaction and the direct Coulomb interaction. The former is replaced with the local density approximation, and the latter is treated with the Hartree potential as described in Sec. II D. The computational complexity of the Hartree potential is , where is linearly proportional to , because the potential in Eq. (3) contains the integral over molecular grid points and has to be evaluated at every single molecular grid point. By introducing the truncation methods described in the Appendix, this complexity can be reduced to . These truncation schemes do not change the quadratic scaling behavior with respect to , but reduce the actual computational time by several times (for example, a factor of two in our following calculations).

Another truncation can be made in the evaluation of one-body matrix elements in Eqs. (10) and (11). Both and are decomposed into atomic contributions by the multicenter integration: and . We define an AO pair and its contribution to each atomic grid, Then we set and to zero if , where ε is a control parameter. The complexity of the integrals in Eqs. (10) and (11) is , where both and are linearly proportional to . By using our truncation scheme described above, we can reduce it to a quadratic behavior with respect to .

Figure 6 shows the size dependence of the computation time of xmolecule with the current truncation schemes. We calculate the HFS ground state of C24H12 molecule (coronene) in its equilibrium molecular geometry taken from Ref. 77. Then, we perform calculations for n such molecules () stacked in the vertical direction with an interlayer separation of 3.3 Å. The minimal NAO basis set is used with N = 20, L = 1 a.u., a.u., and . The y axis is the CPU time per SCF iteration in seconds on a lab workstation (Intel Xeon X5660 2.80 GHz), and the x axis indicates the number of atoms in the stacked (C24H12) molecule. When all truncations are off (blue curve), the computational performance shows close to a cubic dependence. On the other hand, when the truncation method of Eq. (21) is applied with (red curve), the scaling shows a quadratic dependence on the system size. Note that when the truncation of Eq. (21) is used, the complexity of the matrix element calculations is reduced to a quadratic relation, while the Hartree potential calculation becomes the most time-consuming step, which is also governed by a quadratic scaling. The difference in the total energy between the calculations with and without this truncation is less than 0.14 eV/atom, whereas the truncated calculation is about 7.5 times faster than the calculation with no truncation. The calculation with 216 atoms (n = 6) takes 40 second per single SCF iteration on the lab workstation. The whole computation time takes about 14 min including the overhead costs for numerical grid construction and 12 SCF iterations. When additional truncation schemes for the Hartree potential (see the Appendix) are applied with and (green curve), the complexity is a bit reduced towards a linear relation and the errors in the total energy are less than 0.93 eV/atom. The actual computational time per iteration is improved by a factor of two for the 216-atom case.

FIG. 6.

Performance scaling with respect to the molecular size. The y axis is the CPU time per SCF iteration in seconds, and the x axis is the number of atoms in stacked (C24H12) molecules. The dotted lines with and indicate a quadratic behavior and a cubic behavior, respectively, with respect to the number of atoms, .

CONCLUSION

In summary, we present a new method to calculate various multiple-hole electronic states for polyatomic molecules that may be formed by x-ray multiphoton ionization dynamics at high x-ray intensity. The method is based on the Hartree-Fock-Slater method, employing the LCAO scheme, where NAOs are used as a minimal basis set for molecular orbital calculations. Usage of NAOs has two advantages over conventional Gaussian-type basis functions. First, NAOs are obtained from numerical solutions for atomic core-hole states at the same computational level. Second, accuracy and efficiency of numerical integration with NAOs are controllable by grid parameters and truncation schemes. The NAOs presented here are accurately solved by using the numerical grid-based method that is implemented in the xatom toolkit.

Using core-hole-adapted NAOs, molecular orbitals for core-hole states are efficiently calculated. We present benchmark calculations for multiple-core-hole states of N2. The NAO results show consistent accuracy for different charge states, which is not the case for conventional basis sets that are optimized for neutral systems. We demonstrate that our scheme is able to calculate all possible configurations that may be formed by removing zero, one or more electrons from the ground-state configuration of neutral CO molecule. The electronic state during x-ray multiphoton ionization dynamics may visit several of these multiple-hole configurations, which are energetically excited by about 4 keV with respect to the ground-state configuration of neutral CO. For molecular and ionization dynamics during XFEL pulses, we need not only all different multiple-hole states but also potential energy surfaces for individual electronic states. For double-core-hole states of , we calculate potential energy curves with core-hole-adapted NAOs, in good agreement with converged results with respect to the basis-set size. Also, we present single- and double-core-hole ionization potentials for several molecules in comparison with available theoretical and experimental data.

Efficient electronic structure calculations for molecules are essential for dynamical modeling of molecules at high x-ray intensity. We have implemented xmolecule to make a step toward dynamical simulation of molecular imaging with XFELs. Calculations of photoionization cross sections, fluorescence rates, and Auger rates for all possible configurations formed during molecular ionization dynamics are in progress.