Calculation of Core-Excited and Core-Ionized States Using Variational Quantum Deflation Method and Applications to Photocatalyst Modeling.

The possibility of performing quantum-chemical calculations using quantum computers has attracted much interest. Variational quantum deflation (VQD) is a quantum-classical hybrid algorithm for the calculation of excited states with noisy intermediate-scale quantum devices. Although the validity of this method has been demonstrated, there have been few practical applications, primarily because of the uncertain effect of calculation conditions on the results. In the present study, calculations of the core-excited and core-ionized states for common molecules based on the VQD method were examined using a classical computer, focusing on the effects of the weighting coefficients applied in the penalty terms of the cost function. Adopting a simplified procedure for estimating the weighting coefficients based on molecular orbital levels allowed these core-level states to be successfully calculated. The O 1s core-ionized state for a water molecule was calculated with various weighting coefficients, and the resulting ansatz states were systematically examined. The application of this technique to functional materials was demonstrated by calculating the core-level states for titanium dioxide (TiO2) and nitrogen-doped TiO2 models. The results demonstrate that VQD calculations employing an appropriate cost function can be applied to the analysis of functional materials in conjunction with an experimental approach.

Introduction

The concept of performing quantum-chemical calculations using quantum computers has attracted significant attention.[1,2] In this regard, the configuration interaction (CI) framework is the most common means of assessing electron correlation. In this framework, the wave function is described as a linear combination of Slater determinants or spin-adapted configuration state functions that represent the electron configurations.[3] If all of the electron configurations generated from combinations of all the electrons and orbitals of the system are included, the method is termed full CI.[3,4] Although exact wave functions within an adopted basis set can be obtained using full CI calculations, the computational cost increases quite rapidly as the number of orbitals and electrons increases and can easily outgrow the capacity of standard computers. In contrast, quantum computers could potentially perform the required calculations in polynomial time using the quantum-phase estimation (QPE) algorithm.[5−8] Unfortunately, the QPE method requires fault tolerance and cannot be practically executed on current quantum computers (so-called noisy intermediate-scale quantum (NISQ) devices) that operate on the basis of qubits without error correction.[9−11] Accordingly, a quantum-classical hybrid algorithm (termed the variational quantum eigensolver (VQE)) has emerged as a means of allowing quantum-chemical calculations to be performed using NISQ devices.[12,13]

Following the development of the VQE method to compute the ground state for a given system, VQE-based algorithms for excited states have also been proposed, including quantum subspace expansion (QSE),[14] multiscale-contracted VQE (MCVQE),[15] subspace-search VQE (SSVQE),[16] quantum equation of motion VQE (qEOM-VQE),[17] and others. The variational quantum deflation (VQD) method introduced by Higgott et al. is one such algorithm and allows computation of excited states based on the ground state obtained using the VQE method.[18] The superior efficiency and accuracy of the VQD approach compared with the SSVQE and MCVQE methods have been demonstrated based on calculations of the excited states for lithium hydride, diazene, and azobenzene.[19] The VQD method has also been used to investigate the excited states for phenylsulfonylcarbazole compounds during the development of thermally activated delayed fluorescence (TADF) emitters for organic light emitting diode (OLED) applications.[20] These studies suggest that the VQD method could be employed in the research and development of functional materials such as light emitting devices,[21] photocatalysts,[22] photovoltanic cells,[23] and photochromic materials.[24]

In VQD calculations, the excited states are sequentially computed in conjunction with minimization of the defined cost function.[18] Including penalty terms within the cost function is the most general means of ensuring that the computed wave function will be orthogonal to other states and have the desired eigenvalues.[13,18,25,26] Penalty terms are introduced to raise the cost in the case that ansatz states are nonorthogonal to other states or eigenvalues deviate from user-specified values. Thus, the weighting coefficients that determine the magnitudes of the penalty terms dictate whether the desired excited states are obtained or not, and so these coefficients must be carefully selected. Even so, the coefficients are typically chosen heuristically. Although small coefficients are preferable so as to obtain rapid convergence, undesired excited states can possibly be obtained if the coefficients are too small. Conversely, if the coefficients are too large, the optimization may not converge properly. To address this problem, Kuroiwa and Nakagawa proposed a convenient formula to estimate the weighting coefficients as a part of their systematic analysis of the penalty terms.[27] Assuming that this formula is generally applicable to calculations of the excited states for common molecules, it could greatly increase the practical applications of the VQD method. Although this formula was initially assessed on the basis of calculations of simple model systems such as H2 and H4, it should be examined to determine if it is more widely applicable to the excited states for more complex molecules. A practical procedure for utilizing this method should also be explored to prepare for upcoming applications of the VQD method using quantum computers.

On this basis, the present study applied the VQD method to calculate the core-excited and core-ionized states for common molecules based on simulated quantum circuit computations with a classical computer. An analysis using a classical computer is advantageous for this study because the results can be examined without being influenced by noise. The absorption bands derived from the core-excited state calculations were experimentally verified using X-ray absorption near edge spectroscopy (XANES),[28−31] while the core-ionization energies were confirmed using X-ray photoelectron spectroscopy (XPS).[32−35] It should be noted that the core-ionization energy is also referred to as the core–electron binding energy (CEBE) and that XPS is also traditionally known as electron spectroscopy for chemical analysis (ESCA). These spectroscopic techniques are widely utilized for the analysis of molecular structures and to determine the surface compositions of condensed materials, local bonding environments of specific atomic species, and electronic properties such as oxidation states. Computations of core-level states are of interest with regard to fundamental scientific studies and also with respect to the theoretical analyses of spectra.[36,37] If the calculation of core-level states using the VQD method is found to be possible, this technique could have practical applications in conjunction with experimental analyses. As part of this study, the effects of the weighting coefficients used in the penalty terms were systematically analyzed by calculating the O 1s core-ionized state for a water molecule. The core-excitation energies and core-ionization energies for this molecule are on the order of hundreds of eV or higher and so are much larger than those of the more well-studied valence excited states. Therefore, the correlation of the weighting coefficients with the results could be clearly visualized. As examples of applications to functional materials, the O 1s core-excited states and Ti 2p core-ionized states for titanium dioxide (TiO2)[38−41] and nitrogen-doped TiO2 (N-TiO2)[42,43] were also calculated. These materials were selected because anatase TiO2 is well-known for its photocatalytic activity and N-TiO2 is a visible-light-sensitive photocatalyst.

Note that singly core-excited and singly core-ionized states were calculated in this study. Because these states have less multiconfiguration character, their wave functions can be reasonably described using a small number of electron configurations. Thus, even though only a small number of qubits can be currently employed in simulations on classical computers (as with calculations using actual NISQ devices), the excitation and ionization energies in such cases can be evaluated with high accuracy so that the calculation results can be verified by comparison with experimental data.[44] It is noteworthy that Bauman at al. recently calculated the doubly excited core-level states for a water molecule using the QPE algorithm.[45]

Theory

VQD Method and Weighting Coefficients of Penalty Terms

