The ΔSCF method for non-adiabatic dynamics of systems in the liquid phase

The interaction of light with matter is a fundamental process enabling life on Earth. Nature has developed complex pathways to harvest solar energy, perform a plethora of photochemical processes, or protect its compounds against photoradiation damage. Scientists are learning from nature every day and at the same time trying to extend the mechanisms to a vast range of new applications. Light-emitting diodes (LEDs) are a good example of a successful commercial application based on photofunctional materials.¹ The intense electroluminescence of many organic molecules and dyes subsequently inspired the development of organic LEDs. Modern multilayer organic LEDs combine a high efficiency and flexible design with a fast response time to the applied electric field as compared to liquid crystal displays.² Nevertheless, intense investigation into the rational design of new photofunctional materials with improved quantum yield, color purity, and thermal stability for display and lighting technologies continues.³ The conversion of electrical energy into radiative energy is further employed in optical biosensors. Sensitive electrochemiluminescent labels are attached to biological molecules with a high selectivity to detect a myriad of analytes, e.g., in clinical and biomedical diagnostics, food process control, or environmental monitoring.⁴ 

Whereas all of the applications mentioned above focus on the radiative properties of a photofunctional material, non-radiative decay is aimed for in many technologies. Chromophores with phosphorescent light emission are also being developed for photodynamic therapy (PDT).⁵ In PDT, a photosensitizer generates a cell-damaging reactive species after activation by visible light.⁶ This targeted and non-invasive approach, with localized cancer therapy or lesion ablation as the main applications, requires a compound with low dark toxicity, a high extinction coefficient, and intersystem crossing to a long-lived triplet state after illumination.^7–9 Photo-organocatalysis and photoredox catalysis are based on the acceleration of a chemical reaction upon irradiation, in contrast to photochemical reactions, which are initiated by light.^10–12 The excited catalyst enables a non-excited compound to form a highly reactive, high energetic intermediate that would not be accessible thermally, through the formation of an adduct or ion pair, hydrogen bonding, or proton, electron or atom transfer in an ultrafast reaction at ambient conditions. Sustainability issues have played an important role in enhancing photocatalytic research with sustainable hydrogen production, the removal of chemical pollutants, and photochemical synthesis as main applications.^13,14 Hydrogen has the potential to replace fossil fuels as a renewable and clean energy source with a high energy capacity, but also is an important feedstock for the chemical industry.^15–17 The high activation energy for the chemical reaction can also be provided by a photosensitizer, which absorbs and transfers photoenergy to the reaction center instead of actively participating in the chemical transformation. A successful photosynthesizer is stable to photodecomposition, has long-lived excited states, and efficiently absorbs a broad range of visible light, and the redox potential is suited for electron transfer to the catalytic center. The overall yield of a photosensitizer results from the interplay of competing mechanisms, which can impair the catalytic hydrogen production.

Despite the numerous technologies that have been successfully commercialized in the past decades, researchers are challenged to develop more efficient, robust, and sustainable chromophores that can be tailored for a specific application. A deep insight into the underlying physics and chemistry is required to design new compounds and assemblies. In particular, understanding the mechanisms and nature of the excited states is crucial to enhance the field.^10,14 While the spatial and temporal resolution of electrochemical measurements and time-resolved spectroscopy keeps improving, computational simulations are still paramount to fully comprehend the evolution of the geometric and electronic structure on ultrashort time and length scales.

Non-adiabatic molecular dynamics (NAMD) is the computational technique to model the non-radiative photodynamics of a wide variety of systems.^18–20 The approach simultaneously propagates the nuclear and electronic degrees of freedom, since multiple potential energy surfaces (PES) are involved in the deactivation process. Several NAMD methods have been and continue to be developed.^21–23 Ideally, the molecular wave function would be evolved, but including nuclear quantum effects is computationally too expensive for most compounds.²⁴ Instead, a mixed quantum/classical approximation is often applied, meaning that nuclei are propagated classically alongside the quantum mechanical propagation of the electronic structure. Tully’s fewest switches surface hopping (TFSSH) is a popular mixed quantum/classical algorithm, in which multiple independent trajectories are evolved.²⁵ Trajectories are assigned to a particular state and can transition stochastically between excited electronic states, depending on the coupling of the states’ PES. Statistically averaged results can be gained from a large number of trajectories initiated from a relevant distribution of points in the configuration space, while also allowing for uncommon mechanisms. The nuclear dynamics can either be pre-computed approximately based on a ground state PES, as done in the neglect of back-reaction approximation (NBRA) or classical path approximation (CPA),^18,26–29 or determined on-the-fly based on the gradient of the active electronic state. The former is typically used to examine the photodynamics of rigid compounds, e.g., solid-state systems, which are characterized by minimal geometrical deformation following the photoabsorption.¹⁹ 

