Robust ΔSCF calculations with direct energy functional minimization methods and STEP for molecules and materials

Electronically excited states play an important role in various photochemical processes, such as cycloaddition reactions,¹ excited-state dynamics,² and electron transfer processes.^3,4 The excited-state phenomena are used to investigate a variety of physical and chemical properties of, for example, oxide materials,^5,6 organic dye molecules,^7–10 as well as metal-surface charge transfer^11,12 and charge transfer in condensed phase systems.¹³ Calculations of vertical excitation energies and related properties are instrumental in characterizing the electronically excited states. Δ self-consistent field (ΔSCF)^14–17 is a promising method to compute the vertical excitation energy of various systems ranging from molecules to materials, including condensed phase systems.

The ΔSCF calculations are based on a non-Aufbau electronic configuration of the molecular orbitals (MOs) during the variational SCF calculations. To avoid the variational collapse of the ΔSCF calculation to the ground state, ΔSCF algorithms have been developed, such as maximum overlap methods (MOM),^18,19 which are based on the maximum overlap of new and old molecular orbitals during the SCF iterations. Other approaches have been developed to improve convergence issues for difficult cases of SCF convergence for excited states, for instance, the direct energy-targeting method (σ-SCF) by Ye et al.,²⁰ direct optimization of orbitals by Levi et al.,^21,22 the square gradient minimization (SGM) method by Hait and Head-Gordon,²³ or the state-targeted energy projection (STEP) approach by Carter-Fenk and Herbert.²⁴ 

In general, linear response time-dependent density functional theory (TDDFT) is often used for excitation energy calculations due to its balance of computational cost vs accuracy for molecules.^25–27 Alternatively, real-time TDDFT has also been applied for the calculation of vertical excitation energies.^28–31 However, TDDFT has shortcomings for charge transfer states, especially without the use of long-range corrected functionals.^32–35 Additionally, TDDFT does not account for double excitation transitions. For a more comprehensive outlook on ΔSCF and its importance, we refer to a recent paper by Hait and Head-Gordon.,³⁶ especially for cases such as charge transfer states, core-level excitations, and double excitation calculations.

The application of the ΔSCF method has grown beyond obtaining vertical excitation energies. For example, it has been applied to the potential energy surface calculation of excited molecules on surfaces^14,16,37 and nonadiabatic molecular dynamics.^38–43 The ΔSCF calculations are also well suited for the excited-state properties, such as transition dipole moments.⁴⁴ In general, the excited state MOs in ΔSCF calculations are not orthogonal to the ground state MOs, but an orthogonalization scheme can be employed to obtain the transition property.⁴⁴ In terms of implementation, the ΔSCF algorithms can easily be adapted to an existing ground-state SCF code. For example, the excited-state gradients are available at the cost of ground state gradients. The accuracy and convergence performance of ΔSCF methods have been compared for various kindS of molecular systems in Refs. 21–24 and 45. However, a major drawback of ΔSCF methods is the failure in achieving SCF convergence or variational collapse for difficult systems. Levi et al.^21,22 implemented a direct energy optimization method in combination with the MOM methods, which improves the excited-state SCF convergence.

The diagonalization-based SCF works well for most cases, but the SCF convergence might not be achieved in a difficult case. Direct minimization schemes for the energy calculations are more attractive due to their robustness since the direct minimization guarantees convergence if the electronic energy is reduced in each step of SCF iterations.^46–49 The direct minimization methods are based on minimizing the energy functional in each step.^50–52 Exponential parameterization of the orbitals is done by introducing additional rotation parameters. The minimization procedure requires electronic energy and energy gradient calculations. The orbital transformation (OT) method is available, such as steepest descent, conjugate gradient (CG), direct inversion of the iterative subspace (DIIS), and Broyden’s type 2 method^{46,47,53–55} together with golden line search to make the SCF robust in CP2K.^46,56 The computational cost of the OT method for M basis functions and N occupied MOs scales as O(MN²) with a prefactor for large systems,⁴⁶ whereas for the standard diagonalization SCF, the scaling is O(M³). The performance and robustness have been demonstrated for the SCF of molecules in vacuum and liquids for systems up to a few thousand atoms in Refs. 47 and 48.