Here, the VQD algorithm is briefly reviewed[18] as well as the penalty terms related to the cost function.[27] In the VQD algorithm, the kth ansatz state, |ψ(θ)⟩, is calculated by minimizing the cost functionHere Ĥ is a Hamiltonian written as Ĥ = ∑cP̂, where P̂ is a single-qubit Pauli operator and c is the corresponding coefficient. We denote the eigenstates up to the (k – 1) th as |ψ(θ)⟩ (i = 0, 1, ..., k – 1) and assume that |ψ(θ)⟩ has been obtained before the calculation of |ψ(θ)⟩. The ground state, |ψ(θ0)⟩, is initially obtained using the VQE method, after which the |ψ(θ1)⟩ state is calculated. In a similar manner, the |ψ(θ2)⟩, |ψ(θ3)⟩, ... and |ψ(θ)⟩ states are calculated one by one in a sequential fashion. The second term of eq is added to ensure that |ψ(θ)⟩ is orthogonal to |ψ(θ)⟩, while β is a weighting coefficient referred to as the overlap weight. Unless |ψ(θ)⟩ is orthogonal to one of the |ψ(θ)⟩ states, L(θ) is increased by β because of this second term. The third term, Lpenalty, is introduced to impose a constraint on the spin multiplicity, spin quantum number and number of electrons associated with |ψ(θ)⟩. The Lpenalty adopted in this study consisted of three terms, written aswhere Ŝ2 is the total spin-squared operator, Ŝ is the z-component of the total spin operator, Ŝ, N̂ is the particle-number operator, w1, w2, and w3 are weighting coefficients termed the s2 number weight, sz number weight, and particle number weight, respectively, and o1, o2, and o3 are the eigenvalues of Ŝ2, Ŝ, and N̂, respectively, that must be satisfied by |ψ(θ)⟩. Here, the o values (x = 1–3) are input parameters. Lpenalty goes to zero when the eigenvalues of |ψ(θ)⟩ coincide with the user-specified o values. In contrast, if the eigenvalues of |ψ(θ)⟩ deviate from the o values, L(θ) increases because Lpenalty is non-zero. When the operators in eq act on |ψ(θ)⟩, we havewhere S is the total spin angular momentum, MS is the spin quantum number, and N is the number of electrons. If we define Δ as the deviation of the eigenvalues of |ψ(θ)⟩ from o, we obtainwhere Δ1 = S(S + 1) – o1, Δ2 = M – o2, and Δ3 = N – o3.

Estimation of Weighting Coefficients for Excited-State Calculations

Estimation of β

If |ψ(θ)⟩ is orthogonal to |ψ(θ)⟩ and all the Δ values are zero, L(θ) is equal to the energy of |ψ(θ)⟩, meaning that

Next, we assume an ansatz state |ψ(θ*)⟩ that is nonorthogonal to one of the |ψ(θ)⟩ states. Consequently, the Δ values of |ψ(θ*)⟩ are zero and its cost function, L(θ*), can be written asif L(θ) > L(θ*), |ψ(θ*)⟩ is obtained instead of |ψ(θ)⟩. To obtain |ψ(θ)⟩, β must satisfy the conditionso that L(θ*) > L(θ). The right-hand side of eq is the energy gap between |ψ(θ)⟩ and |ψ(θ*)⟩, which is maximized if |ψ(θ*)⟩ is the ground state: |ψ(θ0)⟩. We then have the condition for β written as

This condition was provided in the original paper in which VQD was described.[18] Unfortunately, the right-hand side of eq is the excitation energy of |ψ(θ)⟩, which is a part of the calculation results. As a practical alternative, the excitation energy of the singly excited state can be estimated from the energy gap between the occupied and virtual orbitals related to the excitation aswhere ε and ε are the energy levels of the occupied and virtual orbitals that are relevant to the main configuration of |ψ(θ)⟩, respectively.

In this study, the β values were estimated based on eq , while the ε and ε values were obtained through the Hartree–Fock calculation that was used to obtain an initial state for the ground-state VQE calculation. In general, the gap between ε and ε will be larger than that between E(θ) and E(θ0) when determined using the Hartree–Fock method.[46] Consequently, the β value estimated from eq should satisfy eq . Note that the value of ε – ε based on the use of Kohn–Sham orbitals may be smaller than E(θ) – E(θ0). In addition, calculations using the density functional theory (DFT) method adopting an exchange-correlation functional tend to give a smaller energy gap between the occupied and virtual orbitals compared with the corresponding excitation energy.[46]

In this study, we adopted a procedure using the orbital energies of the ground state to estimate β and other weighting coefficients as mentioned above and below. As a matter of course, more sophisticated methods can also be applied to estimate the excitation energies. For example, Kuroiwa and Nakagawa used the configuration interaction singles and doubles (CISD) method to compute the excitation energies.[27] The estimation based on the orbital energies of cationic species may also be practically useful.[47,48]

Estimation of w1

The values of w (x = 1–3) were estimated using the formula proposed by Kuroiwa and Nakagawa.[27] These calculations assumed |ψ(θ*)⟩ with Δ ≠ 0 and Δ = Δ = 0, where (p, q, r) ∈ x, p ≠ q, p ≠ r and q ≠ r. In addition, |ψ(θ*)⟩ was orthogonal to |ψ(θ)⟩ because Δ ≠ 0. In this case, L(θ*) could be written as

According to eqs and 10, to have L(θ*) > L(θ)must be true. Therefore, the lower limit of w is obtained when |Δ| is at a minimum, |Δmin|, such thatwhere Δmin is the minimum deviation of the eigenvalue from the user-specified o. Equation indicates that the lower limit of w can be estimated based on E(θ) – E(θ*) and Δmin.

We also assumed that, in the case of p = 1, |ψ(θ*)⟩ has Δ1 ≠ 0 and Δ2 = Δ3 = 0 and |ψ(θ*)⟩ is orthogonal to |ψ(θ)⟩ when Δ1 ≠ 0. This gives a L(θ*) value equal to

Because the eigenvalue of Ŝ2 is S(S + 1) and the minimum change in S is 1/2, |Δ1min| is 3/4. Hence, according to eq should be true. The right-hand side of eq is maximized when E(θ*) = E(θ0), such that

The energy gap between E(θ) and E(θ0) is the excitation energy for E(θ), which can be approximated based on ε – ε in a similar manner to eq as

In the present work, eq was adopted as the condition for w1.

Estimation of w2

We assumed the case of p = 2, meaning |ψ(θ*)⟩ with Δ2 ≠ 0 and Δ1 = Δ3 = 0 and |ψ(θ*)⟩ orthogonal to |ψ(θ)⟩ when Δ2 ≠ 0. L(θ*) in this case is

Since the eigenvalue of is M and the minimum change in MS is ±1/2, |Δ2min| is 1/2. According to eq , the w2 value should satisfy w2 > 4(E(θ) – E(θ*)) so that L(θ*) > L(θ). The gap between E(θ) and E(θ*) is maximized when E(θ*) = E(θ0), and the gap between E(θ) and E(θ0) can be approximated by ε – ε as with eq , such that

The w2 values in this study were estimated on the basis of eq .

Estimation of w3

We assumed the case of p = 3, meaning |ψ(θ*)⟩ with Δ3 ≠ 0 and Δ1 = Δ2 = 0 and |ψ(θ*)⟩ orthogonal to |ψ(θ)⟩ when Δ3 ≠ 0. L(θ*) in this case is

Because the eigenvalue of is the number of electrons, its minimum change is ±1 and, accordingly, |Δ3min| is 1. Therefore, the w3 value should satisfy w3 > E(θ) – E(θ*) so that L(θ*) > L(θ). The gap between E(θ) and E(θ*) is maximized when E(θ*) = E(θ0), and the gap between E(θ) and E(θ0) can be approximated by ε – ε. Thus, we have