The NAMD methods are combined with an electronic structure method to provide the ingredients for the propagation of the electron dynamics creating a plethora of possible methodologies, ranging from very accurate and computationally expensive to applicable for large-scale problems.³⁰ Density functional theory (DFT) is the most widely used electronic structure method to optimize a ground state electronic structure. Excited states can be obtained by perturbing the ground state density with a weak external electromagnetic field, either by evolving the perturbed density in time or from its corresponding linear response tensors. The two approaches are better known as real-time (RT-) and linear-response (LR-) time dependent DFT (TDDFT).^31,32 Alternatively, the ground state electronic density can be readapted and variationally optimized until it mimics the density of a chosen excited electronic state. This is the foundation of, e.g., the DFT-based delta self-consistent field (ΔSCF) or n-order constricted variational DFT [CV(n)-DFT] methods.³³ Nevertheless, Ziegler et al. have demonstrated that the main equations of LR-TDDFT can also be obtained variationally.³⁴ In this Perspective, we focus on the state-selective DFT-based ΔSCF electronic structure method.

Since many biologically or technologically interesting compounds operate or are processed in solution, various approaches for including a condensed phase environment in NAMD calculations have been developed.³⁵ The aqueous environment is key to the function and structure of biomolecules.³⁶ Solute–solvent interactions, including hydrogen bonding and electrostatic or hydrophobic effects, also determine the dynamics. Hydrogen bonding and geometrical constraints from the surrounding solvent were shown to considerably influence the excited state lifetime of various proteins, organic, and inorganic compounds.^37–39 The relative energies of excited electronic states, deactivation mechanisms, and excited state lifetimes can be significantly altered in the condensed phase.^35,40

A particular challenge for simulating TFSSH dynamics in solution, with an instantaneous change in electron density following a hop from one electronic state to another, is the sudden shift in the electronic density of the solute induced by the population change to which the solvent needs to adapt. Although the electronic polarization of the solvent quickly synchronizes, the structural adjustments to the new solute density take place more slowly, whereby the rearrangement of the outer solvation shells could take several ps depending on the solvent’s viscosity.^36,41 For example, the rearrangements can include a shift in the solvent orientation and hydrogen bonding network in case of a protic solvent, which in return influences the NA dynamics by effecting the stability of various electronic excited states of the solute in a distinct manner. Hence, it is essential to assess the errors introduced by the solvent model used and how these affect the insights gained from the NAMD simulations.

Each excited electronic state is independently optimized and the correct multiplicity is obtained by restricting the individual electron spin densities. For example, a singlet multiplicity is reproduced by equalizing the alpha and beta spin densities in a restricted KS formulation. Although several groups performed unrestricted calculations, the spin contamination was found to be an important source of error.^44,45 Alternatively, a spin purification procedure can be applied.^46,47 When the MO occupation numbers are rounded, giving one dominant occupied to virtual orbital transition, the wave function of a single reference excited electronic state is approximated hereby using a single Slater determinant.⁴⁸ As ΔSCF is a variational method, the optimized ground state density and MO occupation numbers of the targeted excited state need to be provided. A TDDFT calculation can be performed to indicate the MO occupation and reference character of the excited electronic states if no chemically substantiated guess is available.^48,49

Good agreements with respect to experiments and other computational techniques have been obtained using the ΔSCF method for a variety of systems and problems.^{45,46,50–53} The ΔSCF method is especially attractive to examine Rydberg excitations, double excitations, and charge-transfer states, which are known to be problematic cases for TDDFT.^{49,51,54–59} Cheng et al. linked the reasonable Rydberg energies obtained with semilocal exchange–correlation functionals to the qualitatively correct asymptotic form of the effective potential.⁵⁴ An overview of the excitation energy errors for a set of closed-shell, single-determinant doubly excited states using different exchange–correlation functionals is shown in Fig. 1(a).^49, Figure 1(b) compares the TDDFT and ΔSCF triplet excited state energies of a NH₃–F₂ complex at 1000 Å separation with various functionals. The charge transfer excitation energy is well approximated by ΔSCF with a low functional sensitivity, whereas a wider range of values is obtained with TDDFT. In addition, the ΔSCF electronic structure method achieves good accuracy for core excitations,^60,61 in particular, if the core holes can be localized.⁴⁹ Although the self-interaction error in modeling core-ionized states is notable, hybrid functionals do not outperform semilocal functionals.⁶² Several authors noted that the ΔSCF method generally has a low functional sensitivity.^{49,52,54,56,59,63,64} Functionals with more empirical fitting of parameters appear to be less suited for core excitations or double excitations and corresponding ΔSCF calculations suffer more from variational collapse.^49,62,65 However, the energy estimate of the lowest doubly excited state of the H₂ molecule was reported to improve with the percentage of Fock exchange,⁵⁸ and a correlation of the accuracy of vertical excitation energies for pure double excitations to the functional delocalization error was observed for some compounds.⁶³ Maurer and Reuter demonstrated that semilocal functionals lead to an underestimation of vertical excitation energies for azobenzene.⁶⁶ Hait and Head-Gordon suggested that double hybrid functionals could possibly improve the results, albeit the higher computational cost.⁴⁹ 