In this work, we applied the direct minimization method using the OT^46,47,50,56 scheme to ΔSCF calculations in a modified version of the CP2K program package. In particular, the OT method for non-uniform occupation is applied in the calculations. We have adapted the code to include the fractional contribution of α and β MOs’ occupation for singlet and triplet excited states for ΔSCF calculations. We applied the OT based ΔSCF calculations to molecules, liquids obeying periodic boundary conditions (PBCs), and a periodic surface system to demonstrate the direct minimization OT ΔSCF convergence robustness. Additionally, we implemented another ΔSCF method, the STEP approach, which has been recently proposed by Carter-Fenk and Herbert.²⁴ The STEP implementation has been done for diagonalization based SCF calculation. Section IV shows both OT- and STEP-ΔSCF results for the two lowest singlet excitations of various heteroaromatic systems.

This paper is organized as follows: In Sec. II, we discuss the methodology of energy minimization for SCF calculation using OT schemes in CP2K. The STEP ΔSCF method using the diagonalization SCF method is also discussed in Sec. II. In Sec. III, the computational details are given for the calculations. The vertical excitation energy results of various heteroaromatic gas-phase molecules and condensed-phase systems, namely, ethylene solvated in water, uracil solvated in water, and a Re–phosphate complex attached to the anatase (TiO₂) surface, using ΔSCF are discussed in Sec. IV.

For the standard ground-state SCF calculation, the energy minimization is performed using the OT scheme in the QUICKSTEP program⁵⁶ of the CP2K program package. The QUICKSTEP program is an implementation of a hybrid Gaussian and plane wave (GPW) method.^56,57 In the GPW method, two representations of the electron density can be employed, allowing an efficient treatment of the electrostatic interaction for large or periodic system. The first representation is described in atom-centered contracted Gaussian functions, and the second representation of the electronic density is in plane waves. The interconversion between the Gaussian-function based one and the plane wave density is done using the mapping of a product of Gaussians on the real space grid and fast Fourier transforms (FFTs).⁵⁶ The PBCs are taken care of by FFT based treatment of the Poisson equation.

The MOs {ϕ_i} in the atomic-orbital (AO) basis {χ_μ} can be written as

where C_μi is the i-th MO coefficient and M is the number of AO basis functions. The electron density is defined as

and P is the one-electron density matrix

Here, f_i is the occupation number of the ith MO. In the case of restricted calculation, f_i = 2 for occupied orbitals and N is the number of MOs. For all the MOs, occupied and virtual, the exponential orbital transformation is given, in general, as $C′=C0′⁡exp(A)$⁠, where

Here, C′, $C0′$⁠, and A are M × M matrices. $C0′$ is any initial guess MO coefficient matrix. For uniformly occupied MOs, the occupied–occupied (oo) block of A can be set to 0, since the electronic energy E[C] is invariant with respect to rotation in oo subspace. However, in the case of fractional occupation, the orbitals are not uniformly occupied anymore and E[C] is not invariant with respect to the rotation of MO coefficients in the oo subspace. Thus, the occupied–occupied block of the transformation matrix is non-zero. The energy remains invariant with respect to orbital rotation in the virtual–virtual (vv) block for both uniform and non-uniform occupation of MOs and can be set to 0 without loss of generality,

The direct energy minimization requires the minimization of the energy functional E[C] with respect to MO coefficients C subject to the orthogonality constraint of the MOs,^46,48

where S is the overlap matrix with dimension M × M and I′ is the M × M identity matrix.

For the set of N uniformly occupied MOs C₀, the occupied-virtual block of A is parameterized using variable X for minimization of the energy functional E[C], where the variable X is linearly constraint as X^TSC₀ = 0.^46,58 A detailed derivation of generalized parameterization of exp(A) is given in Ref. 58. The parameterization of matrix C (with dimension M × N) for uniformly occupied MOs is given as⁴⁶ 

where the unitary matrix U is defined as

X is an M × N matrix, U is an N × N matrix, and X is linearly constrained as X^TSC₀ = 0. In the case of the non-uniform occupation of the MOs, the rotation of the oo block shall be considered due to non-invariance of the total energy. Therefore, the parameterization of A_oo block is included in the orbital transformation as