The w3 values in this study were estimated on the basis of eq . Overall, the weighting coefficients for the calculations of core-excited states were estimated so that

Estimation of Weighting Coefficients for Core-Ionized State Calculations

In the present study, each core-ionized state was calculated as a high-energy excited state for a molecule in its ionized state. In this procedure, the valence-ionized state was obtained as |ψ(θ0)⟩ via the VQE calculations. Accordingly, the charge-neutral ground state having a lower energy than |ψ(θ0)⟩ could be obtained as |ψ(θ*)⟩ and L(θ*) in this case was

If |ψ(θ)⟩ is the core-ionized state, Lpenalty should satisfy the conditionsuch that L(θ*) > L(θ).

The right-hand side of eq is the core-ionization energy resulting from the VQD calculation. This energy can be approximated as the negative of the core-orbital energy (ε) according to Koopman’s theorem[49]

Equation indicates that the weighting coefficients used in calculations of core-ionized states should be estimated based on – ε instead of ε – ε, in contrast to the calculations for the excited states. Accordingly, this work adopted the conditionand the ε values obtained from the Hartree–Fock calculations were used.

Computational Details

Optimizations of molecular geometries via coupled cluster singles, doubles, and perturbative triples (CCSD(T))[50−52] calculations and DFT calculations with the B3LYP functional[53,54] were carried out using the Gaussian 09 program.[55] The VQE and VQD calculations with the optimized structures were performed using the Qamuy program version 0.26.1[56] in conjunction with a hardware-efficient ansatz.[57] The Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm was used to update variational parameters on a classical computer,[58] while the Jordan–Wigner transformation was used to map Hamiltonians into a qubit circuit.[59,60] The depth of the simulated quantum circuit (D) was given a value of 10. Noise was not considered with regard to the simulated quantum circuit, meaning that exact expected values were used in the simulations.

Because twice as many qubits as the number of active orbitals were required within the framework detailed above, the computational cost of the simulated calculations on a classical computer increased rapidly with increasing number of active orbitals. Accordingly, the active space of the VQD calculations had to be constructed using as few active orbitals as possible. In addition, the molecular orbitals were greatly perturbed by core hole formation. A core hole attracts electrons, leading to radial contraction of the outer orbitals, such that the core-hole state becomes highly stable. This effect is referred to as orbital relaxation,[61] and the effects of orbital relaxation were not fully incorporated into the CI calculations using a limited number of electron configurations. In such cases, orbital optimization can be an effective means of improving the accuracy of the calculations,[44] and so VQD calculations with and without molecular orbital optimization were compared. The calculations with and without orbital optimization corresponded to conventional complete active space self-consistent field (CASSCF)[62,63] and complete active space configuration interaction (CASCI) calculations,[64] respectively. A flexible basis set was also required to deal with the significant changes of the orbitals and Dunning’s cc-VXZ basis set (X = D, T, Q) was adopted.[65−67] The weighting coefficients were estimated based on the Hartree–Fock orbital energies with the cc-pVDZ basis set.

Calculations of Core-Excited States for CO, H2CO, and HCN

The 1s → π* core-excited states for CO, H2CO, and HCN molecules were calculated. The molecular geometries were optimized using the CCSD(T)/cc-pVQZ approach and the core-excited states were calculated using the VQD method with the cc-pVXZ basis set. The 1s → π* core-excitation energies were obtained as the energy gaps between the ground and core-excited states. In the case of CO and H2CO, both the C 1s and O 1s core-excited states were calculated. Similarly, both the C 1s and N 1s core-excited states were calculated for HCN. The 1s and π* orbitals were selected as the active orbitals, and two electrons in these orbitals were treated as active. In the following discussion, the CAS constructed from these active orbitals and electrons is denoted as CAS(2e, 2o). During the calculations for CO and HCN, two 1s → π* core-excited states were found to be degenerate because two π* orbitals had the same energy level, and so only one state was calculated. For H2CO, the π* orbitals were split into different energy values because of the lower symmetry of the molecule. The weighting coefficients were estimated from the energy gaps between the 1s and π* orbitals in conjunction with eq .

Calculations of Core-Ionized States for CH4, NH3, H2O, and FH

The 1s core-ionized states for CH4, NH3, H2O, and FH molecules were calculated. The molecular geometries were optimized using the CCSD(T)/cc-pVQZ method, and the core-excited states were calculated using the VQD method with the cc-pVXZ basis set. The 1s core-ionization energies were obtained as the energy gaps between the charge-neutral (closed-shell) singlet ground state and doublet core-ionized state. The 1s orbital and the highest occupied orbital (HOMO) were selected as the active orbitals and three electrons in these orbitals were treated as active electrons. Henceforth, the CAS constructed from these active orbitals and electrons is denoted as CAS(3e, 2o). When adopting CAS(3e, 2o) as the configuration space, the valence-ionized state having the configuration (1s)2(HOMO)1 was the ground state and the core-ionized state was the excited state having the primary configuration (1s)1(HOMO)2. Accordingly, the core-ionized states could be calculated using the VQD algorithm. In a full CI calculation using a quantum computer with a larger number of qubits, all the occupied and virtual orbitals would be treated as active orbitals, such that the core-ionized states would be obtained as the excited states. The calculations presented herein based on CAS(3e, 2o) simulated this type of full CI calculation. The weighting coefficients for the VQD calculations were estimated based on the 1s orbital energies and using eq .

Effects of the Weighting Coefficients

The relationship between the weighting coefficients and the ansatz states was examined by calculating the 1s core-ionized state for a H2O molecule while varying the weighting coefficients. The molecular structure was optimized via a CCSD(T)/cc-pVQZ calculation while VQD calculations were carried out using the cc-pVDZ basis set.

Application to the TiO2 and N-TiO2 Models

Preparation of the Models

In the quantum-chemical calculations of solid materials such as metals and metal oxides, a small cluster of atoms is often employed as a model.[68] The calculation results obtained using this model cluster approach are significantly affected by the size and shape of the cluster, as well as the manner in which dangling bonds are terminated, all of which are arbitrarily assigned. Nakatsuji et al. proposed three principles for the construction of clusters: neutrality, stoichiometry, and coordination.[69] In this work, we adopted a (TiO2)16 cluster, which is known to represent the smallest cluster of anatase-type TiO2 that satisfies these three conditions.[70] Specifically, this cluster is charge neutral, has a stoichiometric Ti:O ratio of 1:2 (equal to that of the bulk system), and has no dangling bonds. The coordination numbers of the titanium atoms in this cluster were four or more, while those of the oxygen atoms were two or three. We initially prepared a (TiO2)16 cluster and optimized the geometrical structure using the B3LYP/cc-pVDZ basis set. Following this, a smaller cluster with one titanium atom and four oxygen atoms was cut out from the surface of this (TiO2)16 cluster and its dangling bonds were terminated with hydrogen atoms to obtain Ti(OH)4. The positions of the four hydrogen atoms in the Ti(OH)4 were reoptimized using the B3LYP/cc-pVDZ basis set, holding the other atoms at fixed positions during the reoptimization. Finally, one of the oxygen atoms of Ti(OH)4 was replaced with a nitrogen atom and the newly generated dangling bond on the nitrogen atom was saturated with a hydrogen atom to obtain Ti(OH)3(NH2). The positions of the hydrogen atoms of Ti(OH)3(NH2) were reoptimized using the B3LYP/cc-pVDZ basis set with the positions of the other atoms fixed. The resulting Ti(OH)4 and Ti(OH)3(NH2) clusters were employed as models of TiO2 and N-TiO2, respectively.