Another challenge for TDDFT is the description of conical intersections (CIs).⁶⁷ To the best of our knowledge, no thorough theoretical study on the accuracy of ΔSCF in regions of conical intersections has been published yet. Barca et al. investigated a conical intersection of the H₃ molecule and of retinal with their initial maximum overlap method, which gave satisfactory results.⁵⁸ The singlet excited state PES topologies of azobenzene were also consistent with the coupled cluster singles and doubles (RI-CC2) results, even though long-range corrected hybrid functionals were required for accurate excitation energies.⁶⁶ As presented in Fig. 2, the CI in the rotation pathway and lack of degeneracies in the inversion pathway of azobenzene’s molecular switching were correctly reproduced. Moreover, ΔSCF proved to be robust and numerically efficient in handling state crossings, whereas an unidentified state-crossing in TDDFT previously led to wrong conclusions.⁶⁸ Pradhan et al. recorded an approximately zero curvature of ethylene’s excited state surfaces near the conical intersection calculated using ΔSCF with fractional occupation numbers, in contrast to discontinuous derivatives computed with TDDFT.⁶⁹ Albeit improved topologies as compared to TDDFT, the PES were downshifted and a localized negative energy gap near the conical intersection region was still detected with a pure exchange–correlation density functional. The use of hybrid functional breaks the degeneracy by lifting the excited state but preserves the overall PES topology.⁴⁸ The direct optimization maximum overlap method (MOM) was proposed to improve the convergence at conical intersections.⁷⁰ 

Whereas the first studies focused on single point calculations of small molecules, the method was soon applied to small nanoclusters, spectroscopic calculations, ab initio MD, and supramolecular complexes.^42,71–75 Experimental core excitation energies could be reproduced, while TDDFT underestimated the experimental values, as well as pre-edge features in x-ray absorption spectra.⁷⁶ However, the basis set and exchange–correlation density functional used were found to significantly impact pump–probe near edge x-ray absorption fine structure (PP-NEXAFS) spectra computed with the transition-potential ΔSCF method.⁷⁷ The lowest excited state energies of multiple-resonance type thermally activated delayed fluorescence (MR-type TADF) molecules, which can be applied as organic light emitting diode (OLED) emitters, approximated experimental values, which encouraged the authors to successfully simulate the reverse intersystem crossing process of an emitter using spin–orbit coupling (SOC) constants calculated with TDDFT.⁶⁴ The charge-transfer states of a benchmark set of TADF emitters in solutions were further surveyed using ΔSCF in combination with a polarizable-continuum model.⁵⁹ Whereas the performance of different range-separated functionals was surprisingly similar, giving robust adiabatic energy gaps between the relaxed lowest singlet and triplet excited states with respect to the amount of exact exchange and the range-separation parameter, restricted open Kohn-Sham (ROKS) outperformed the unrestricted Kohn-Sham (UKS) approach. ROKS ΔSCF with semilocal hybrid functionals even proved to be accurate for this type of compounds.⁷⁸ Gavnholt et al. extended the ΔSCF method to linear combinations of KS orbitals thereby improving the accuracy of excitation energies for molecules near surfaces with highly hybridized MOs, as the charge can be localized on the adsorbed molecules.⁷⁹ Their method was further generalized for large adsorbates.⁸⁰ Behler et al. developed a similar method, called locally constrained DFT, which also converges to the ΔSCF result in the limit of infinite separation between the molecule and the surface.⁸¹ ΔSCF with an additional orthogonality condition gave reliable charge transfer excitation energies for several hydrocarbon complexes of varying sizes, including a carotenoid–porphyrin–buckyball triad, which is typically underestimated by TDDFT using standard functionals due to the negligible overlap between donor and acceptor orbitals.^82,83 Liu et al. developed the vertical self-consistent reaction field solvation model based on the ΔSCF methodology to predict solvatochromism of solvent-sensitive dyes.⁸⁴ 