where (Y = −Y^†). exp(Y) is additional parameterization for the orbital rotation of the oo block, A_oo, and has a dimension of N × N. exp(Y) can be diagonalized using the spectral theorem for matrices as QVQ^†, where V is a complex diagonal matrix with elements V_ii = exp(−ιv_i), where v_i are the eigenvalues corresponding to ith column of the unitary matrix Q, and columns of matrix Q are the eigenvectors of ιY. Here, † in superscript represents the complex transpose.

As noted, in the case of fractional occupation, there are two subspaces (A_oo and A_ov) where the energy is not invariant with respect to orbital rotation; thus, the energy functional needs to be minimized with respect to both subspaces. For the ΔSCF calculations, the energy gradient calculation is done in two parts: one for A_ov and another for A_oo. The energy gradient for A_ov subspace, $∂E[C(X)]∂X$⁠, is the standard ground state energy formalism for uniformly occupied MOs given in terms of constrained gradient^46,58 and unconstrained gradient $(∂E∂C)$⁠.^46,58 The energy gradient term is given by⁴⁶ 

where

In Eqs. (10) and (11), HC represents $∂E∂C$⁠, H is the Kohn–Sham matrix expressed in AO basis functions, and R contains the eigenvectors of X^TSX. In Eq. (11), the D¹ and D² matrices are functions of eigenvalues of R and ⊗ is a matrix product by elements. The energy gradient term for additional parameterization of occ–occ subspace [exp(Y)] is given as

where

and D⁽³⁾ has elements

For OT ΔSCF calculation, the additional gradient term is evaluated together with standard energy gradient in the same loop during SCF iteration.

For unrestricted OT ΔSCF, the general exponential orbital transformation $(C′=C0′⁡exp(A))$ formulation is done analogously for each spin channel, the density, density matrix, and orbital transformation,

where $fiσ$ is the occupation number of the ith σ MO. During the SCF procedure, the SCF iteration loops over each spin channel and the total gradient is a sum of α and β contributions.

For the ΔSCF calculations, additional orbitals are added to include the fractional MOs by extending the C₀ matrix. The CG and DIIS minimization methods are available when the additional orbitals are included for the OT ΔSCF calculation. Due to the parameterization variable in the case of non-uniform occupied orbitals, the gradient term $∂E∂Y$ is added to the standard energy gradient [Eq. (10)]. For unrestricted ΔSCF, the occupation number f_i values are taken as 1 for completely occupied orbitals and between 0 and 1 for partially occupied orbitals. The values of fractional occupation are chosen based on desired electronic excitation in the ΔSCF calculations. For example, in electronic excitations involving the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO), the HOMO and LUMO orbitals’ occupation are provided between 0 and 1. The fractional occupation also allows us to model a single excitation, which is a combination of HOMO − 1 → LUMO and HOMO → LUMO + 1. Currently, only a limited number of excited states can be calculated using the available method since we cannot perform an excited state calculation with holes below HOMO. The implementation is available for both restricted and unrestricted calculations.

The STEP approach has been recently proposed by Carter-Fenk and Herbert.²⁴ In this approach, a level-shift constraint is applied on the Fock matrix by the projection operator,

on the virtual space, where {ϕ_a} are the virtual MOs. In the AO basis, the projector operator [Eq. (18)] can be written as

where the matrix V is the virtual block of the density matrix P. By adding the projector operator to the Fock matrix F, the modified Fock matrix becomes

Here, η is a shift parameter and SVS is the level-shift operator. The role of ηSVS is to increase the electronic orbital energy of involved virtual orbitals in constructing the projector operator [Eq. (18)] by η amount. The η parameter is defined as

Here, ɛ′ is an empirical parameter, ɛ_HOMO is the HOMO energy and ɛ_LUMO is the LUMO energy. As discussed in Ref. 24, the energy shift parameter given in Eq. (21) is a preferred choice for the balance of the number of iterations required for convergence of SCF. ɛ′ is added to avoid the HOMO and LUMO energy degeneracy, and in our calculations, it is taken as 0.1 Hartree. Our implementation of the STEP approach in CP2K is for the diagonalization based SCF.