Calculation of Valence Excited States, O 1s Core-Excited States, and Ti 2p Core-Ionized States

The HOMO → LUMO excited states, O 1s → LUMO core-excited states and Ti 2p core-ionized states for Ti(OH)4 and Ti(OH)3NH2 were calculated utilizing the VQD algorithm. The HOMO and LUMO of a solid material correspond to the valence band maximum (VBM) and the conduction band minimum (CBM), respectively. Therefore, absorption band edge shifts can be estimated from HOMO → LUMO excitation energies. The peak with the lowest energy in the O K-edge XANES spectrum was derived from the O 1s → LUMO excitation, and the corresponding peak in the XPS spectrum was situated at the Ti 2p core-ionization energy. The effects of nitrogen doping on these spectroscopic features could be examined based on calculations involving the Ti(OH)4 and Ti(OH)3(NH2) cluster models.

The CAS for the calculations of the HOMO → LUMO excited states was constructed using these two orbitals and the two electrons in the HOMO. Similarly, the CAS for the calculations of the O 1s → LUMO excited states was constructed using these two orbitals and the two electrons in the O 1s orbital. The CAS adopted for these calculations is denoted herein as CAS(2e, 2o). During the calculations of the Ti 2p core-ionized states, one of the Ti 2p orbitals and the HOMO were selected as the active orbitals and three electrons in these orbitals were treated as active electrons. This CAS is denoted as CAS(3e, 2o). Since there are three Ti 2p orbitals, three VQD calculations were made for each model, employing a different Ti 2p orbital as an active orbital.

The energy levels of the Ti 2p core-ionized states are split into 2p1/2 and 2p3/2 levels because of spin–orbit interactions.[71] Unfortunately, the ab initio algorithm for relativistic effects is not presently implemented in the Qamuy program, and so this splitting had to be approximated. Although the gap between the 2p1/2 and 2p3/2 states (ΔSO) of Ti varies depending on the material, the value is typically within the narrow range of 5.60–6.13 eV,[72−74] and so the ΔSO value was assumed to be 5.7 eV in the present work.[73] In the case of the core-ionized doublet state for an atom, the 2p1/2 energy level will be lower by (3/2)ΔSO than the 2p level, while the 2p3/2 energy level will be higher by (1/2)ΔSO (Figure S1). Consequently, the core-ionization energies of Ti(OH)4 and Ti(OH)3NH2, E2p1/2 and E2p3/2 derived from the 2p1/2 and 2p3/2 states, respectively, were estimated asandwhere E2p is the Ti 2p core-ionization energy calculated using the VQD method.

It should be noted that calculations using cluster models provide values for gas phase molecules, while the experimental values for TiO2 and N-TiO2 have been obtained using solid materials. When determining the ionization energies of finite systems such as molecules, the vacuum level is usually chosen as the reference. In contrast, the ionization energies of solids and surfaces are typically reported relative to the Fermi levels.[75−77] Thus, the standard for the ionization energy will be different between molecules in the gas phase and solid materials. In addition, a core hole in a solid can be greatly stabilized by charge redistributions in the surrounding continuum, while such effects are absent in the gas phase.[78,79] Therefore, the core-ionization energies that were calculated based on the cluster models could not be directly compared to the experimental data obtained from solids. For this reason, we also adjusted the calculated values so that the simulated E2p value for Ti(OH)4 coincided with the experimental E2p value for bulk TiO2 as estimated from the E2p1/2 and E2p3/2 values. It is noteworthy that theoretical methods for the calculation of absolute core ionization energies in solids from first-principles have been investigated in detail.[76−80]

Results and Discussion

Calculations of Core-Excited States for CO, H2CO, and HCN

The weighting coefficients used in the VQD calculations and the calculated core-excitation energies for CO, H2CO, and HCN molecules are presented in Table with available experimental values.[81−84] The 1s (ε1s) and π* (επ*) orbital energies which were used for the estimation for weighting coefficients are summarized in Table S1. Because the ε1s energies increased in the order of C 1s > N 1s > O 1s, the weighting coefficients were in the reverse order of O 1s > N 1s > C 1s. The 1s → π* core-excited states could be successfully calculated using the VQD method with these weighting coefficients. The calculated core-excitation energies deviated from the experimental values by only 3–5% in the case that Hartree–Fock orbitals were employed. Remarkably, these deviations were further reduced to less than 1.1% by incorporating orbital relaxations. In particular, discrepancies of less than 0.4% were achieved using the cc-pVTZ or cc-pVQZ basis set, meaning that the core-excitation energies calculated using the VQD method were in quantitatively good agreement with the experimental data. The CAS contained only the minimum configurations describing the wave function of the core-excited state, and so the results were highly dependent on whether or not orbital optimizations were applied. The calculated core-excitation energies of H2CO were lower than those of CO as well as the experimental values. The calculated core-excitation energies were also correlated with the local chemical-bonding states, indicating that the VQD calculations could be applied to structural analyses based on comparisons with experimental XANES spectra.

Table 1

Weighting Coefficients in the Cost Function for the VQD Calculations: Overlap Weights (β), s2 Number Weights (w1), sz Number Weights (w2) and Particle Number Weights (w3), and Calculated 1s → π* Core Excitation Energies for CO, H2CO, and HCNa

		weighting coefficients (hartrees)				core excitation energy (eV)
molecule	core orbital	β	w1	w2	w3		basis set	without orbital optimization		with orbital optimization
CO	C 1s	12.0	21.4	48.0	12.0	calc	cc-pVDZ	298.39	(+3.82%)	290.47	(+1.07%)
							cc-pVTZ	298.79	(+3.96%)	288.39	(+0.34%)
							cc-pVQZ	299.50	(+4.21%)	288.10	(+0.24%)
						exptlb				287.4
	O 1s	21.0	37.4	84.0	21.0	calc	cc-pVDZ	552.92	(+3.52%)	536.14	(+0.38%)
							cc-pVTZ	553.25	(+3.59%)	533.97	(−0.03%)
							cc-pVQZ	553.74	(+3.68%)	533.57	(−0.10%)
						exptlc				534.1
H2CO	C 1s	12.0	21.4	48.0	12.0	calc	cc-pVDZ	297.99	(+4.34%)	288.54	(+1.03%)
							cc-pVTZ	298.20	(+4.41%)	286.58	(+0.34%)
							cc-pVQZ	298.73	(+4.60%)	286.37	(+0.27%)
						exptld				285.6
	O 1s	21.0	37.4	84.0	21.0	calc	cc-pVDZ	548.83	(+3.40%)	532.95	(+0.41%)
							cc-pVTZ	549.35	(+3.50%)	530.82	(+0.00%)
							cc-pVQZ	549.88	(+3.60%)	530.45	(−0.07%)
						exptld				530.8
HCN	C 1s	12.0	21.4	48.0	12.0	calc	cc-pVDZ	299.36	(+4.53%)	289.38	(+1.04%)
							cc-pVTZ	299.55	(+4.59%)	287.27	(+0.30%)
							cc-pVQZ	300.15	(+4.80%)	287.01	(+0.21%)
						exptle				286.4
	N 1s	16.0	28.5	64.0	16.0	calc	cc-pVDZ	415.33	(+3.91%)	402.53	(+0.71%)
							cc-pVTZ	415.78	(+4.02%)	400.28	(+0.15%)
							cc-pVQZ	416.49	(+4.20%)	399.97	(+0.07%)
						exptle				399.7