ΔSCF was first applied to model a photogenerated electron–hole pair in rutile using the Ehrenfest dynamics approach.^85,86 Later, trajectory surface hopping studies of small molecular systems were performed giving good agreement with the experiment.²⁸ In a successive paper, ΔSCF was shown to be capable of modeling photoinduced nuclear reorganization dynamics with improved PES topologies as compared to TDDFT, while the accuracy of the energy values was found to be exchange–correlation functional dependent.⁶⁹ Excited state lifetimes and photoisomerization quantum yields were consistent with previous studies. Fractional occupation numbers were used in that study, and the orbital populations were updated periodically during the SCF procedure.⁶⁹ Environmental effects were added by means of a Nosé–Hoover chain thermostat in combination with a uniform velocity rescaling scheme. This only accounts for energy dissipation or redistribution, but not for polarization effects between the system and the environment. Experimental photolysis data of oxirane and the corresponding Gomer–Noyes mechanism were also confirmed with ΔSCF TFSSH simulations.⁸⁷ Discrepancies were assigned to the inaccurate energy dissipation description, which significantly impacted the dissociation yield.

The photo-induced non-radiative energy dissipation in a MoSe₂ bilayer has been studied using ultrafast electron diffraction and NAMD.⁸⁸ Although LR-TDDFT was used for the propagation of the trajectories, phonon dispersion curves for different electronic excitations were calculated using the ΔSCF electronic structure method, whereby one and two electrons were promoted from the valence to the conduction band. An efficient energy transfer of the photoenergy to the lattice vibrations in the case of a high charge carrier density was observed experimentally and reproduced computationally. The NAMD also revealed a subsequent increase in lattice temperature via softened phonon modes in the excited state, as the electronic energy is converted into atomic motion with a high efficiency. The light-induced structural phase transitions in a MoSe₂ monolayer following the promotion of eight electrons to the conduction band were further examined using an analog methodology.⁸⁹ The ΔSCF method was also used to test the neglect of back-reaction (NBRA) approximation as applied to model the photoinduced interfacial charge transfer from zinc phthalocyanine to ultrathin graphene.⁹⁰ More applications of the ΔSCF approach for the NAMD of solid-state systems are expected with the development of software packages.^91–93 

Despite the good performance of ΔSCF, there is no underlying theory to justify the method as used in applications. Gunnarsson and Lundqvist generalized the proof of Hohenberg and Kohn to the energetically lowest state of each symmetry.⁹⁴ A formal justification in the KS formalism was given by Görling, but special orbital-dependent exchange–correlation potentials for excited states should be used, which are currently still unknown.⁹⁵ Cheng et al. argued that the excited states are rigorous approximations to the true excited states if they deliver the minimum noninteracting kinetic energy for the corresponding electronic density.⁵⁴ In addition, Kowalczyk et al. demonstrated that the exact stationary density within the adiabatic approximation is formally obtained in ΔSCF,⁴⁶ while Park et al. underlined that ΔSCF is a higher order extension of adiabatic time dependent theory.⁹⁶ 

There are currently several shortcomings that limit the applicability of ΔSCF. First, the use of electronic configurations represented by a single Slater determinant limits the method to the description of excited states with single reference character. Mixed states can be approximated by a linear combination of single Slater determinants, called the sum-method.⁴³ Fractional occupation numbers used for symmetrical molecules with (nearly) degenerate orbitals also aid in the convergence.⁶⁹ In particular, geometries at the CI region or large conjugated compounds cannot be accurately described.

Second, the convergence of the ΔSCF method can be problematic. The method is especially not suited for systems with many close-lying excited states.⁴³ The regular DFT convergence schemes can also be extended to the excited state.⁹⁷ In addition, the MOM method can be used to prevent variational collapse to the ground state^{44,58,74,76,98} Levi et al. further extended the MOM method with a direct optimization approach,^70,99 while Macetti and Genoni proposed the initial-MOM approach and coupled it to the QM/ELMO (quantum mechanics/extremely localized molecular orbitals) embedding method to study larger biological systems, also in the condensed phase.¹⁰⁰ Ehlert and Klamroth only applied the overlap criterion to a subset of user-defined orbitals.⁷⁷ Several alternative methods have been developed to avoid the variational collapse, including constricted variational theory (CV-DFT),¹⁰¹ its extensions SCF-CV(∞)-DFT³³ and relaxed SCF-CV(n)-DFT,¹⁰² orthogonality constrained DFT (OCDFT),¹⁰³ single orbital replacement relaxed CV(∞)-DFT,^96, σ-SCF,¹⁰⁴ and excited constrained DFT (XCDFT).⁶⁵ The aforementioned methods are related to ΔSCF and do not significantly improve the performance if the excitation corresponds to a single orbital displacement.^33,102

Third, whereas TDDFT simultaneously provides all lowest excited states of a specified multiplicity solely based on the ground state density, targeted states are individually optimized with ΔSCF. As mentioned above, the nature of the states has to be known upfront and their convergence depends on the quality of the initial guess.¹⁰⁵ Specification of fragment charges also proved to be necessary to predict the charge transfer state of a supersystem with large interfragment separation.⁴⁹ Optimization of excited state densities dissimilar from the ground state density is challenging. The separate optimization of excited states might be cumbersome if many states are required. A final consideration to be made is the nonorthogonality of the states, which can lead to a partial variational collapse and affect, e.g., the oscillator strengths.^49,106 Nonetheless, a well-selected initial guess and restoring the translational invariance of the transition dipole moment can alleviate these issues, respectively.⁴⁵ 