All implementations, ΔSCF calculations, and density functional theory (DFT) based molecular dynamics (MD) calculations are done at the Kohn–Sham (KS) DFT level using the GPW method as implemented in CP2K.⁵⁶ The OT ΔSCF calculations are performed using the CG method available in QUICKSTEP.⁵⁶ The generalized gradient approximation (GGA) Perdew–Burke–Ernzerhof (PBE)^56,59 functional has been employed for the calculations. The Goedecker–Teter–Hutter (GTH) pseudopotentials^60–62 are used in all calculations. For the molecules in vacuum, the aug-TZVP-GTH basis set is used. For the solvated ethylene and uracil system, the TZVP-GTH basis set is used, and for the Re–phosphate complex system, the DZVP-MOLOPT-SR-GTH basis set. The STEP ΔSCF calculations for the heteroaromatic systems are done using diagonalization based SCF. The ATOMIC guess density is used as the initial SCF guess for both OT and STEP ΔSCF. The linear response time-dependent density functional perturbation theory (TDDFPT)⁶³ calculations for the heteroaromatic systems are done using CP2K. The nuclear structures of the heteroaromatic systems in Sec. IV are taken from the Ref. 64. The cutoff value for the SCF convergence, which is square root of the largest element of energy gradient vectors, is taken at 1.0 × 10⁻⁶ a.u. The excitation energy is computed by taking the difference of ground-state SCF energy and excited-state SCF energy for both OT ΔSCF and STEP ΔSCF methods.

The low-lying ππ^⋆ singlet excited states of ten heteroaromatic systems, shown in Fig. 1, are calculated using the OT ΔSCF and STEP ΔSCF methods. These singlet states are labeled L_a and L_b, where L_a corresponds to the HOMO → LUMO transition and L_b corresponds to nearly equal contributions of HOMO − 1 → LUMO and HOMO → LUMO + 1 transitions, and have been previously investigated at the TDDFT level.⁶⁴ The heteroaromatic systems are of interest in organic electronics and considered to be challenging test cases for excited-state electronic structure methods due to bi-reference character of ππ^⋆ states; see, for example, Refs. 64–67. The natural transition orbitals of the TDDFT calculations show that the second singlet excited state, L_b, is of a bi-reference character in these molecules.⁶⁴ 

During the OT ΔSCF calculations, the mixed contribution of α and β orbitals are considered using fractional smearing. In the fractional smearing approach, a list of fractional occupation numbers of α and β orbitals is employed. For instance, a singlet electronic transition is performed by changing both α and β orbitals occupations at equal footing in the list during unrestricted SCF calculations. For the OT ΔSCF L_a state calculation, the occupation number of the HOMO orbital for both spin channels is changed from 1 to 0.5. The occupation number of the LUMO orbital becomes 0.5 for both spin channels. Similarly, for the OT ΔSCF L_b state calculation, the occupation number of the HOMO − 1 and HOMO orbitals for both spin channels is changed from 1 to 0.75. The occupation number of the LUMO and LUMO + 1 orbitals becomes 0.25 for both spin channels. The smearing of all molecules is reported in supplementary material.

In our diagonalization based STEP ΔSCF calculations, only a single reference transition is included. The orbital occupation change is restricted to an integer. STEP ΔSCF calculations for the L_a state are based on a HOMO to LUMO orbital transition. In this case, the HOMO occupation changes from 1 to 0, and the LUMO occupation is changed from 0 to 1 for the same spin. For the L_b state, HOMO −1 to LUMO single orbital transition is considered. To calculate the corresponding triplet states using STEP ΔSCF, the occupation of α occupied orbital changes from 1 to 0, and then the occupation of β unoccupied orbitals becomes 1 or vice versa.

The OT ΔSCF and STEP ΔSCF results are reported in Table I. The excitation energy results for the L_a and L_b states of the compounds I–X (Fig. 1) using TDDFPT are also given in Table I. Based on the integrated spin density values in the STEP ΔSCF, we observe that the singlet excitation energies results have spin-contamination. Here, the spin energy purification^68,69 has been applied to the STEP ΔSCF calculations for the vertical excitation energy,