Available experimental values are also shown. The values in parentheses are deviations from the experimental values.

Reference (81).

Reference (82).

Reference (83).

Reference (84).

The calculated core-excitation energies were reduced depending on the basis set in the order of cc-pVDZ > cc-pVTZ > cc-pVQZ when orbital optimizations were carried out. A larger basis set was advantageous because it allowed orbital relaxation effects on core hole formation to be incorporated, thus reducing the energy of the core-excited state relative to that of the ground state. This obvious effect of the basis set size on the core-excitation energy was not observed in calculations performed without orbital optimization. Thus, the effects of enlarging the basis set were also fully incorporated in association with orbital optimization. Nevertheless, the calculated core-excitation energies deviated from the experimental values even when the larger cc-pVTZ and cc-pVQZ basis sets were used. The remaining discrepancies can be attributed to the lack of dynamical correlation within the present calculations.

Calculations of Core-Ionized States for CH4, NH3, H2O, and FH

The weighting coefficients used in the VQD calculations and the calculated core-ionization energies are collected in Table along with the available experimental data.[85] The ε1s energies for CH4, NH3, H2O, and FH molecules used for the estimation of the weighting coefficients are presented in Table S2. The weighting coefficients decreased in the order of FH > H2O > NH3 > CH4 because the ε1s energies were in the order of C 1s > N 1s > O 1s > F 1s. By adopting these weighting coefficients, the 1s core-ionized states were successfully calculated as |ψ(θ1)⟩. The resulting energies were found to deviate from the experimental values by 3.01–4.24% when orbital optimization was not applied, although these discrepancies were significantly reduced (to less than 1%) by carrying out orbital optimizations. The core-ionization energies obtained using the VQD method were in quantitatively good agreement with the experimental data in the case that orbital optimizations were employed, as was also the case for the core-excited state calculations. These results indicate that the peak position in an XPS spectrum can be predicted with reasonable accuracy using VQD calculations.

Table 2

Weighting Coefficients in the Cost Function for the VQD Calculations: Overlap Weights (β), s2 Number Weights (w1), sz Number Weights (w2), and Particle Number Weights (w3), and Calculated Core Ionization Energies for CH4, NH3, H2O, and FHa

	weighting coefficients (hartrees)				core ionization energy
molecule	β	w1	w2	w3		basis set	without orbital optimization		with orbital optimization
CH4	12.0	21.4	48.0	12.0	calc	cc-pVDZ	302.23	(+3.91%)	290.51	(−0.12%)
						cc-pVTZ	302.05	(+3.85%)	288.35	(−0.86%)
						cc-pVQZ	302.15	(+3.88%)	288.20	(−0.92%)
					exptlb				290.86
NH3	16.0	28.5	64.0	16.0	calc	cc-pVDZ	422.77	(+4.24%)	407.63	(+0.51%)
						cc-pVTZ	422.70	(+4.22%)	405.64	(+0.02%)
						cc-pVQZ	422.78	(+4.24%)	405.41	(−0.04%)
					exptlb				405.57
H2O	21.0	37.4	84.0	21.0	calc	cc-pVDZ	559.21	(+3.58%)	541.64	(+0.33%)
						cc-pVTZ	559.32	(+3.61%)	539.60	(−0.05%)
						cc-pVQZ	559.46	(+3.63%)	539.27	(−0.11%)
					exptlb				539.86
FH	27.0	48.0	108.0	27.0	calc	cc-pVDZ	715.05	(+3.01%)	695.84	(+0.24%)
						cc-pVTZ	715.28	(+3.04%)	693.63	(−0.08%)
						cc-pVQZ	715.41	(+3.06%)	693.20	(−0.14%)
					exptlb				694.18

Available experimental values are also shown. The values in parentheses are deviations from the experimental values.

Reference (86).

The calculated core-ionization energies decreased with increasing size of the adopted basis set when applying orbital optimization, in the order of cc-pVDZ > cc-pVTZ > cc-pVQZ. This ordering suggests that orbital relaxation, along with core-hole formation, could be more accurately described using a more flexible basis set, leading to a downward shift of the energy level of the core-ionized state relative to that of the neutral ground state.[44] The core-ionization energies calculated without orbital optimization were not clearly correlated with the size of the basis set because the effects of the orbital relaxations were only insufficiently considered. The core-ionization energies were found to be generally underestimated in the case that larger basis sets such as the cc-pVTZ or cc-pVQZ sets were employed. It has also been reported that the core-ionization energy is typically underestimated in the framework of the ΔSCF method.[86] Therefore, the deviations from the experimental values can be attributed to a lack of incorporated dynamical correlation. The charge-neutral ground state has one more electron than the core-ionized state, and so it is reasonable to expect that the extent of dynamical electron correlation in the charge-neutral ground state will be larger than that in the core-ionized state. Therefore, the gap between these states is underestimated when the dynamical correlation is not fully considered.

Effect of the Weighting Coefficients

The present calculations were carried out while varying the weighting coefficients, and Figure shows the primary configurations of the electronic states obtained for the H2O molecule. The calculated electronic state energies, eigenvalues, and imposed penalties are provided in Table . Table summarizes the calculated energies and the assignments of the first, second, and third electronic states: |ψ(θ0)⟩, |ψ(θ1)⟩, and |ψ(θ2)⟩ The costs of these states are summarized in Table . Below, we review the results with respect to each coefficient.

Figure 1

Primary configurations of the electronic states determined from VQD calculations of the core-ionized states for H2O. Here, ψ1 and ψ2 correspond to the O 1s orbital and HOMO, respectively.

Table 3

Calculated Energies, Eigenvalues of the Observables, and Imposed Penalty Terms of |ψ⟩ states (k ≥ 1) Obtained from VQD Calculations of the Core-Ionized States for H2Oa

		eigenvalue
electronic state	energy	S(S + 1)	MS	N	imposed penalty terms
|122̅⟩	–55.4760	3/4	1/2	3	0
|11̅2⟩	–75.5335	3/4	1/2	3	β
|11̅22̅⟩	–76.0268	0	0	4	(9/16)w1 + (1/4)w2 + w3
|11̅2̅⟩	–75.5335	3/4	–1/2	3	w2
|11̅⟩	–74.2798	0	0	2	(9/16)w1 + (1/4)w2 + w3

Each electronic state is denoted by a Slater determinant that corresponds to its main configuration as shown in Figure .

Table 4

Electronic State Energies (E(θ)) and Assignments of |ψ(θ)⟩ for H2O as Calculated Using the VQD Method (k = 0, 1, and 2)a

					|ψ(θ0)⟩		|ψ(θ1)⟩		|ψ(θ2)⟩
condition	β	w1	w2	w3	E(θ0)	assignment	E(θ1)	assignment	E(θ2)	assignment
(i)	15.0	37.4	84.0	21.0	–75.5335	|11̅2⟩	–75.5335	|11̅2⟩	–55.4760	|122̅⟩
	18.0	37.4	84.0	21.0	–75.5335	|11̅2⟩	–75.5335	|11̅2⟩	–55.4760	|122̅⟩
	21.0	37.4	84.0	21.0	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–75.5335	|11̅2⟩
	24.0	37.4	84.0	21.0	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–75.5335	|11̅2⟩
	27.0	37.4	84.0	21.0	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–75.5335	|11̅2⟩