The ΔSCF electronic structure method is particularly interesting for condensed phase NAMD simulations, since excited electronic states can be obtained at the cost of a ground-state calculation with an accuracy often analog to TDDFT.^83,107 Moreover, specific excited states of the photoactive system can be selected, and their energy evolution can be monitored throughout the dynamics, in contrast to the often-used TDDFT method, where only the lowest energy excited states are optimized at each point of the trajectory. However, as a consequence, while TDDFT clearly reveals when states cross or when the character of a state changes, this can be overlooked with the ΔSCF method. Even though the single reference character of a state might change during the dynamics, good results were obtained for several systems (as described in Sec. III), also at CI regions, where states have a multi-reference character and TDDFT typically does not converge.

In contrast to solvent continuum models, which allow us to approximately include solvent effects in an implicit manner,^35,114 the quantum mechanics/molecular mechanics (QM/MM) approach is a commonly applied methodology to explicitly include solvent molecules for the modeling of non-adiabatic processes in solution, in particular, for biological systems.^{20,35,115–117} Whereas the observables to propagate the dynamics of the solute are obtained from a quantum-mechanically based method, the interactions of the solvent molecules are described classically using a force-field approach.¹¹⁸ The explicit inclusion of a condensed phase environment naturally adds dynamical solvation effects to the NAMD. Several approaches to couple the QM and MM regions exist, e.g., mechanical embedding, electrostatic embedding, or polarizable embedding.^119–122 

The ΔSCF electronic structure method has been used in a QM/MM NAMD study of a diplatinum complex in aqueous solution using the electrostatic embedding scheme.¹²³ [Pt₂(P₂O₅H₂)₄]⁴⁻ is a prototype system to study the photoreactivity of d8–d8 binuclear complexes. The UV–vis absorption spectrum displays a clear HOMO to LUMO transition, corresponding to a dσ* antibonding to pσ bonding metal–metal excitation. Both the first singlet and triplet excited states are of this character, characterized by reduced Pt–Pt bond lengths. In addition, the vibronic progression of the UV–vis absorption bands established that the PES of the first singlet (S1) and triplet (T1) excited states are parallel along the Pt–Pt coordinate. The lifetime of the T1 state is in the microseconds range, while the experimental intersystem crossing time is in the picoseconds range, depending on the temperature, making this compound interesting for photocatalytic applications. Nevertheless, since the mechanism of the photoevolution cannot be determined experimentally, the ΔSCF/MM NAMD provided valuable insights into the condensed phase photodynamics in the S1 state. A Gaussian smearing of the discrete orbital occupation numbers was applied to improve the convergence. The vibrational cooling in an aqueous environment observed in experiments and parallel S1–T1 PES along the Pt–Pt coordinate were reproduced by the simulations using the grid-based projector-augmented wave (GPAW) density functional theory method. The singlet excited state geometry exhibited a symmetry lowering induced by the distortion of the PtP₄ moieties toward a quasitrigonal bipyramidal geometry. Vibrational analyses based on the decomposition of time-dependent vibrational kinetic energy in generalized normal modes indicated that the excess vibrational energy of the transition-metal atoms is distributed to the ligands, before reaching the solvent. This intramolecular vibrational energy redistribution to the ligand deformation modes was found to be favored by the surrounding solvent molecules, which increased the anharmonic couplings between the Pt–Pt pinching motion and other intramolecular vibrational modes. However, simulation times were too short to examine any change in the electronic excited state.

The TFSSH method in combination with ΔSCF and the Fireball/Amber QM/MM approach with electrostatic embedding was applied to study the photochemistry of the thymine dimerization.¹²⁴ The Fireball method speeds up DFT computations by storing and interpolating orbital integrals.¹²⁵ Since the cyclobutane thymine dimer is the most frequent DNA photolesion product, the study gave new insights into the absorption effects of ultraviolet radiation by DNA. More specifically, exposure of DNA to sunlight can lead to the formation of two C–C bonds between adjacent thymine bases, resulting in the dimer. Experiments on oligonucleotides revealed that the photochemical reaction is completed in 1 ps, while the yield depends on the thymine conformation. Previous computational studies on simplified compounds identified a barrierless pathway from the first singlet excited state to the ground state. By using the methodology mentioned above, a more realistic model for DNA at physiological conditions could be examined, including ten DNA base pairs in an aqueous environment with Na⁺ counterions. A ΔSCF QM region of 52 atoms was defined, containing the thymine nucleobases and corresponding deoxyriboses. The MM region contained ∼11 000 atoms and was evolved using the Amber and TIP3P force fields. An ensemble of 2560 trajectories was propagated for 1500 fs, starting from the region of reactive conformations. The simulations revealed that the double helix structure of DNA induces a free energy barrier between the reactive and nonreactive conformations, as shown in Fig. 4. At the conical intersection seam, one covalent C–C bond is formed. The subsequent ultrafast nonradiative decay to the ground state mostly (i.e., 96%) resulted in the standard B-DNA conformation due to the sloped topography of the CI region and the nonequilibrium nature of the process. Hence, only a very small fraction (i.e., 4%) formed the cyclobutane thymine dimer, explaining the remarkable stability of DNA against photodamage.