E_S is the purified singlet excited state energy, E_↑↓ is the energy of ΔSCF calculation for the singlet transition, and E_↑↑ is the energy of ΔSCF calculation for the triplet transition for the same occupied and unoccupied orbitals. Comparing the vertical excitation energies of the L_a and L_b states using STEP ΔSCF and OT ΔSCF, they vary by 0.02–0.50 eV for the compounds. The fractional smearing approach in OT ΔSCF calculation allows us to change the occupation of α and β MOs simultaneously for a singlet excited state. Thus, the singlet excited state energies obtained from the OT ΔSCF calculations do not require spin energy purification.^68,69

Vertical excitation energy results from ΔSCF differ from the TDDFPT excitation energies. The electronic energies of L_a states are systematically lower using the OT ΔSCF methods when compared to the TDDFPT results, ranging from 0.05 to 0.95 eV for compounds I–X. For compound V, the OT ΔSCF method gives a 0.95 eV lower energy than the TDDFPT result, whereas it is 0.05 eV for compound II. While comparing the vertical excitation energies of L_a using STEP ΔSCF with TDDFPT, the energies are lowered by 0.06–0.46 eV. Calculations with correlated methods, such as ADC(2) (algebraic diagrammatic construction up to second-order) and configuration interaction (CI) methods, as well as the TDDFT at PBE/aug-cc-pVTZ level are reported in Ref. 64, the reference excitation energies of which have been added to Table I. The CC2/aug-cc-pVTZ results⁶⁴ show that TDDFPT and ΔSCF excitation energies are systematically lower for both single-reference L_a state and bi-reference L_b states. In addition, L_a and L_b states’ energies are differently ordered compared to TDDFPT and ΔSCF results for various systems, such as compounds I, II, III, IX, and X.

Convergence: The OT ΔSCF calculations converge for all studied heteroaromatic molecules. However, the second singlet L_b state of compounds V, VII, and X did not converge using the STEP ΔSCF. While comparing the robustness of the direct energy minimization method with the diagonalization based ΔSCF, it is clear that direct minimization is advantageous for the L_b state calculations. There were also convergence issues for the lowest L_a singlet state using default SCF parameters in diagonalization based STEP ΔSCF. We recognized that the extent of direct mixing of new and old density matrices plays an essential role in the convergence of the diagonalization based ΔSCF. In some SCF calculations, the default value of 0.4 for the Pulay mixing (DIIS)⁵⁴ in CP2K did not converge, whereas lowering the mixing extent from 0.4 to 0.2 improved the SCF convergence. Thus, changing the default parameters for the ΔSCF calculations can be advantageous to improve convergence cases.

To demonstrate the applicability of the OT ΔSCF for charge transfer, hydrogen chloride (HCl) molecule’s first singlet excited state has been calculated. The singlet electronic charge transfer excitation of HCl has been investigated in Ref. 70. The HOMO to LUMO transition gives the excitation energy of 8.13 eV using OT ΔSCF. In another study,⁷⁰ the reported value of the HCl charge-transfer singlet state at the CC2/cc-pVTZ level is 8.23 eV, whereas TDDFT gives 7.55 and 7.79 eV using the PBE and CAM-B3LYP functionals, respectively, and cc-pVTZ basis set. The OT ΔSCF calculation for HCl converges smoothly. We also investigated the excitation energy profile of HOMO to LUMO singlet transition using OT ΔSCF and TDDFPT. The excitation energy value decreases when the HCl bond is elongated compared to the equilibrium bond length and increases when shortening the bond length. As shown in Fig. 2, the trend remains the same for OT ΔSCF and TDDFPT. Based on the excitation amplitude, the TDDFPT excitation results are dominated by a single HOMO to LUMO transition for all bond lengths. The convergence of OT ΔSCF results remains robust for all investigated bond distances of HCl.

We investigated the trans-azobenzene (Fig. 3) ground state energy and first singlet excitation energy profile (HOMO to LUMO transition) along the C-N-N-C (CNNC) dihedral angle. Azobenzene system has been of interest in many previous studies due to its photochemical properties, such as photoisomerization of cis-, trans-azobenzene.^71–76 The constraint geometry optimizations have been done for CNNC dihedral starting from trans geometry (CNNC angle of 180°) up to 20° at the PBE/aug-TZVP-GTH level. Figure 4 shows the OT ΔSCF and TDDFPT profile for the first singlet excited state. The relative ground state energy of each dihedral angle (θ) geometry ((E)_CNNC=θ) is calculated with respect to the highest energy point $((E)CNNC=80°)$ of the different CNNC angle investigated, which is at a dihedral angle of 80°,