(ii)	21.0	24.0	0.0	0.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–76.0268	|11̅22̅⟩
	21.0	32.0	0.0	0.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–76.0268	|11̅22̅⟩
	21.0	37.4	0.0	0.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–55.4760	|122̅⟩
	21.0	40.0	0.0	0.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–55.4760	|122̅⟩
	21.0	48.0	0.0	0.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–55.4760	|122̅⟩
(iii)	21.0	0.0	60.0	0.0	–75.5335	|11̅2⟩	–76.0268	|11̅22̅⟩	–74.2798	|11̅⟩
	21.0	0.0	72.0	0.0	–75.5335	|11̅2⟩	–76.0268	|11̅22̅⟩	–74.2798	|11̅⟩
	21.0	0.0	84.0	0.0	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–76.0268	|11̅22̅⟩
	21.0	0.0	96.0	0.0	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–75.5335	|11̅2⟩
	21.0	0.0	108.0	0.0	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–75.5335	|11̅2⟩
(iv)	21.0	0.0	0.0	15.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–76.0268	|11̅22̅⟩
	21.0	0.0	0.0	18.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–76.0268	|11̅22̅⟩
	21.0	0.0	0.0	21.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–55.4760	|122̅⟩
	21.0	0.0	0.0	24.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–55.4760	|122̅⟩
	21.0	0.0	0.0	27.0	–75.5335	|11̅2⟩	–75.5335	|11̅2̅⟩	–55.4760	|122̅⟩
(v)	21.0	0.0	21.0	15.75	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–76.0268	|11̅22̅⟩
	21.0	7.0	21.0	11.8125	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–76.0268	|11̅22̅⟩
	21.0	14.0	21.0	7.875	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–76.0268	|11̅22̅⟩
	21.0	21.0	21.0	3.9375	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–76.0268	|11̅22̅⟩
	21.0	28.0	21.0	0.0	–75.5335	|11̅2⟩	–55.4760	|122̅⟩	–76.0268	|11̅22̅⟩

All values are in hartrees.

Table 5

Calculated Electronic State Energies for |ψ(θ)⟩ (E(θ)), Cost Functions (L(θ)), and L(θ) – E(θ) Values (k = 1, 2)a

					|ψ(θ1)⟩			|ψ(θ2)⟩
condition	β	w1	w2	w3	E(θ1)	L(θ1)	L(θ1) – E(θ1)	E(θ2)	L(θ2)	L(θ2) – E(θ2)
(i)	15.0	37.4	84.0	21.0	–75.5335	–60.5335	15.0	–55.4760	–55.4760	0.0
	18.0	37.4	84.0	21.0	–75.5335	–57.5335	18.0	–55.4760	–55.4760	0.0
	21.0	37.4	84.0	21.0	–55.4760	–55.4760	0.0	–75.5335	–54.5335	21.0
	24.0	37.4	84.0	21.0	–55.4760	–55.4760	0.0	–75.5335	–51.5335	24.0
	27.0	37.4	84.0	21.0	–55.4760	–55.4760	0.0	–75.5335	–48.5335	27.0
(ii)	21.0	24.0	0.0	0.0	–75.5335	–75.5335	0.0	–76.0268	–62.5268	13.5
	21.0	32.0	0.0	0.0	–75.5335	–75.5335	0.0	–76.0268	–58.0268	18.0
	21.0	37.4	0.0	0.0	–75.5335	–75.5335	0.0	–55.4760	–55.4760	0.0
	21.0	40.0	0.0	0.0	–75.5335	–75.5335	0.0	–55.4760	–55.4760	0.0
	21.0	48.0	0.0	0.0	–75.5335	–75.5335	0.0	–55.4760	–55.4760	0.0
(iii)	21.0	0.0	60.0	0.0	–76.0268	–61.0268	15.0	–74.2798	–59.2798	15.0
	21.0	0.0	72.0	0.0	–76.0268	–58.0268	18.0	–74.2798	–56.2798	18.0
	21.0	0.0	84.0	0.0	–55.4760	–55.4760	0.0	–76.0268	–55.0268	21.0
	21.0	0.0	96.0	0.0	–55.4760	–55.4760	0.0	–75.5335	–54.5335	21.0
	21.0	0.0	108.0	0.0	–55.4760	–55.4760	0.0	–75.5335	–54.5335	21.0
(iv)	21.0	0.0	0.0	15.0	–75.5335	–75.5335	0.0	–76.0268	–61.0268	15.0
	21.0	0.0	0.0	18.0	–75.5335	–75.5335	0.0	–76.0268	–58.0268	18.0
	21.0	0.0	0.0	21.0	–75.5335	–75.5335	0.0	–55.4760	–55.4760	0.0
	21.0	0.0	0.0	24.0	–75.5335	–75.5335	0.0	–55.4760	–55.4760	0.0
	21.0	0.0	0.0	27.0	–75.5335	–75.5335	0.0	–55.4760	–55.4760	0.0
(v)	21.0	0.0	21.0	15.75	–55.4760	–55.4760	0.0	–76.0268	–55.0268	21.0
	21.0	7.0	21.0	11.8125	–55.4760	–55.4760	0.0	–76.0268	–55.0268	21.0
	21.0	14.0	21.0	7.875	–55.4760	–55.4760	0.0	–76.0268	–55.0268	21.0
	21.0	21.0	21.0	3.9375	–55.4760	–55.4760	0.0	–76.0268	–55.0268	21.0
	21.0	28.0	21.0	0.0	–55.4760	–55.4760	0.0	–76.0268	–55.0268	21.0

All values are in hartrees.

With regard to condition (i), the core-ionized state represented by |122̅⟩ was obtained as the ansatz state |ψ(θ1)⟩ with β ≥ 21.0. In contrast, the state was |ψ(θ2)⟩ and the ground state |11̅2⟩ was obtained as both |ψ(θ0)⟩ and |ψ(θ1)⟩ with β < 21.0. In the case of β < 21.0, |ψ(θ1)⟩ was nonorthogonal to |ψ(θ0)⟩ and the penalty of β was imposed on its cost function. However, ⟨11̅2|Ĥ|11̅2⟩ + β was still less than the cost function for the core-ionized state |122̅⟩. Therefore, the ground state |11̅2⟩ was obtained again as |ψ(θ1)⟩. If |ψ(θ2)⟩ is also |11̅2⟩, its cost function will be ⟨11̅2|Ĥ|11̅2⟩ + 2β because |ψ(θ2)⟩ is nonorthogonal to both |ψ(θ0)⟩ and |ψ(θ1)⟩, and so the cost function will be greater than that of the core-ionized state |122̅⟩. Thus, the core-ionized state |122̅⟩ was obtained as |ψ(θ2)⟩. In the case with β ≥ 21.0, the cost function for |11̅2⟩ became greater than that for |122̅⟩. Consequently, |122̅⟩ was obtained as |ψ(θ1)⟩ and |11̅2⟩ was obtained as |ψ(θ2)⟩ in ascending order of the cost function.