Implicit solvation methods and QM/MM cannot describe the changes in the electronic structure of the solvent. Consequently, no electron or proton transfers can be observed during the photodynamics. In addition, the inclusion of complex (dynamic) interactions, such as the hydrogen bonding network of protic solvents and solute–solvent polarization effects, is the major challenge. As a result, important state intersections might be overlooked, yielding incorrect lifetimes and mechanisms insignificant for the condensed phase. Full-atomistic studies on the radiationless relaxation of chromophores surrounded by small water clusters provided valuable insights into the important role of solvent molecules on the deactivation process.^126,127 Albeit more correct than implicit or classical approaches, cluster models having a site-specific addition of solvent molecules still cannot take into account the effects of bulk solvent or solvation shells further away from the solute in a sophisticated way. Moreover, the number of solvent molecules added and the solute–solvent configurations investigated critically influence the deactivation mechanisms. Periodic, explicit quantum-mechanical solvation studies for a highly sophisticated treatment of the solvation environment are, in general, hindered by the computational costs of the accurate electronic structure methods needed to properly treat NAMD. However, the photochemistry of small molecules in the condensed phase has been modeled full-atomistically using the efficient ΔSCF electronic structure method.

The NAMD propagation using ΔSCF was shown to be an efficient method for studying condensed-phase systems with periodic boundary conditions. Nieber and Doltsinis performed a primordial TFSSH study in the Born–Oppenheimer mode of uracil in fully explicit aqueous solution using a plane wave code.¹²⁸ The periodic box contained 39 water molecules. A monoexponential fit of the excited state population returned lifetimes of 608 and 359 fs in the gas phase and solution, respectively. Although the solvent generally did not affect the energetic and geometrical parameters of uracil, the out-of-plane torsional motion along the C=C bond was damped by the surrounding water. Internal conversion was further accelerated by the hydrogen bonds between the solute and the solvent.

Mališ and Luber examined the photodynamics of trans-diimide surrounded by 62 water molecules and cis-diimide in a periodic box containing 29 solvent molecules using the Landau–Zener procedure in combination with the restricted Kohn–Sham orbital formulation of the ΔSCF method relying on the efficient Gaussian and plane waves approach.^48,129 As shown in Fig. 5, the ΔSCF method accurately reproduces the topology of the S1 PES with respect to the torsional angle between the two NH groups in comparison to TDDFT and complete active space self-consistent field (CASSCF) values. Since the first solvation shell was determined to consist of 14–18 solvent molecules, the system size can capture the condensed phase environment of the solutes well. Ground-state DFT-MD calculations revealed that the diimide geometry and state energies were only slightly perturbed after solvation with weak hydrogen bonding to the surrounding solvent molecules. Nevertheless, a significant change in the photodynamics in solution with respect to vacuum was observed for both conformers (see Fig. 6). The trans-diimide showed a reduced hopping rate as compared to vacuum NAMD results. The fast deactivation mechanism was absent in the condensed phase, since the hydrogen bonding network hindered the torsional motion necessary to reach the CI. The cis-diimide trajectories, which were evolved five times longer than the bigger trans-diimide system, exhibited similar deactivation mechanisms to those in vacuum, but with a prolonged lifetime. This was explained by the thermostatic effect of the aqueous environment, in combination with a lower initial kinetic energy. Initial conditions for the solvated compounds were randomly sampled from a 15 ps NVT DFT-MD trajectory. Hence, the trajectories originated from a narrower part of the configuration space as compared to the gas-phase, where initial configurations were sampled from the thermalized Wigner distribution.^130–132 The reduced nuclear velocities result in a smaller coupling between the electronic states and consequently longer decay times. Whereas the hydrogen-bonding pattern is preserved in the first singlet excited state, a disruption following the hop to the ground state was found from evolutions of corresponding radial distribution functions (see Figs. 7 and 8). However, the first solvation shell remained intact during the cooling of the vibrationally excited state.