The conical interaction region is around the energy point at 90°. Although a single reference method cannot describe well the excited electronic character near the conical interaction region, the OT ΔSCF converges smoothly and produces a well-defined topology near the conical region. However, the relative energy at 90° dihedral angle is below the ground state by 0.15 eV. The conical topology of TDDFPT is different compared to OT ΔSCF, since the relative energy of the studied singlet excited state is lower at three points, namely, 70, 80 and 90°, by 0.42, 0.09, and 0.38 eV respectively. In another study of the azobenzene’s excitation energy profile,⁷⁶ the XCDFT/B3LYP/TZP(Slater-type) results show the crossing point at 90° and have similar trend near the conical topology as in OT ΔSCF [Fig. 4(a)].

The vertical excitation energy of trans-azobenzene at the OT ΔSCF/PBE/aug-TZVP-GTH level is 2.353 eV, whereas spin-polarized TDDFPT/PBE/aug-TZVP-GTH gives a much lower value of 1.596 eV. In previous studies, the vertical excitation energy of trans-azobenzene has been reported at different levels of theory, such as CC2/aug-cc-pVTZ⁷⁴ (2.84 eV) and CAS–SCF/cc-pVTZ⁷⁵ (3.10 eV). The experimental vertical excitation energy in the gas phase is reported⁷⁷ to be 2.82 eV.

Several works about ethylene have been published in gas-phase environment; see, for example, Refs. 78–82. In addition, nonadiabatic molecular dynamics studies have been done for the ethylene gas-phase system using the ΔSCF methods.^38,40 Here, we focus on the ethylene’s vertical excitation energy in a condensed phase environment with PBC using OT ΔSCF to test the convergence robustness. An ethylene molecule was solvated with 27 water molecules in a periodic cubic box, shown in Fig. 5. The box size was first equilibrated using DFT-MD to a pressure of 1 atm in the isothermal–isobaric (NPT) ensemble at 300 K, followed by a more extended (14 ps) canonical (NVT) ensemble dynamics to thermalize the system to 300 K. For the solvated ethylene, a side box length of 9.85 Å is sufficient to contain the complete first solvation shell around the molecule. The displayed structure in Fig. 5 is randomly sampled from the last four ps of the NVT trajectory. The dynamics simulations were performed at the PBE/TZVP-GTH level of theory. The OT ΔSCF calculation with smeared occupation set up converges with the default SCF parameters in CP2K, whereas the diagonalization based ΔSCF is not able to converge.

The singlet vertical excitation energy for the HOMO to LUMO transition from the randomly sampled single configuration is 5.338 eV. The triplet excitation energy for HOMO to LUMO transition is 4.049 eV. While analyzing the electronic density difference between the excited states and the ground state, it becomes apparent that the excitations are primarily localized at the ethylene molecule and some near water molecules; see Fig. 1 in the supplementary material. For the gas-phase ethylene, the lowest singlet excitation energy is 5.891 eV and the corresponding triplet energy is 4.247 eV. The shift in the singlet excited state is due to the solvent effects and is an example of solvatochromism.⁸³ 

There have been several benchmark vertical excitation energy calculations of uracil systems.^85–87 Uracil is a nucleobase and a planar system in the gas-phase.^86,88 Here, similar to the ethylene condensed phase ΔSCF calculation, the uracil molecule has been investigated in solution. The solvated system is prepared with 126 water molecules in a periodic cubic box, with a box length of 15.40 Å. The structure shown in Fig. 6 is a randomly sampled configuration from a 30 ps long NVT trajectory. The OT ΔSCF calculation with smearing converges without any issue for the investigated states.