For condition (ii), the ansatz state |ψ(θ1)⟩ was calculated to be |11̅2̅⟩ in all cases. Because w2 = 0 in this case, no penalty was imposed on the deviation of the spin quantum number. In addition, the |11̅2̅⟩ state is orthogonal to |11̅2⟩, resulting in no penalty being imposed by β. As a result, the cost function for |11̅2̅⟩ was the same as that for the ground state |11̅2⟩. Therefore, |11̅2̅⟩ was consistently obtained as |ψ(θ1)⟩ and other states were obtained as |ψ(θ2)⟩. In the case of w1 < 37.4, the charge neutral ground state |11̅22̅⟩ with four electrons was obtained as |ψ(θ2)⟩. Since the charge neutral ground state |11̅22̅⟩ is orthogonal to the ground state of the ionized state |11̅2⟩ and w2 = w3 = 0, only an L of (9/16)w1 was imposed on its cost function. Therefore, ⟨11̅22̅|Ĥ|11̅22̅⟩ + (9/16)w1 was less than the energy of the core ionized state |122̅⟩ for w1 < 37.4. The ordering of the |11̅22̅⟩ and |122̅⟩ cost functions was inverted in the case of w1 ≥ 37.4 and |122̅⟩ was obtained as |ψ(θ2)⟩.

The calculation results for condition (iii) were complicated. The calculations with w2 < 84.0 provided the charge-neutral ground state |11̅22̅⟩ as |ψ(θ1)⟩ and the valence-double-ionized state |11̅⟩ as |ψ(θ2)⟩. Both states are orthogonal to the ground state |11̅2⟩ and w1 = w3 = 0, such that an Lpenalty of (1/4)w2 was imposed on |11̅22̅⟩. However, for w2 < 84.0, the Lpenalty value was less than the energy gap between |11̅2⟩ and |122̅⟩ and so |11̅22̅⟩ was obtained as |ψ(θ1)⟩. The same penalty was imposed on |11̅⟩ for w1 = w3 = 0, and was less than the energy gap between |11̅⟩ and |122̅⟩. Consequently, |11̅⟩ was obtained as |ψ(θ2)⟩. In the case with w2 = 84.0, the core-ionized state |122̅⟩ was obtained as |ψ(θ1)⟩ and |11̅22̅⟩ was obtained as |ψ(θ2)⟩ because L(θ1) = ⟨122̅|Ĥ|122̅⟩ was less than L(θ2) = ⟨11̅22̅|Ĥ|11̅22̅⟩ + (1/4)w2. In the case of w2 > 84.0, the ground state |11̅2⟩ was obtained as |ψ(θ2)⟩ because ⟨11̅22̅|Ĥ|11̅22̅⟩ + (1/4)w2 became greater than ⟨11̅2|Ĥ|11̅2⟩ + β.

The results obtained using condition (iv) were similar to those observed for condition (ii). The |ψ(θ1)⟩ state was calculated to be |11̅2̅⟩ in all cases and, although |Δ2| = 1, no penalty was imposed because w2 = 0. In the case of w3 < 21.0, the neutral ground state |11̅22̅⟩ was obtained as |ψ(θ2)⟩ because its Lpenalty of w3 was less than the energy gap between |11̅22̅⟩ and |122̅⟩. The core-ionized state |122̅⟩ was obtained as |ψ(θ2)⟩ for w3 ≥ 21.0, suggesting that the ordering of the cost function was inverted due to the increase in w3.

In conclusion, the weighting coefficients estimated from eqs and 25 were adequate to raise the cost of undesired ansatz states so that the target state was obtained. In addition, the penalty terms imposed on the ansatz states shown in Table suggest that fewer weighting coefficients were effective for this system. A penalty of w2 was imposed on |11̅2̅⟩, and so w2 = −ε1s was sufficient to raise the cost function of |11̅2̅⟩ above that of |122̅⟩. Accordingly, (9/16)w1 + w3 > – (3/4)ε1s was suitable for increasing the cost functions for |11̅22̅⟩ and |11̅⟩ such that they were higher than that for |122̅⟩. The calculations using the weighting coefficients meeting these requirements were carried out as condition (v). The core-ionized state |122̅⟩ and the neutral ground state |11̅22̅⟩ were obtained as |ψ(θ1)⟩ and |ψ(θ2)⟩, respectively. Thus, the lower limit of each weighting coefficient was dependent on the system. The weighting coefficients estimated from eqs and 25 were based on the minimum changes in the eigenvalues;[27] therefore, the estimated values were expected to be consistently equal to or larger than the lower limits that were required. The estimations based on eqs or 25 are potentially effective for any construction of the active space. However, the results in Table indicate that various electronic states other than the target state can be found in the calculations; the number of electronic states that can be found in the calculations rapidly increases with the active orbitals and electrons. Those states can be local minima and disturb the minimization of the cost function and its convergence. The influence of an expansion of the active space on the calculation results and the minimization process of the cost function should be clarified through further applications.

Calculation Results for the TiO2 and N-TiO2 Models

The optimized Cartesian coordinates for the models are summarized in Tables S3 and S4. The calculation results for the optimized structures are summarized in Table and discussed below.

Table 6

Calculated HOMO → LUMO Excitation Energies, O 1s → LUMO Core-Excitation Energies, Calculated Ti 2p Core-Ionization Energies (E2p), and Estimated E2p1/2 and E2p3/2 Values for the Ti(OH)4 and Ti(OH)3(NH2) Models; Experimental Values Obtained for Bulk TiO2 and N-TiO2 Are Also Presented

		HOMO → LUMO excitation energy			Ti 2p core-ionized (eV)
model	orbital	(eV)	(nm)	O 1s core-excitation energy (eV)	E2p	E2p1/2	E2p3/2
Ti(OH)4	not optimized	9.25	134.1	555.92	492.01	495.81	490.11
	optimized	3.21	386.3	532.41	477.42	481.22	475.52
Ti(OH)4 (shifted)					460.50a	464.30	458.60
TiO2 (exptl)		3.2b	387.4b	531.1c	460.5d	464.3,e 464.5f	458.6,e 458.7f
Ti(OH)3NH2	not optimized	7.67	161.7	556.36	491.06	494.86	489.16
	optimized	2.92	424.5	533.01	475.98	479.78	474.08
Ti(OH)3NH2 (shifted)					459.06g	462.86	457.16
N-TiO2 (exptl)		2.5b	495.9b	531.0–531.4h		464.0,i 464.1j	458.0,e 458.3j

Shifted downward by 16.92 eV from the calculated values after orbital optimizations so that the calculated E2p for Ti(OH)4 coincides with the experimental value of 460.5 eV.

Data obtained from ref (42) based on the absorption band edge.

Data taken from ref (87).