The same methodology in combination with the TFSSH algorithm was also applied to study cyclopropanone and its hydrate in an aqueous environment.¹¹³ The periodic simulation box contained 25 water molecules, while the 250 initial configurations were evenly sampled from a 45 ps NVT DFT-MD trajectory. The NAMD was evolved for 250 fs. In vacuum, a cyclopropanone molecule is known to dissociate forming ethylene and carbon monoxide after photoabsorption.^133–135 A static analysis of the symmetric dissociative PES gave an analog electronic energy evolution as compared to CASSCF, CASPT2, and TDDFT results.

In solution, however, a reduced photodissociation and shorter lifetime were computed. Moreover, whereas the ground state of cyclopropanone is stabilized by the surrounding solvent molecules, the first singlet excited state of the n → π* character is destabilized. In water, cyclopropanone is in equilibrium with cyclopropanone hydrate, the product of the nucleophilic addition reaction of a water molecule resulting in a reduced ring-strain. Although cyclopropanone hydrate has been studied in the role of an enzyme inhibitor,^136–139 its photochemistry has remained unresolved. NAMD using ΔSCF disclosed a range of photoproducts, as summarized in Fig. 9. After excitation to the first singlet excited state, one C–C bond is broken. Following a hop to the ground state, the bond was restored in 59% of the trajectories. Alternatively, an inter- or intramolecular proton transfer yielded a linear diol or carboxylic acid. Finally, the second C–C bond was broken in 15% of the NAMD trajectories. Subsequent proton transfer resulted in ethylene and formic acid. Interestingly, several trajectories showed a proton transfer from the diol to the surrounding solvent, initiating the Grotthuss proton transfer mechanism in water.¹⁴⁰ 

The computational cost of the trajectory surface hopping method applied to full-atomistic condensed phase NAMD simulations can be reduced by combining the variational ΔSCF method with subsystem density embedding.¹⁴¹ In subsystem DFT, the electron density of a system is approximated as a sum of subsystem densities.^142–145 The total energy of each of the subsystems can be minimized independently and simultaneously,^142,143 whereby the systems are coupled by an effective potential that is included in the optimization with the subsystem densities. The subsystems are often comprised of one or multiple molecules, in order to have only weak interactions between the units. Since subsystem DFT can significantly reduce the computational cost, while retaining a good accuracy, it has gained much interest for the modeling of complex systems. For example, subsystem DFT is well suited to study large condensed phase systems. The reduced computational cost also allows for the inclusion of more solvent molecules as compared to explicit quantum mechanical solvation studies. In addition to the computational efficiency, the analysis is simplified by the focus on the excitations of the solute subsystem.¹⁴⁵ 

Until now, only one study has been published dealing with ΔSCF in the context of subsystem DFT. This periodic subsystem density embedding in a ΔSCF formulation for NAMD was tested on a solvated diimide box using various partitioning schemes (shown in Fig. 10) in combination with the Landau–Zener and TFSSH algorithms.¹⁴¹ While the solute was always treated separately, the solvent was randomly split into groups containing 9, 3, or 1 water molecule each.

The excitation character and overall energy profile of the electronic states along the trajectory were preserved for the embedded systems as compared to full-atomistic results obtained at the same level of theory but without the application of an embedding scheme. Although an upward shift of the separate energy curves of a representative NAMD trajectory was found for all embedded systems, the energy differences between the ground and first singlet excited states were analog to the reference results. Values for the non-adiabatic coupling between the electronic states were slightly smaller for the embedded systems with a steeper diminishment. On average, the gradients from embedded and non-embedded trajectories followed the same trend, while the individual gradient vectors of the diimide molecule in the ground state showed a close agreement to the corresponding gradients of the non-embedded calculations. Finally, the propagation of NAMD trajectories using the periodic subsystem density embedding methodology revealed a slight change in the hydrogen bonding interactions and distribution of the solute conformations as compared to non-embedded trajectories, but analog lifetimes as shown in Fig. 11.

Several random selection schemes yielded equivalent results. Moreover, the size and number of the subgroups hardly influenced the observables detailed above. In addition, a partitioning based on the inner and outer solvation shells of the central diimide molecule was tested. Since no water molecules were found to change shells, the division was kept fixed. The PES and its properties were again shown to correspond well to the values obtained from the non-embedded NAMD evolution. A reduction in computational time by a factor of two to three was obtained when using the subsystem DFT with three water molecules per subsystem. Hence, the accuracy achieved for diimide in aqueous solution was more or less comparable to the results obtained with conventional Kohn–Sham DFT for the evaluated configuration space, but with a reduced simulation time. This promising methodology has not yet been tested for more complex systems, albeit an increased gain in computational efficiency with system size can be expected.