For this snapshot of solvated uracil, the singlet vertical excitation energy for HOMO to LUMO transition is 3.842 eV. The electronic density difference between the singlet excited state and the ground state shows that the excitation is localized at the uracil molecule; see Fig. 3 in the supplementary material. The HOMO to LUMO singlet transition in gas-phase uracil molecule is 4.128 eV. The difference between the HOMO to LUMO singlet excitation of uracil in the gas phase and the solvated uracil is 0.286 eV. The excitation energy shift of uracil in the solvent for the sampled structure can be mainly attributed to polarization due to the solvent water molecules and change in geometry of uracil.⁸⁹ The triplet excitation energy for the HOMO to LUMO transition in gas-phase uracil is 3.378 eV, whereas the triplet excitation energy for HOMO to LUMO transition in the solvent is 3.198 eV.

There has been an increasing interest in studying TiO₂ surface-based photocatalysis using computational chemistry methods,^90,91 in particular, the adsorption of organic dye or photosensitizing molecules on anatase or rutile.^92–94 Here, we demonstrate the applicability of OT ΔSCF by calculating singlet and triplet vertical excitation energies for a HOMO → LUMO transition.

In the model molecule–(anatase)surface system, the anatase slab has three layers, and the phosphate group’s oxygen atoms of the molecules are attached to the surface as depicted in Fig. 7. One of the phosphate groups is attached bidentate, and the other is monodentate through the oxygen atoms. The vertical excitation energy for the Re–phosphate complex on the anatase TiO₂ surface system is calculated at the PBE/DZVP-MOLOPT-SR level. These calculations are done using PBC and a cell size of 20.46 × 22.68 × 45.00 ų.

In Table II, the singlet vertical excitation energy is reported for the HOMO → LUMO transition. By looking at the electronic density difference plot, it is clear that the singlet excitation is localized near the molecular complex’s metal center; see Fig. 5 in the supplementary material. The triplet state for the HOMO → LUMO transition is at 1.576 eV. The OT ΔSCF calculations for the singlet and triplet calculation converge smoothly. The results show that the OT ΔSCF method can efficiently be applied to a molecule-material system for excited-state properties’ calculations.

We have demonstrated the use of OT for the ΔSCF method and investigated its behavior as well as the one of STEP ΔSCF as implemented in an extended version of the CP2K package. The robustness of the direct energy minimization approach is evident since the OT ΔSCF with smearing scheme converges for all the considered transitions of a set of ten heteroaromatic molecules, such as HOMO to LUMO, and transitions of bi-reference characters, such as a mixed contribution of HOMO to LUMO and HOMO − 1 to LUMO transitions. The ΔSCF vertical excitation energy results are mostly comparable with TDDFPT results although they are systematically shifted from the TDDFPT results for both L_a and L_b states.

The excitation energy profile of the HCl molecule’s HOMO to LUMO singlet transition is calculated using OT ΔSCF. The OT ΔSCF results are comparable with the TDDFPT results for all investigated bond distances of HCl for this charge transfer transition. In addition, trans-azobenzene’s excitation energy profile is computed to show further the robustness of the OT ΔSCF method for a more challenging case. The OT ΔSCF results near the conical topology of azobenzene’s singlet excited state are well defined, and it converges smoothly in the conical region.

We investigated periodic condensed phase systems’ excited states using the direct minimization OT ΔSCF method with smeared MO occupation, confirming the robustness of the SCF convergence. Moreover, by applying the OT ΔSCF to both liquids and molecule–metal surface system, we demonstrated the ΔSCF application to large periodic systems. Earlier, diagonalization based SCF was not able to converge for these structures. In the future, the OT ΔSCF method might efficiently be applied to perform non-adiabatic molecular dynamics for complex systems since it alleviates ΔSCF’s failure to converge for certain snapshot geometries obtained during the dynamics.

The supplementary material includes the plots of electronic density difference between the HOMO → LUMO singlet state and the ground state for ethylene, uracil, and the Re–phosphate complex; the fractional occupation list of MOs for the ten aromatic compounds, ethylene, uracil, and the Re–phosphate complex; and the coordinates of solvated ethylene, solvated uracil, and the Re–phosphate complex.

We acknowledge the Swiss National Science Foundation (Grant No. PP00P2 170667) and the University of Zurich (UZH) for funding support. We thank Dr. Momir Mališ (UZH) for discussion and providing the solvated systems’ structures and Eva Vandaele (UZH) for creating the structure of the model Re–phosphate complex on the anatase system. We also thank the Swiss National Supercomputing Center (CSCS) for computing resources (accounts: s1001 and s875).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.