Estimated from the experimental E2p1/2 and E2p3/2 values in ref (73) utilizing eqs and (27.

Data taken from ref (73).

Data taken from ref (74).

Shifted downward by 16.92 eV from the calculated values for Ti(OH)3(NH2) after orbital optimizations.

Data taken from ref (87). The values were shifted upward with increasing amounts of nitrogen doping.

Extracted from Figure 8b of ref (73).

Extracted from Figure 2a of ref (74).

HOMO → LUMO Excited States

The HOMO and LUMO energy levels and weighting coefficients used in the calculations of HOMO → LUMO excited states are summarized in Table S5. The εHOMO–LUMO values were rounded up to the first decimal place prior to being used for the estimation of weighting coefficients. Since the HOMO → LUMO excitation energies were much smaller than the core-excitation energies and core-ionization energies, the weighting coefficients for the calculations of HOMO → LUMO excited states were also smaller.

The calculated excitation energies are presented in Table along with the wavelengths converted from the excitation energies; the total electronic state energies are shown in Table S6. The HOMO → LUMO excited states corresponding to the S1 states were obtained as the |ψ(θ1)⟩ ansatz states using the weighting coefficients in Table S5. The calculated absorption wavelength for Ti(OH)3(NH2) was longer than that for Ti(OH)4, indicating that the absorption wavelength for TiO2 should be red-shifted following nitrogen doping, in qualitatively good agreement with experimental findings.[42,43]Figure presents images of the HOMOs and LUMOs for the models as calculated using the Hartree–Fock method. The HOMO and LUMO for Ti(OH)4 were mainly derived from O 2p and Ti 3d orbitals, respectively. In contrast, the HOMO for Ti(OH)3(NH2) originated from the N 2p orbital. Because the energy level of the N 2p orbital is higher than that of the O 2p orbital (Table S5), the HOMO–LUMO gap is decreased by the substitution of an NH2 group for an OH group, which is analogous to the redshift of the absorption wavelength for bulk TiO2 as a result of nitrogen doping. The calculated HOMO → LUMO excitation energies were also greatly decreased following orbital optimization. The calculated energies for Ti(OH)4 and Ti(OH)3(NH2) were 3.21 and 2.92 eV, respectively, in good agreement with the experimentally observed band gaps of 3.2 eV for TiO2 and 2.5 eV for N-TiO2.[42] These results suggest that the adopted models had similar electronic properties to those of the bulk systems, even though they were much smaller in size. It should also be noted that the calculation results obtained using a cluster model depend on the properties of the cluster and that dynamical correlations were not considered in the present VQD calculations. The errors derived from these conditions might counteract one another to provide successful results.

Figure 2

Images showing the HOMO and LUMO for (a) Ti(OH)4 and (b) Ti(OH)3(NH2) as calculated using the Hartree–Fock method.

O 1s Core-Excited States

The O 1s and LUMO energy levels and weighting coefficients used in the calculations of the O 1s → LUMO excited states are summarized in Table S7. Note that the ε1s–LUMO value was rounded up to the first decimal place prior to being used for the estimation of the weighting coefficients. The calculated O 1s → LUMO excited state energies for the Ti(OH)4 and Ti(OH)3(NH2) models are provided in Table S8. These energies were averaged for the calculation of the core excitation energy value of each model, and the core excitation energy of Ti(OH)3(NH2) was found to be higher than that of Ti(OH)4 by approximately 0.6 eV (Table ). These results suggest that the peak in the XANES spectrum of TiO2 originating from the O 1s → LUMO excitation was shifted toward higher energy as a result of nitrogen doping, in agreement with experimental observations.[87] The core-excitation energies calculated with orbital optimization were generally in good agreement with the corresponding experimental values. The LUMOs shown in Figure suggest that these orbitals primarily originated from Ti 3d orbitals and that their energy levels were shifted upward due to the substitution of a N atom for an O atom of Ti(OH)4 (Table S8). In this oxide, the Ti–O bond is polarized to Tiδ+–Oδ− because of the high electronegativity of the oxygen atom.[88] Consequently, the electron density in the Ti atom is increased by the substitution of a N atom for an O atom because of the lower electronegativity of N compared with O. This effect, in turn, is responsible for the upward shift of the Ti 3d level.[74] As a result, the core excitation energy of Ti(OH)3(NH2) can be higher than that of Ti(OH)4.

The calculated core-excitation energies were decreased by approximately 13 eV following orbital optimization, suggesting that significant orbital relaxation was associated with core-hole formation compared with HOMO–LUMO excitation. Although a hole in the HOMO primarily affects the electrons excited to the virtual orbitals, a core hole additionally affects the valence occupied orbitals, leading to significant relaxation. As a result, the calculated O 1s → LUMO core-excitation energies were greatly reduced through orbital optimization.

The calculated shift of the core-excitation energy due to the substitution of a N atom for an O atom was larger than the experimentally observed value. However, the substitution of a N atom for one of the O atoms in the Ti(OH)4 model corresponded to a nitrogen doping level of 25 atom %, which was much larger than the concentration of less than 0.2 atom % applied in the experiments.[87] Therefore, the effects of nitrogen doping could be overestimated in the calculations results. A larger model capable of accurately simulating lower levels of doping would therefore be necessary for quantitatively accurate calculations.

Ti 2p Core-Ionized States

Ti 2p orbital energies calculated using the Hartree–Fock method and the weighting coefficients estimated from these orbital energies are summarized in Table S9, and the total electronic state energies are shown in Table S10. The three orbital energies were determined to be approximately −18.08 hartree for Ti(OH)4 and −18.05 hartree for Ti(OH)3(NH2). Note that all of the −ε2p values, when rounded up to the first decimal place, were 18.1 hartree. As a result, common weighting coefficients were adopted for both models. The calculation results for the Ti 2p core-ionized states are provided in Table . The calculated core-ionization energies derived from the three Ti 2p orbitals were averaged to obtain the E2p value for each model and the resulting E2p1/2 and E2p3/2 values for Ti(OH)3(NH2) were found to be lower than those for Ti(OH)4. The results suggest that the Ti 2p peaks in the XPS spectrum should be shifted to lower energy due to nitrogen doping, in agreement with the experimental findings.[73,74]

Since the Ti 2p level was shifted upward as a consequence of the substitution of a N atom for an O atom (Table S9), the Ti 2p core-ionization energy was decreased. However, the calculated core-ionized energies of Ti(OH)3(NH2) were lower than those of Ti(OH)4 by 0.9–1.5 eV, and this shift was larger than the experimental values of 0.3–0.6 eV.[73,74] This discrepancy can be attributed to the overestimated extent of nitrogen doping in the models, as with the discrepancy for the calculations of O 1s → LUMO core-excited states. The calculated core-ionization energies were decreased by approximately 15 eV following orbital optimization, suggesting that orbital relaxation had a significant effect.

Overall, the results of the VQD calculations suggest that nitrogen doping was responsible for the red-shift in the absorption band edge, the increase in the O 1s → Ti 3d excitation energy and the upward shift of the Ti 2p core-level, all of which were in qualitatively good agreement with the experimental data.

Conclusions

The VQD method has attracted attention as a promising quantum-classical hybrid algorithm for the calculation of excited states. In the present study, calculations of the core-excited states and core-ionized states for common molecules utilizing the VQD method were simulated on a classical computer, focusing on the penalty terms in the cost function. The weighting coefficients for these penalty terms were estimated on the basis of molecular orbital levels and the minimum deviations of the eigenvalues of the wave function regarding the spin and the number of electrons. Adopting this simple procedure allowed the core-level states to be successfully calculated. The relationship between the weighting coefficients and the resulting ansatz states was systematically examined on the basis of calculations of the O 1s core-ionized state for a water molecule. The results indicate that the weighting coefficients estimated by our procedure allowed the target states to be obtained by raising the cost of undesired states. The O 1s core-excited states and the Ti 2p core-ionized states for TiO2 and N-TiO2 photocatalysts serving as model compounds could also be calculated in the same manner as the valence excited states. The results of such calculations demonstrated that nitrogen doping induced a red-shift in the absorption band edge, increased the O 1s → Ti 3d excitation energy, and raised the Ti 2p core-level. These results were consistent with experimental findings, suggesting that VQD calculations can be applied to the analysis of functional materials in collaboration with experimental work. The results of this study are expected to provide valuable guidelines for future applications of the VQD method using quantum computers and also to inspire further development of the algorithm.