We have shown that NAMD in combination with the ΔSCF electronic structure method gave unique insights into the non-radiative decay of various systems in solution. Inclusion of the condensed phase environment is pivotal for the predictive modeling of photochemical processes for scientific purposes and commercial applications. The relative order of excited states, deactivation mechanisms, and lifetimes can be altered by the solvent through thermostatic, electrostatic, and steric interactions. Due to the computational efficiency of ΔSCF, explicit quantum mechanical solvation studies can be performed. In contrast to the often performed studies relying on solvent continuum models or various QM/MM schemes, the consistent description of the solute and solvent naturally includes polarization effects, which can drastically influence the photodynamics. Moreover, solvent molecules or atoms can participate in the deactivation mechanism, e.g., through electron or proton transfer or the hydrogen bonding network. Full-atomistic studies on the radiationless relaxation of chromophores surrounded by small water clusters already provided valuable insights into the important role of solvent molecules on the deactivation process.^126,127 However, the number of solvent molecules added and the solute–solvent configurations investigated critically influence the deactivation mechanisms that can be observed and cluster models do not take into account effects from bulk solvent or solvation shells further away. Hence, it is often difficult to judge the accuracy of cluster models. Compared to those approaches, explicit quantum mechanical solvation studies with periodic boundary conditions are the most complete and consistent approach to model the photochemistry of solvated compounds.

Subsystem embedding methods can additionally reduce the computational cost thereby enabling NAMD of large compounds and long timescales. Several embedding methods have been tested for excited state electronic structure calculations,^{100,144,146,147} but their application to ΔSCF and NAMD is still very limited with only one TFSSH ΔSCF study in solution so far.⁴⁸ The comparison of various embedding schemes would give valuable information on their relative applicability, accuracy, and speed-up for full quantum mechanical NAMD studies of large systems in the condensed phase.

Various techniques to perform or analyze a ground state DFT calculation could be extended to the excited state with ΔSCF. Nevertheless, optimization of orbitals specific for the excited state is still challenging and convergence issues are currently limiting the applications.⁴⁹ Different theoretical developments address this problem and successfully mitigate the variational collapse to the ground state.^{52,57,62,63,70,100,104,148,149} Alternatively, a many body expansion method to improve the SCF initial guess and, therefore, reduce the computational time was recently proposed.¹⁰⁵ The development of a robust scheme to converge higher excited states would be beneficial, especially since ΔSCF is a powerful method to track a specific state within the manifold of excited electronic states in a condensed phase system.¹⁵⁰ An interest in orbital-dependent functionals for exchange–correlation energies has been repeatedly expressed, but this has not yet gained traction in the ΔSCF research.^95,151,152 The development of functionals targeting excited states would contribute greatly. Moreover, additional studies focusing on compounds with degenerate states, extending the method beyond single Slater determinants and the NAMD to non-singlet multiplicities,¹⁰⁶ are necessary for various applications. In addition, a profound investigation of the performance in the region of conical intersections would endorse its application in NAMD.

The implementation of polarizable force fields for ΔSCF could significantly impact the dynamics and further elucidate the effect a condensed phase environment has on the photoevolution of large systems. Polarizable embedding models were developed to not only consider the polarization of the inner region by the solvent, but also consider the influence of the QM charge distribution on the MM region.¹⁵³ The polarization can be included separately for each region or with a self-consistent algorithm. However, the implementation for excited state processes is not straightforward. In linear response models, accounting for the polarization of the chromophore’s transition density based on the surrounding solvent is complicated.¹²² As ΔSCF is a variational method, the implementation of the reverse polarization is more clear as compared to standard linear response electronic structure methods.

External electromagnetic fields are instrumental in many applications and experiments, e.g., photoelectronics or attosecond spectroscopy.^154,155 The photodynamics and resulting photoproducts of numerous quantum systems can be manipulated or even controlled with laser pulses.¹⁵⁶ Hence, the inclusion of external electromagnetic fields in computational NAMD studies is indispensable to reproduce and predict experimental results.^157–159 Since ground state DFT algorithms can be easily extended to the excited state using the ΔSCF method, molecule–field interactions can be directly simulated. In addition, computing the photodynamics of a solvated molecule after excitation by a shaped laser pulse is a major challenge.¹⁶⁰ NAMD trajectories using the ΔSCF method could help in shining new light on this topic and elucidating the role of the solvent using explicit solvation. Finally, the efficient ΔSCF method is ideally suited to train machine learning (ML) models. ML can improve current NA techniques and their analyses. Longer length and time scales could efficiently be simulated, while also generating more universal insights with new data processing approaches or even enabling systemic screening of the photoactivity of a group of related compounds. The computational cost of generating a training set with the currently often used multireference methods is preventing ML of reaching its full potential. With ΔSCF, training data for a more comprehensive part of the configuration space or larger systems could potentially be collected. ML models trained with ΔSCF data could revolutionize the field due to their unparalleled speed enabling the fast screening of photoactive systems.

This work was supported by the University of Zurich.

The authors have no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.