Computational studies of ultrafast photoinduced processes give valuable insights into the photochemical mechanisms of a broad range of compounds. In order to accurately reproduce, interpret, and predict experimental results, which are typically obtained in a condensed phase, it is indispensable to include the condensed phase environment in the computational model. However, most studies are still performed in vacuum due to the high computational cost of state-of-the-art non-adiabatic molecular dynamics (NAMD) simulations. The quantum mechanical/molecular mechanical (QM/MM) solvation method has been a popular model to perform photodynamics in the liquid phase. Nevertheless, the currently used QM/MM embedding techniques cannot sufficiently capture all solute–solvent interactions. In this Perspective, we will discuss the efficient ΔSCF electronic structure method and its applications with respect to the NAMD of solvated compounds, with a particular focus on explicit quantum mechanical solvation. As more research is required for this method to reach its full potential, some challenges and possible directions for future research are presented as well.

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.1 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.2 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.3 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.4 

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).5 In PDT, a photosensitizer generates a cell-damaging reactive species after activation by visible light.6 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.24 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.25 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.19 

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.30 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.33 Nevertheless, Ziegler et al. have demonstrated that the main equations of LR-TDDFT can also be obtained variationally.34 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.35 The aqueous environment is key to the function and structure of biomolecules.36 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.

ΔSCF is an efficient method to calculate the excitation energies and excited state properties of systems with discrete eigenstates.42 In DFT-based ΔSCF, the density of an excited electronic state i (ρi(r)) is obtained by fixing the Kohn–Sham (KS) non-Aufbau molecular orbital (MO) occupations (njσi) during the DFT optimization of the KS MOs (φjσi(r)) with spin state σ using the ground-state algorithms,43 
ρi(r)=jσ={α,β}njσi|φjσi(r)|2=ρiα(r)+ρiβ(r).
(1)
The sum of the orbital occupation numbers equals the total number of electrons Ne, i.e., j,σnjσi=Ne. Every KS orbital j is an eigenfunction of the KS equation with corresponding eigenenergy εji,
222+w(r)+ρi(r)|rr|dr+vxc[ρiα,ρiβ](r)φjσi(r)=εjiφjσi(r),
(2)
where the left-hand side is constituted by a sum of the kinetic energy operator for non-interacting electrons, the external electrostatic interaction, electron–electron repulsion, and exchange–correlation potentials, respectively.

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.48 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.54 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 NH3–F2 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.49 Although the self-interaction error in modeling core-ionized states is notable, hybrid functionals do not outperform semilocal functionals.62 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 H2 molecule was reported to improve with the percentage of Fock exchange,58 and a correlation of the accuracy of vertical excitation energies for pure double excitations to the functional delocalization error was observed for some compounds.63 Maurer and Reuter demonstrated that semilocal functionals lead to an underestimation of vertical excitation energies for azobenzene.66 Hait and Head-Gordon suggested that double hybrid functionals could possibly improve the results, albeit the higher computational cost.49 

FIG. 1.

(a) Comparison of the ΔSCF performance to predict doubly excited states using various functionals in combination with the aug-cc-pVTZ basis set as compared to CC3 results. (b) Lowest triplet excited state charge transfer energies from a NH3 molecule to a F2 molecule 1000 Å away calculated with TDDFT and ΔSCF using various functionals and the def2-QZVPPD basis set. The dotted line indicates the estimate from the vertical ionization potential and electron affinity computed with CCSD(T) at the complete basis set limit. Reprinted with permission from D. Hait and M. Head-Gordon, J. Phys. Chem. Lett. 12, 4517–4529 (2021). Copyright 2021 American Chemical Society.

FIG. 1.

(a) Comparison of the ΔSCF performance to predict doubly excited states using various functionals in combination with the aug-cc-pVTZ basis set as compared to CC3 results. (b) Lowest triplet excited state charge transfer energies from a NH3 molecule to a F2 molecule 1000 Å away calculated with TDDFT and ΔSCF using various functionals and the def2-QZVPPD basis set. The dotted line indicates the estimate from the vertical ionization potential and electron affinity computed with CCSD(T) at the complete basis set limit. Reprinted with permission from D. Hait and M. Head-Gordon, J. Phys. Chem. Lett. 12, 4517–4529 (2021). Copyright 2021 American Chemical Society.

Close modal

Another challenge for TDDFT is the description of conical intersections (CIs).67 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 H3 molecule and of retinal with their initial maximum overlap method, which gave satisfactory results.58 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.66 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.68 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.69 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.48 The direct optimization maximum overlap method (MOM) was proposed to improve the convergence at conical intersections.70 

FIG. 2.

PES scans of the rotational (left) and inversion (right) pathways for the photoisomerization of azobenzene. The black, red, and blue colors indicate the ground state, first and second singlet excited states, respectively, calculated with the ΔSCF (solid line), TDDFT (dotted line), and RI-CC2 (dashed line) using the PBE functional for ΔSCF and TDDFT. Reprinted with permission from R. J. Maurer and K. Reuter, J. Chem. Phys. 135, 224303 (2011). Copyright 2011 AIP Publishing LLC.

FIG. 2.

PES scans of the rotational (left) and inversion (right) pathways for the photoisomerization of azobenzene. The black, red, and blue colors indicate the ground state, first and second singlet excited states, respectively, calculated with the ΔSCF (solid line), TDDFT (dotted line), and RI-CC2 (dashed line) using the PBE functional for ΔSCF and TDDFT. Reprinted with permission from R. J. Maurer and K. Reuter, J. Chem. Phys. 135, 224303 (2011). Copyright 2011 AIP Publishing LLC.

Close modal

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.76 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.77 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.64 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.59 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.78 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.79 Their method was further generalized for large adsorbates.80 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.81 Δ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.84 

Δ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.28 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.69 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.69 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.87 Discrepancies were assigned to the inaccurate energy dissipation description, which significantly impacted the dissociation yield.

The photo-induced non-radiative energy dissipation in a MoSe2 bilayer has been studied using ultrafast electron diffraction and NAMD.88 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 MoSe2 monolayer following the promotion of eight electrons to the conduction band were further examined using an analog methodology.89 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.90 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.94 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.95 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.54 In addition, Kowalczyk et al. demonstrated that the exact stationary density within the adiabatic approximation is formally obtained in ΔSCF,46 while Park et al. underlined that ΔSCF is a higher order extension of adiabatic time dependent theory.96 

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.43 Fractional occupation numbers used for symmetrical molecules with (nearly) degenerate orbitals also aid in the convergence.69 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.43 The regular DFT convergence schemes can also be extended to the excited state.97 In addition, the MOM method can be used to prevent variational collapse to the ground state44,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.100 Ehlert and Klamroth only applied the overlap criterion to a subset of user-defined orbitals.77 Several alternative methods have been developed to avoid the variational collapse, including constricted variational theory (CV-DFT),101 its extensions SCF-CV()-DFT33 and relaxed SCF-CV(n)-DFT,102 orthogonality constrained DFT (OCDFT),103 single orbital replacement relaxed CV()-DFT,96, σ-SCF,104 and excited constrained DFT (XCDFT).65 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.105 Specification of fragment charges also proved to be necessary to predict the charge transfer state of a supersystem with large interfragment separation.49 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.45 

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.

Although the states are optimized independently and, therefore, not orthogonal, which, in principle, results in a non-Hermitian spin–orbit coupling operator, a recently developed methodology to project the excited state Slater-determinants of corresponding ΔSCF states onto a linear combination of singly excited Slater determinants allows for the calculation of observables where two excited electronic states are involved.106 The ΔSCF KS determinant (|ψi⟩) mapping the electronic density of state i can be expanded as a linear combination of configuration state functions (|ϕov,o′→v′,…⟩),
|ψi=ovCovi|ϕov+oovvCoovvi|ϕov,ov+.
(3)
The configuration state functions can be made of determinants composed of ground electronic state KS orbitals (φS0) where some electrons are excited from occupied MOs o, o′, … to virtual MOs v, v′, …. The coefficients Ci are projections corresponding to
Coo,vv,i=ϕov,ov,|ψi.
(4)
Practically, the expansion is truncated at the first order containing only singly excited determinants [configuration interaction singles (CIS) expansion]. The framework resembles the auxiliary wave function technique used in linear-response TDDFT with the Tamm–Dancoff approximation.108 This configuration interaction representation is complete in the infinite basis set limit.
Orthogonalizing the excited states with respect to the ground state simplifies the calculation of the intersystem crossings in NAMD and establishes a Hermitian spin–orbit coupling operator. Starting from the nuclear–electron magnetic interaction term in the Breit–Pauli Hamiltonian,109,110 the one-electron spin–orbit coupling operator in atomic units is
ĤSOC=α22i=1Nen=1NnZn*|rin|3r̂in×p̂iŝi,
(5)
with Ne number of electrons, Nn number of nuclei, Zn* being the nth nucleus effective nuclear charge, and α being the fine structure constant. The distance vector between the nth nucleus and electron i, with momentum operator p̂i and spin operator ŝi, is rin. Insertion of the CIS expansion for the singlet and triply degenerate triplet excited states gives the following expressions for the SOC terms:
ψSi|ĤSOC|ψTj(0)ıα22ovvCovSiCovTjφvS0|ξ̂z|φvS0oovCovSiCovTjφoS0|ξ̂z|φoS0
(6)
and
ψSi|ĤSOC|ψTj(±1)ıα223/2ovvCovSiCovTjφvS0|ξ̂x±ıξ̂y|φvS0oovCovSiCovTjφoS0|ξ̂x±ıξ̂y|φoS0
(7)
with
ξ̂k=n=1NnZn*|rn|3r̂n×k,
(8)
where k is the component of the orbital angular momentum term (k = x, y, z). Operator designates the electron’s spatial partial derivative along all three space directions. The overall SOC term can be obtained from combining the contributions of the three degenerate triplet states using the root sum of squares,111,112
ψSi|ĤSOC|ψTj=Ms=+1,0,1ψSi|ĤSOC|ψTj(Ms)2.
(9)
The ΔSCF values for the one-electron spin–orbit coupling between the singlet and triplet excited states of formaldehyde using a restricted open KS formulation were shown to match the corresponding TDDFT results well, as shown in Fig. 3.106 
FIG. 3.

Electronic energy and corresponding SOC evolution during a NAMD trajectory of formaldehyde excited to the first singlet excited state (S1) obtained using the ΔSCF (solid lines) and TDDFT (dashed lines) methods. The blue, red, orange, and purple curves correspond to the ground state (S0), first triplet (T1), S1, and second triplet (T2) excited states, respectively. The pink and brown curves denote the coupling between S0 and T1 and S1 and T2 states, respectively.

FIG. 3.

Electronic energy and corresponding SOC evolution during a NAMD trajectory of formaldehyde excited to the first singlet excited state (S1) obtained using the ΔSCF (solid lines) and TDDFT (dashed lines) methods. The blue, red, orange, and purple curves correspond to the ground state (S0), first triplet (T1), S1, and second triplet (T2) excited states, respectively. The pink and brown curves denote the coupling between S0 and T1 and S1 and T2 states, respectively.

Close modal
Analogously, the UV-absorption spectra of the first singlet excited state of a cyclopropanone molecule in vacuum and in aqueous solution were calculated at the ΔSCF level of theory by approximating the transition dipole moment (with electric dipole operator μ̂) between the ground and S1 states using the CIS expansion of the S1 excited electronic state density,
ψS0|μ̂|ψS1=2ovCovS1φoS0|μ̂|φvS0.
(10)
The experimental gas-phase spectrum was found to be well reproduced in both the position and velocity representations for the electric dipole operator.113 

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.118 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.123 [Pt2(P2O5H2)4]4− 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 PtP4 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.124 The Fireball method speeds up DFT computations by storing and interpolating orbital integrals.125 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.

FIG. 4.

PES showing the free energy barrier and photoproducts obtained from the ultrafast nonradiative decay. Reprinted with permission from Mendieta-Moreno et al., J. Phys. Chem. Lett. 7, 4391–4397 (2016). Copyright 2016 American Chemical Society.

FIG. 4.

PES showing the free energy barrier and photoproducts obtained from the ultrafast nonradiative decay. Reprinted with permission from Mendieta-Moreno et al., J. Phys. Chem. Lett. 7, 4391–4397 (2016). Copyright 2016 American Chemical Society.

Close modal

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.128 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.

FIG. 5.

Evolution of the ground state (blue) and first singlet excited state (orange) energies of diimide with respect to the torsional angle θ between the two NH groups using the ΔSCF (solid lines), TDDFT (dashed lines), and SA-CASSCF(6,4) methods. The cis isomer conformation at its ground state minimum (left) is indicated with the angle of 0°. Energies are relative to the corresponding electronic ground state value of the trans conformer (right). Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 16, 4071–4086 (2020). Copyright 2020 American Chemical Society.

FIG. 5.

Evolution of the ground state (blue) and first singlet excited state (orange) energies of diimide with respect to the torsional angle θ between the two NH groups using the ΔSCF (solid lines), TDDFT (dashed lines), and SA-CASSCF(6,4) methods. The cis isomer conformation at its ground state minimum (left) is indicated with the angle of 0°. Energies are relative to the corresponding electronic ground state value of the trans conformer (right). Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 16, 4071–4086 (2020). Copyright 2020 American Chemical Society.

Close modal
FIG. 6.

Evolution of NAMD trajectories for trans-diimide (left) and cis-diimide (right), respectively, in aqueous solution. (a) and (d) Energy difference between the ground and S1 states. (b) and (e) Average population of the S1 state. (c) and (f) The torsional angle θ between the two NH groups along each trajectory. Lines are colored red (black) if the corresponding trajectory is in the first excited (ground) electronic state. Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 16, 4071–4086 (2020). Copyright 2020 American Chemical Society.

FIG. 6.

Evolution of NAMD trajectories for trans-diimide (left) and cis-diimide (right), respectively, in aqueous solution. (a) and (d) Energy difference between the ground and S1 states. (b) and (e) Average population of the S1 state. (c) and (f) The torsional angle θ between the two NH groups along each trajectory. Lines are colored red (black) if the corresponding trajectory is in the first excited (ground) electronic state. Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 16, 4071–4086 (2020). Copyright 2020 American Chemical Society.

Close modal
FIG. 7.

Radial distribution functions for the solvated cis-diimide computed at the first (black line) and last (dashed red line) step of NAMD simulations. Images (a) and (b) show the Hdii.Owat., while (c) and (d) display the Ndii.Hwat. radial distribution functions. The distance between atom pairs (dij) designates the distance for Hdii.Owat. in images (a) and (b), and for Ndii.Hwat. in images (c) and (d). Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 16, 4071–4086 (2020). Copyright 2020 American Chemical Society.

FIG. 7.

Radial distribution functions for the solvated cis-diimide computed at the first (black line) and last (dashed red line) step of NAMD simulations. Images (a) and (b) show the Hdii.Owat., while (c) and (d) display the Ndii.Hwat. radial distribution functions. The distance between atom pairs (dij) designates the distance for Hdii.Owat. in images (a) and (b), and for Ndii.Hwat. in images (c) and (d). Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 16, 4071–4086 (2020). Copyright 2020 American Chemical Society.

Close modal
FIG. 8.

Average distribution of water oxygen atoms around the NH diimide group in diimide cis and trans conformers shown in (a) and (b), respectively. (c) Oxygen atom distribution difference between the cis and trans conformers. The position and orientation of the NH group are displayed by the ball-and-stick model. Units of box axes are in Å. For better visibility, a region [0, −2, 0, −2, 0, −4 Å] was cut from all three boxes and the density is shown on the two vertical planes. Regions colored in red (gray) indicate positive (negative) values whose absolute magnitudes exceed 0.05 Å−3. Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 16, 4071–4086 (2020). Copyright 2020 American Chemical Society.

FIG. 8.

Average distribution of water oxygen atoms around the NH diimide group in diimide cis and trans conformers shown in (a) and (b), respectively. (c) Oxygen atom distribution difference between the cis and trans conformers. The position and orientation of the NH group are displayed by the ball-and-stick model. Units of box axes are in Å. For better visibility, a region [0, −2, 0, −2, 0, −4 Å] was cut from all three boxes and the density is shown on the two vertical planes. Regions colored in red (gray) indicate positive (negative) values whose absolute magnitudes exceed 0.05 Å−3. Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 16, 4071–4086 (2020). Copyright 2020 American Chemical Society.

Close modal

The same methodology in combination with the TFSSH algorithm was also applied to study cyclopropanone and its hydrate in an aqueous environment.113 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.140 

FIG. 9.

Photoproducts of the cyclopropanone hydrate NAMD calculations with their corresponding yields. The remaining trajectories (25%) decayed to the ground state, but the intermediate structures did not form a stable product by the end of the simulation time. Proton transfer will result in one of the above products.

FIG. 9.

Photoproducts of the cyclopropanone hydrate NAMD calculations with their corresponding yields. The remaining trajectories (25%) decayed to the ground state, but the intermediate structures did not form a stable product by the end of the simulation time. Proton transfer will result in one of the above products.

Close modal

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.141 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.145 

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.141 While the solute was always treated separately, the solvent was randomly split into groups containing 9, 3, or 1 water molecule each.

FIG. 10.

Simulation boxes showing different partitions of the solvated cis-diimide system tested. The water molecules were randomly assigned to 3, 9, or 27 subsystems of equal size, respectively, while the diimide molecule (dark blue) always uniquely defined a subsystem. Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 17, 1653–1661 (2021). Copyright 2021 American Chemical Society.

FIG. 10.

Simulation boxes showing different partitions of the solvated cis-diimide system tested. The water molecules were randomly assigned to 3, 9, or 27 subsystems of equal size, respectively, while the diimide molecule (dark blue) always uniquely defined a subsystem. Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 17, 1653–1661 (2021). Copyright 2021 American Chemical Society.

Close modal

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.

FIG. 11.

Population evolution of cis-diimide’s first singlet excited state in an aqueous environment. The black curves correspond to the non-embedded system, while the red curves result from the partitioning in ten subsystems. The full and dashed lines denote the Tully fewest switches surface hopping and Landau–Zener algorithms, respectively. Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 17, 1653–1661 (2021). Copyright 2021 American Chemical Society.

FIG. 11.

Population evolution of cis-diimide’s first singlet excited state in an aqueous environment. The black curves correspond to the non-embedded system, while the red curves result from the partitioning in ten subsystems. The full and dashed lines denote the Tully fewest switches surface hopping and Landau–Zener algorithms, respectively. Adapted with permission from M. Mališ and S. Luber, J. Chem. Theory Comput. 17, 1653–1661 (2021). Copyright 2021 American Chemical Society.

Close modal

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.48 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.49 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.105 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.150 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,106 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.153 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.122 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.156 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.160 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.

1.
B.
Minaev
,
G.
Baryshnikov
, and
H.
Agren
, “
Principles of phosphorescent organic light emitting devices
,”
Phys. Chem. Chem. Phys.
16
,
1719
(
2014
).
2.
L.
Xiao
,
Z.
Chen
,
B.
Qu
,
J.
Luo
,
S.
Kong
,
Q.
Gong
, and
J.
Kido
, “
Recent progresses on materials for electrophosphorescent organic light-emitting devices
,”
Adv. Mater.
23
,
926
952
(
2011
).
3.
M.
Mauro
, “
Phosphorescent multinuclear complexes for optoelectronics: Tuning of the excited-state dynamics
,”
Chem. Commun.
57
,
5857
5870
(
2021
).
4.
M. M.
Richter
, “
Electrochemiluminescence (ECL)
,”
Chem. Rev.
104
,
3003
3036
(
2004
).
5.
Y.
Dong
,
P.
Kumar
,
P.
Maity
,
I.
Kurganskii
,
S.
Li
,
A.
Elmali
,
J.
Zhao
,
D.
Escudero
,
H.
Wu
,
A.
Karatay
,
O. F.
Mohammed
, and
M.
Fedin
, “
Twisted BODIPY derivative: Intersystem crossing, electron spin polarization and application as a novel photodynamic therapy reagent
,”
Phys. Chem. Chem. Phys.
23
,
8641
8652
(
2021
).
6.
O. J.
Stacey
and
S. J. A.
Pope
, “
New avenues in the design and potential application of metal complexes for photodynamic therapy
,”
RSC Adv.
3
,
25550
25564
(
2013
).
7.
R. R.
Allison
and
C. H.
Sibata
, “
Oncologic photodynamic therapy photosensitizers: A clinical review
,”
Photodiagn. Photodyn. Ther.
7
,
61
75
(
2010
).
8.
A.
Kamkaew
,
S. H.
Lim
,
H. B.
Lee
,
L. V.
Kiew
,
L. Y.
Chung
, and
K.
Burgess
, “
BODIPY dyes in photodynamic therapy
,”
Chem. Soc. Rev.
42
,
77
88
(
2013
).
9.
C.
Imberti
,
P.
Zhang
,
H.
Huang
, and
P. J.
Sadler
, “
New designs for phototherapeutic transition metal complexes
,”
Angew. Chem., Int. Ed.
59
,
61
(
2020
).
10.
J.
Xuan
and
W.-J.
Xiao
, “
Visible-light photoredox catalysis
,”
Angew. Chem., Int. Ed.
51
,
6828
6838
(
2012
).
11.
D.
Ravelli
,
M.
Fagnoni
, and
A.
Albini
, “
Photoorganocatalysis. What for?
,”
Chem. Soc. Rev.
42
,
97
113
(
2013
).
12.
J.
Großkopf
,
T.
Kratz
,
T.
Rigotti
, and
T.
Bach
, “
Enantioselective photochemical reactions enabled by triplet energy transfer
,”
Chem. Rev.
122
,
1626
(
2021
).
13.
D.
Ravelli
,
D.
Dondi
,
M.
Fagnoni
, and
A.
Albini
, “
Photocatalysis. A multi-faceted concept for green chemistry
,”
Chem. Soc. Rev.
38
,
1999
2011
(
2009
).
14.
L. I.
Granone
,
F.
Sieland
,
N.
Zheng
,
R.
Dillert
, and
D. W.
Bahnemann
, “
Photocatalytic conversion of biomass into valuable products: A meaningful approach?
,”
Green Chem.
20
,
1169
1192
(
2018
).
15.
A. J.
Esswein
and
D. G.
Nocera
, “
Hydrogen production by molecular photocatalysis
,”
Chem. Rev.
107
,
4022
4047
(
2007
).
16.
M.
Wang
,
Y.
Na
,
M.
Gorlov
, and
L.
Sun
, “
Light-driven hydrogen production catalysed by transition metal complexes in homogeneous systems
,”
Dalton Trans.
2009
,
6458
6467
.
17.
J.
Corredor
,
M. J.
Rivero
,
C. M.
Rangel
,
F.
Gloaguen
, and
I.
Ortiz
, “
Comprehensive review and future perspectives on the photocatalytic hydrogen production
,”
J. Chem. Technol. Biotechnol.
94
,
3049
3063
(
2019
).
18.
L.
Wang
,
A.
Akimov
, and
O. V.
Prezhdo
, “
Recent progress in surface hopping: 2011–2015
,”
J. Phys. Chem. Lett.
7
,
2100
2112
(
2016
).
19.
B.
Smith
and
A. V.
Akimov
, “
Modeling nonadiabatic dynamics in condensed matter materials: Some recent advances and applications
,”
J. Phys.: Condens. Matter
32
,
073001
(
2020
).
20.
J. P.
Zobel
and
L.
González
, “
The quest to simulate excited-state dynamics of transition metal complexes
,”
JACS Au
1
,
1116
1140
(
2021
).
21.
B. F. E.
Curchod
and
T. J.
Martínez
, “
Ab initio nonadiabatic quantum molecular dynamics
,”
Chem. Rev.
118
,
3305
3336
(
2018
).
22.
F.
Agostini
and
B. F.
Curchod
, “
Different flavors of nonadiabatic molecular dynamics
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
9
,
e1417
(
2019
).
23.
T. R.
Nelson
,
A. J.
White
,
J. A.
Bjorgaard
,
A. E.
Sifain
,
Y.
Zhang
,
B.
Nebgen
,
S.
Fernandez-Alberti
,
D.
Mozyrsky
,
A. E.
Roitberg
, and
S.
Tretiak
, “
Non-adiabatic excited-state molecular dynamics: Theory and applications for modeling photophysics in extended molecular materials
,”
Chem. Rev.
120
,
2215
2287
(
2020
).
24.
L. M.
Ibele
and
B. F. E.
Curchod
, “
A molecular perspective on Tully models for nonadiabatic dynamics
,”
Phys. Chem. Chem. Phys.
22
,
15183
15196
(
2020
).
25.
J. C.
Tully
, “
Non-adiabatic molecular dynamics alkali molecules
,”
Int. J. Quantum Chem.
40
,
299
(
1991
).
26.
A. V.
Akimov
and
O. V.
Prezhdo
, “
The PYXAID program for non-adiabatic molecular dynamics in condensed matter systems
,”
J. Chem. Theory Comput.
9
,
4959
4972
(
2013
).
27.
A. V.
Akimov
,
A. J.
Neukirch
, and
O. V.
Prezhdo
, “
Theoretical insights into photoinduced charge transfer and catalysis at oxide interfaces
,”
Chem. Rev.
113
,
4496
4565
(
2013
).
28.
A. V.
Akimov
, “
Libra: An open-source ‘methodology discovery’ library for quantum and classical dynamics simulations
,”
J. Comput. Chem.
37
,
1626
1649
(
2016
).
29.
B.
Smith
and
A. V.
Akimov
, “
A comparative analysis of surface hopping acceptance and decoherence algorithms within the neglect of back-reaction approximation
,”
J. Chem. Phys.
151
,
124107
(
2019
).
30.
E.
Tapavicza
,
G. D.
Bellchambers
,
J. C.
Vincent
, and
F.
Furche
, “
Ab initio non-adiabatic molecular dynamics
,”
Phys. Chem. Chem. Phys.
15
,
18336
18348
(
2013
).
31.
N. T.
Maitra
, “
Perspective: Fundamental aspects of time-dependent density functional theory
,”
J. Chem. Phys.
144
,
220901
(
2016
).
32.
J.
Mattiat
and
S.
Luber
, “
Electronic circular dichroism with real time time dependent density functional theory: Propagator formalism and gauge dependence
,”
Chem. Phys.
527
,
110464
(
2019
).
33.
T.
Ziegler
,
M.
Krykunov
, and
J.
Cullen
, “
The implementation of a self-consistent constricted variational density functional theory for the description of excited states
,”
J. Chem. Phys.
136
,
124107
(
2012
).
34.
T.
Ziegler
,
M.
Seth
,
M.
Krykunov
,
J.
Autschbach
, and
F.
Wang
, “
On the relation between time-dependent and variational density functional theory approaches for the determination of excitation energies and transition moments
,”
J. Chem. Phys.
130
,
154102
(
2009
).
35.
F.
Santoro
,
J. A.
Green
,
L.
Martinez-Fernandez
,
J.
Cerezo
, and
R.
Improta
, “
Quantum and semiclassical dynamical studies of nonadiabatic processes in solution: Achievements and perspectives
,”
Phys. Chem. Chem. Phys.
23
,
8181
8199
(
2021
).
36.
D.
Laage
,
T.
Elsaesser
, and
J. T.
Hynes
, “
Water dynamics in the hydration shells of biomolecules
,”
Chem. Rev.
117
,
10694
10725
(
2017
).
37.
M.
Boggio-Pasqua
,
M. A.
Robb
, and
G.
Groenhof
, “
Hydrogen bonding controls excited-state decay of the photoactive yellow protein chromophore
,”
J. Am. Chem. Soc.
131
,
13580
13581
(
2009
).
38.
M.
Ruckenbauer
,
M.
Barbatti
,
T.
Müller
, and
H.
Lischka
, “
Nonadiabatic excited-state dynamics with hybrid ab initio quantum-mechanical/molecular-mechanical methods: Solvation of the pentadieniminium cation in apolar media
,”
J. Phys. Chem. A
114
,
6757
6765
(
2010
).
39.
X.
Guo
,
Y.
Zhao
, and
Z.
Cao
, “
A QM/MM MD insight into photodynamics of hypoxanthine: Distinct nonadiabatic decay behaviors between keto-N7H and keto-N9H tautomers in aqueous solution
,”
Phys. Chem. Chem. Phys.
16
,
15381
15388
(
2014
).
40.
M.
Wohlgemuth
,
V.
Bonacić-Koutecký
, and
R.
Mitrić
, “
Time-dependent density functional theory excited state nonadiabatic dynamics combined with quantum mechanical/molecular mechanical approach: Photodynamics of indole in water
,”
J. Chem. Phys.
135
,
054105
(
2011
).
41.
G.
Prampolini
,
F.
Ingrosso
,
J.
Cerezo
,
A.
Iagatti
,
P.
Foggi
, and
M.
Pastore
, “
Short- and long-range solvation effects on the transient UV–vis absorption spectra of a Ru(II)–polypyridine complex disentangled by nonequilibrium molecular dynamics
,”
J. Phys. Chem. Lett.
10
,
2885
2891
(
2019
).
42.
A.
Hellman
,
B.
Razaznejad
, and
B. I.
Lundqvist
, “
Potential-energy surfaces for excited states in extended systems
,”
J. Chem. Phys.
120
,
4593
4602
(
2004
).
43.
R. O.
Jones
and
O.
Gunnarsson
, “
The density functional formalism, its applications and prospects
,”
Rev. Mod. Phys.
61
,
689
746
(
1989
).
44.
I.
Ljubić
,
A.
Kivimäki
, and
M.
Coreno
, “
An experimental NEXAFS and computational TDDFT and ΔDFT study of the gas-phase core excitation spectra of nitroxide free radical TEMPO and its analogues
,”
Phys. Chem. Chem. Phys.
18
,
10207
10217
(
2016
).
45.
S. B.
Worster
,
O.
Feighan
, and
F. R.
Manby
, “
Reliable transition properties from excited-state mean-field calculations
,”
J. Chem. Phys.
154
,
124106
(
2021
).
46.
T.
Kowalczyk
,
S. R.
Yost
, and
T.
Van Voorhis
, “
Assessment of the ΔSCF density functional theory approach for electronic excitations in organic dyes
,”
J. Chem. Phys.
134
,
054128
(
2011
).
47.
M. W.
Hanson-Heine
,
M. W.
George
, and
N. A.
Besley
, “
Calculating excited state properties using Kohn-Sham density functional theory
,”
J. Chem. Phys.
138
,
064101
(
2013
).
48.
M.
Mališ
and
S.
Luber
, “
Trajectory surface hopping nonadiabatic molecular dynamics with Kohn–Sham ΔSCF for condensed-phase systems
,”
J. Chem. Theory Comput.
16
,
4071
4086
(
2020
).
49.
D.
Hait
and
M.
Head-Gordon
, “
Orbital optimized density functional theory for electronic excited states
,”
J. Phys. Chem. Lett.
12
,
4517
4529
(
2021
); arXiv:2103.04573.
50.
T.
Ziegler
,
A.
Rauk
, and
E. J.
Baerends
, “
On the calculation of multiplet energies by the Hartree-Fock-Slater method
,”
Theor. Chim. Acta
43
,
261
271
(
1977
).
51.
D. J.
Tozer
and
N. C.
Handy
, “
On the determination of excitation energies using density functional theory
,”
Phys. Chem. Chem. Phys.
2
,
2117
2121
(
2000
).
52.
B.
Peng
,
B. E.
Van Kuiken
,
F.
Ding
, and
X.
Li
, “
A guided self-consistent-field method for excited-state wave function optimization: Applications to ligand-field transitions in transition-metal complexes
,”
J. Chem. Theory Comput.
9
,
3933
3938
(
2013
).
53.
W. C.
Mackrodt
,
S.
Salustro
,
B.
Civalleri
, and
R.
Dovesi
, “
Low energy excitations in NiO based on a direct Δ-SCF approach
,”
J. Phys.: Condens. Matter
30
,
495901
(
2018
).
54.
C. L.
Cheng
,
Q.
Wu
, and
T.
Van Voorhis
, “
Rydberg energies using excited state density functional theory
,”
J. Chem. Phys.
129
,
124112
(
2008
).
55.
K.
Yang
,
R.
Peverati
,
D. G.
Truhlar
, and
R.
Valero
, “
Density functional study of multiplicity-changing valence and Rydberg excitations of p-block elements: Delta self-consistent field, collinear spin-flip time-dependent density functional theory (DFT), and conventional time-dependent DFT
,”
J. Chem. Phys.
135
,
044118
(
2011
).
56.
I.
Seidu
,
M.
Krykunov
, and
T.
Ziegler
, “
Applications of time-dependent and time-independent density functional theory to Rydberg transitions
,”
J. Phys. Chem. A
119
,
5107
5116
(
2015
).
57.
J.
Liu
,
Y.
Zhang
,
P.
Bao
, and
Y.
Yi
, “
Evaluating electronic couplings for excited state charge transfer based on maximum occupation method ΔSCF quasi-adiabatic states
,”
J. Chem. Theory Comput.
13
,
843
851
(
2017
).
58.
G. M. J.
Barca
,
A. T. B.
Gilbert
, and
P. M. W.
Gill
, “
Simple models for difficult electronic excitations
,”
J. Chem. Theory Comput.
14
,
1501
1509
(
2018
).
59.
L.
Kunze
,
A.
Hansen
,
S.
Grimme
, and
J.-M.
Mewes
, “
PCM-ROKS for the description of charge-transfer states in solution: Singlet-triplet gaps with chemical accuracy from open-shell Kohn–Sham reaction-field calculations
,”
J. Phys. Chem. Lett.
12
,
8470
8480
(
2021
).
60.
N. A.
Besley
and
F. A.
Asmuruf
, “
Time-dependent density functional theory calculations of the spectroscopy of core electrons
,”
Phys. Chem. Chem. Phys.
12
,
12024
12039
(
2010
).
61.
D.
Hait
and
M.
Head-Gordon
, “
Highly accurate prediction of core spectra of molecules at density functional theory cost: Attaining sub-electronvolt error from a restricted open-shell Kohn–Sham approach
,”
J. Phys. Chem. Lett.
11
,
775
786
(
2020
); arXiv:1912.05249.
62.
K.
Carter-Fenk
and
J. M.
Herbert
, “
State-targeted energy projection: A simple and robust approach to orbital relaxation of non-Aufbau self-consistent field solutions
,”
J. Chem. Theory Comput.
16
,
5067
5082
(
2020
).
63.
D.
Hait
and
M.
Head-Gordon
, “
Excited state orbital optimization via minimizing the square of the gradient: General approach and application to singly and doubly excited states via density functional theory
,”
J. Chem. Theory Comput.
16
,
1699
1710
(
2020
).
64.
W.
Sotoyama
, “
Simulation of low-lying singlet and triplet excited states of multiple-resonance-type thermally activated delayed fluorescence emitters by delta self-consistent field (ΔSCF) method
,”
J. Phys. Chem. A
125
,
10373
10378
(
2021
).
65.
P.
Ramos
and
M.
Pavanello
, “
Low-lying excited states by constrained DFT
,”
J. Chem. Phys.
148
,
144103
(
2018
).
66.
R. J.
Maurer
and
K.
Reuter
, “
Assessing computationally efficient isomerization dynamics: ΔSCF density-functional theory study of azobenzene molecular switching
,”
J. Chem. Phys.
135
,
224303
(
2011
).
67.
B. G.
Levine
,
C.
Ko
,
J.
Quenneville
, and
T. J.
MartÍnez
, “
Conical intersections and double excitations in time-dependent density functional theory
,”
Mol. Phys.
104
,
1039
1051
(
2006
).
68.
M. L.
Tiago
,
S.
Ismail-Beigi
, and
S. G.
Louie
, “
Photoisomerization of azobenzene from first-principles constrained density-functional calculations
,”
J. Chem. Phys.
122
,
094311
(
2005
).
69.
E.
Pradhan
,
K.
Sato
, and
A. V.
Akimov
, “
Non-adiabatic molecular dynamics with ΔSCF excited states
,”
J. Phys.: Condens. Matter
30
,
484002
(
2018
).
70.
G.
Levi
,
A. V.
Ivanov
, and
H.
Jónsson
, “
Variational density functional calculations of excited states via direct optimization
,”
J. Chem. Theory Comput.
16
,
6968
6982
(
2020
); arXiv:2006.15922.
71.
D.
Ceresoli
,
E.
Tosatti
,
S.
Scandolo
,
G.
Santoro
, and
S.
Serra
, “
Trapping of excitons at chemical defects in polyethylene
,”
J. Chem. Phys.
121
,
6478
6484
(
2004
).
72.
E.
Degoli
,
G.
Cantele
,
E.
Luppi
,
R.
Magri
,
D.
Ninno
,
O.
Bisi
, and
S.
Ossicini
, “
Ab initio structural and electronic properties of hydrogenated silicon nanoclusters in the ground and excited state
,”
Phys. Rev. B
69
,
155411
(
2004
).
73.
R. R.
Zope
,
T.
Baruah
,
S. L.
Richardson
,
M. R.
Pederson
, and
B. I.
Dunlap
, “
Optical excitation energies, Stokes shift, and spin-splitting of C24H72Si14
,”
J. Chem. Phys.
133
,
034301
(
2010
).
74.
E. A.
Briggs
,
N. A.
Besley
, and
D.
Robinson
, “
QM/MM excited state molecular dynamics and fluorescence spectroscopy of BODIPY
,”
J. Phys. Chem. A
117
,
2644
2650
(
2013
).
75.
G.
Levi
,
E.
Biasin
,
A. O.
Dohn
, and
H.
Jónsson
, “
On the interplay of solvent and conformational effects in simulated excited-state dynamics of a copper phenanthroline photosensitizer
,”
Phys. Chem. Chem. Phys.
22
,
748
757
(
2020
).
76.
N. A.
Besley
,
A. T.
Gilbert
, and
P. M.
Gill
, “
Self-consistent-field calculations of core excited states
,”
J. Chem. Phys.
130
,
124308
(
2009
).
77.
C.
Ehlert
and
T.
Klamroth
, “
PSIXAS: A Psi4 plugin for efficient simulations of X-ray absorption spectra based on the transition-potential and Δ-Kohn–Sham method
,”
J. Comput. Chem.
41
,
1781
1789
(
2020
).
78.
D.
Hait
,
T.
Zhu
,
D. P.
McMahon
, and
T.
Van Voorhis
, “
Prediction of excited-state energies and singlet–triplet gaps of charge-transfer states using a restricted open-shell Kohn–Sham approach
,”
J. Chem. Theory Comput.
12
,
3353
3359
(
2016
).
79.
J.
Gavnholt
,
T.
Olsen
,
M.
Engelund
, and
J.
Schiøtz
, “
Δ self-consistent field method to obtain potential energy surfaces of excited molecules on surfaces
,”
Phys. Rev. B
78
,
075441
(
2008
).
80.
R. J.
Maurer
and
K.
Reuter
, “
Excited-state potential-energy surfaces of metal-adsorbed organic molecules from linear expansion Δ-self-consistent field density-functional theory (ΔSCF-DFT)
,”
J. Chem. Phys.
139
,
014708
(
2013
).
81.
J.
Behler
,
B.
Delley
,
K.
Reuter
, and
M.
Scheffler
, “
Nonadiabatic potential-energy surfaces by constrained density-functional theory
,”
Phys. Rev. B
75
,
115409
(
2007
).
82.
T.
Baruah
and
M. R.
Pederson
, “
DFT calculations on charge-transfer states of a carotenoid-porphyrin-C60 molecular triad
,”
J. Chem. Theory Comput.
5
,
834
843
(
2009
).
83.
T.
Baruah
,
M.
Olguin
, and
R. R.
Zope
, “
Charge transfer excited state energies by perturbative delta self consistent field method
,”
J. Chem. Phys.
137
,
084316
(
2012
).
84.
T.
Liu
,
W.-G.
Han
,
F.
Himo
,
G. M.
Ullmann
,
D.
Bashford
,
A.
Toutchkine
,
K. M.
Hahn
, and
L.
Noodleman
, “
Density functional vertical self-consistent reaction field theory for solvatochromism studies of solvent-sensitive dyes
,”
J. Phys. Chem. A
108
,
3545
3555
(
2004
).
85.
G.
Kolesov
,
D.
Vinichenko
,
G. A.
Tritsaris
,
C. M.
Friend
, and
E.
Kaxiras
, “
Anatomy of the photochemical reaction: Excited-state dynamics reveals the C–H acidity mechanism of methoxy photo-oxidation on titania
,”
J. Phys. Chem. Lett.
6
,
1624
1627
(
2015
).
86.
G. A.
Tritsaris
,
D.
Vinichenko
,
G.
Kolesov
,
C. M.
Friend
, and
E.
Kaxiras
, “
Dynamics of the photogenerated hole at the rutile TiO2(110)/water interface: A nonadiabatic simulation study
,”
J. Phys. Chem. C
118
,
27393
27401
(
2014
).
87.
M.
Krenz
,
U.
Gerstmann
, and
W. G.
Schmidt
, “
Photochemical ring opening of oxirane modeled by constrained density functional theory
,”
ACS Omega
5
,
24057
24063
(
2020
).
88.
M. F.
Lin
,
V.
Kochat
,
A.
Krishnamoorthy
,
L.
Bassman
,
C.
Weninger
,
Q.
Zheng
,
X.
Zhang
,
A.
Apte
,
C. S.
Tiwary
,
X.
Shen
,
R.
Li
,
R.
Kalia
,
P.
Ajayan
,
A.
Nakano
,
P.
Vashishta
,
F.
Shimojo
,
X.
Wang
,
D. M.
Fritz
, and
U.
Bergmann
, “
Ultrafast non-radiative dynamics of atomically thin MoSe2
,”
Nat. Commun.
8
,
1745
(
2017
).
89.
L.
Bassman
,
A.
Krishnamoorthy
,
H.
Kumazoe
,
M.
Misawa
,
F.
Shimojo
,
R. K.
Kalia
,
A.
Nakano
, and
P.
Vashishta
, “
Electronic origin of optically-induced sub-picosecond lattice dynamics in MoSe2 monolayer
,”
Nano Lett.
18
,
4653
4658
(
2018
).
90.
H.
Mehdipour
,
B. A.
Smith
,
A. T.
Rezakhani
,
S. S.
Tafreshi
,
N. H.
de Leeuw
,
O. V.
Prezhdo
,
A. Z.
Moshfegh
, and
A. V.
Akimov
, “
Dependence of electron transfer dynamics on the number of graphene layers in π-stacked 2D materials: Insights from ab initio nonadiabatic molecular dynamics
,”
Phys. Chem. Chem. Phys.
21
,
23198
23208
(
2019
).
91.
K.
Sato
,
E.
Pradhan
,
R.
Asahi
, and
A. V.
Akimov
, “
Charge transfer dynamics at the boron subphthalocyanine chloride/C60 interface: Non-adiabatic dynamics study with Libra-X
,”
Phys. Chem. Chem. Phys.
20
,
25275
25294
(
2018
).
92.
F.
Shimojo
,
S.
Fukushima
,
H.
Kumazoe
,
M.
Misawa
,
S.
Ohmura
,
P.
Rajak
,
K.
Shimamura
,
L.
Bassman
,
S.
Tiwari
,
R. K.
Kalia
,
A.
Nakano
, and
P.
Vashishta
, “
QXMD: An open-source program for nonadiabatic quantum molecular dynamics
,”
SoftwareX
10
,
100307
(
2019
).
93.
Zagreb surface hopping code; available from https://e-cam.readthedocs.io/en/latest/quantum-dynamics-modules/index.html,
2018
.
94.
O.
Gunnarsson
and
B. I.
Lundqvist
, “
Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism
,”
Phys. Rev. B
13
,
4274
4298
(
1976
).
95.
A.
Görling
, “
Density-functional theory beyond the Hohenberg-Kohn theorem
,”
Phys. Rev. A
59
,
3359
3374
(
1999
).
96.
Y. C.
Park
,
M.
Krykunov
, and
T.
Ziegler
, “
On the relation between adiabatic time dependent density functional theory (TDDFT) and the ΔSCF-DFT method. Introducing a numerically stable ΔSCF-DFT scheme for local functionals based on constricted variational DFT
,”
Mol. Phys.
113
,
1636
1647
(
2015
).
97.
A.
Baldes
,
W.
Klopper
,
J.
Šimunek
,
J.
Noga
, and
F.
Weigend
, “
Acceleration of self-consistent-field convergence by combining conventional diagonalization and a diagonalization-free procedure
,”
J. Comput. Chem.
32
,
3129
3134
(
2011
).
98.
A. T. B.
Gilbert
,
N. A.
Besley
, and
P. M. W.
Gill
, “
Self-consistent field calculations of excited states using the maximum overlap method (MOM)
,”
J. Phys. Chem. A
112
,
13164
13171
(
2008
).
99.
G.
Levi
,
A. V.
Ivanov
, and
H.
Jónsson
, “
Variational calculations of excited states: Via direct optimization of the orbitals in DFT
,”
Faraday Discuss.
224
,
448
466
(
2020
).
100.
G.
Macetti
and
A.
Genoni
, “
Initial maximum overlap method for large systems by the quantum mechanics/extremely localized molecular orbital embedding technique
,”
J. Chem. Theory Comput.
17
,
4169
4182
(
2021
).
101.
J.
Cullen
,
M.
Krykunov
, and
T.
Ziegler
, “
The formulation of a self-consistent constricted variational density functional theory for the description of excited states
,”
Chem. Phys.
391
,
11
18
(
2011
).
102.
M.
Krykunov
and
T.
Ziegler
, “
Self-consistent formulation of constricted variational density functional theory with orbital relaxation. Implementation and applications
,”
J. Chem. Theory Comput.
9
,
2761
2773
(
2013
).
103.
F. A.
Evangelista
,
P.
Shushkov
, and
J. C.
Tully
, “
Orthogonality constrained density functional theory for electronic excited states
,”
J. Phys. Chem. A
117
,
7378
7392
(
2013
).
104.
H. Z.
Ye
,
M.
Welborn
,
N. D.
Ricke
, and
T.
Van Voorhis
, “
σ-SCF: A direct energy-targeting method to mean-field excited states
,”
J. Chem. Phys.
147
,
214104
(
2017
).
105.
F.
Ballesteros
and
K. U.
Lao
, “
Accelerating the convergence of self-consistent field calculations using the many-body expansion
,”
J. Chem. Theory Comput.
18
,
179
(
2021
).
106.
M.
Mališ
,
E.
Vandaele
, and
S.
Luber
, “
Spin-orbit couplings for nonadiabatic molecular dynamics at the ΔSCF level
” (unpublished).
107.
O. V.
Prezhdo
, “
Modeling non-adiabatic dynamics in nanoscale and condensed matter systems
,”
Acc. Chem. Res.
54
,
4239
(
2021
).
108.
I.
Tavernelli
,
B. F. E.
Curchod
,
A.
Laktionov
, and
U.
Rothlisberger
, “
Nonadiabatic coupling vectors for excited states within time-dependent density functional theory in the Tamm–Dancoff approximation and beyond
,”
J. Chem. Phys.
133
,
194104
(
2010
).
109.
G.
Breit
, “
Dirac’s equation and the spin-spin interactions of two electrons
,”
Phys. Rev.
39
,
616
624
(
1932
).
110.
R.
McWeeny
,
Spins in Chemistry
(
Dover Publications
,
Mineola, NY
,
2004
).
111.
X.
Gao
,
S.
Bai
,
D.
Fazzi
,
T.
Niehaus
,
M.
Barbatti
, and
W.
Thiel
, “
Evaluation of spin-orbit couplings with linear-response time-dependent density functional methods
,”
J. Chem. Theory Comput.
13
,
515
524
(
2017
).
112.
F.
Franco De Carvalho
,
C. A.
Pignedoli
, and
I.
Tavernelli
, “
TDDFT-based spin–orbit couplings of 0D, 1D, and 2D carbon nanostructures: Static and dynamical effects
,”
J. Phys. Chem. C
121
,
10140
10152
(
2017
).
113.
E.
Vandaele
,
M.
Mališ
, and
S.
Luber
, “
The photodissociation of solvated cyclopropanone and its hydrate explored via non-adiabatic molecular dynamics using ΔSCF
,”
Phys. Chem. Chem. Phys.
24
,
5669
(
2022
).
114.
J.
Tomasi
,
B.
Mennucci
, and
R.
Cammi
, “
Quantum mechanical continuum solvation models
,”
Chem. Rev.
105
,
2999
3093
(
2005
).
115.
M.
Boggio-Pasqua
,
C. F.
Burmeister
,
M. A.
Robb
, and
G.
Groenhof
, “
Photochemical reactions in biological systems: Probing the effect of the environment by means of hybrid quantum chemistry/molecular mechanics simulations
,”
Phys. Chem. Chem. Phys.
14
,
7912
7928
(
2012
).
116.
L. W.
Chung
,
W. M. C.
Sameera
,
R.
Ramozzi
,
A. J.
Page
,
M.
Hatanaka
,
G. P.
Petrova
,
T. V.
Harris
,
X.
Li
,
Z.
Ke
,
F.
Liu
,
H.-B.
Li
,
L.
Ding
, and
K.
Morokuma
, “
The ONIOM method and its applications
,”
Chem. Rev.
115
,
5678
5796
(
2015
).
117.
E.
Brunk
and
U.
Rothlisberger
, “
Mixed quantum mechanical/molecular mechanical molecular dynamics simulations of biological systems in ground and electronically excited states
,”
Chem. Rev.
115
,
6217
6263
(
2015
).
118.
A.
Warshel
and
M.
Levitt
, “
Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme
,”
J. Mol. Biol.
103
,
227
249
(
1976
).
119.
D.
Bakowies
and
W.
Thiel
, “
Hybrid models for combined quantum mechanical and molecular mechanical approaches
,”
J. Phys. Chem.
100
,
10580
10594
(
1996
).
120.
H.
Lin
and
D. G.
Truhlar
, “
QM/MM: What have we learned, where are we, and where do we go from here?
,”
Theor. Chem. Acc.
117
,
185
199
(
2007
).
121.
H. M.
Senn
and
W.
Thiel
, “
QM/MM methods for biomolecular systems
,”
Angew. Chem., Int. Ed.
48
,
1198
1229
(
2009
).
122.
M.
Bondanza
,
M.
Nottoli
,
L.
Cupellini
,
F.
Lipparini
, and
B.
Mennucci
, “
Polarizable embedding QM/MM: The future gold standard for complex (bio)systems?
,”
Phys. Chem. Chem. Phys.
22
,
14433
14448
(
2020
).
123.
G.
Levi
,
M.
Pápai
,
N. E.
Henriksen
,
A. O.
Dohn
, and
K. B.
Møller
, “
Solution structure and ultrafast vibrational relaxation of the PtPOP complex revealed by ΔSCF-QM/MM direct dynamics simulations
,”
J. Phys. Chem. C
122
,
7100
7119
(
2018
).
124.
J. I.
Mendieta-Moreno
,
D. G.
Trabada
,
J.
Mendieta
,
J. P.
Lewis
,
P.
Gómez-Puertas
, and
J.
Ortega
, “
Quantum mechanics/molecular mechanics free energy maps and nonadiabatic simulations for a photochemical reaction in DNA: Cyclobutane thymine dimer
,”
J. Phys. Chem. Lett.
7
,
4391
4397
(
2016
).
125.
J. I.
Mendieta-Moreno
,
R. C.
Walker
,
J. P.
Lewis
,
P.
Gómez-Puertas
,
J.
Mendieta
, and
J.
Ortega
, “
Fireball/amber: An efficient local-orbital DFT QM/MM method for biomolecular systems
,”
J. Chem. Theory Comput.
10
,
2185
2193
(
2014
).
126.
M.
Barbatti
, “
Photorelaxation induced by water-chromophore electron transfer
,”
J. Am. Chem. Soc.
136
,
10246
10249
(
2014
).
127.
J.
Thisuwan
,
S.
Chaiwongwattana
,
M.
Sapunar
,
K.
Sagarik
, and
N.
Došlić
, “
Photochemical deactivation pathways of microsolvated hydroxylamine
,”
J. Photochem. Photobiol., A
328
,
10
15
(
2016
).
128.
H.
Nieber
and
N. L.
Doltsinis
, “
Elucidating ultrafast nonradiative decay of photoexcited uracil in aqueous solution by ab initio molecular dynamics
,”
Chem. Phys.
347
,
405
412
(
2008
).
129.
J.
Vandevondele
,
M.
Krack
,
F.
Mohamed
,
M.
Parrinello
,
T.
Chassaing
, and
J.
Hutter
, “
Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach
,”
Comput. Phys. Commun.
167
,
103
128
(
2005
).
130.
M.
Barbatti
and
K.
Sen
, “
Effects of different initial condition samplings on photodynamics and spectrum of pyrrole
,”
Int. J. Quantum Chem.
116
,
762
771
(
2016
).
131.
J. J.
Nogueira
and
L.
González
, “
Computational photophysics in the presence of an environment
,”
Annu. Rev. Phys. Chem.
69
,
473
497
(
2018
).
132.
J. P.
Zobel
,
M.
Heindl
,
J. J.
Nogueira
, and
L.
González
, “
Vibrational sampling and solvent effects on the electronic structure of the absorption spectrum of 2-nitronaphthalene
,”
J. Chem. Theory Comput.
14
,
3205
3217
(
2018
).
133.
H. J.
Rodriguez
,
J.-C.
Chang
, and
T. F.
Thomas
, “
Thermal, photochemical, and photophysical processes in cyclopropanone vapor
,”
J. Am. Chem. Soc.
98
,
2027
2034
(
1976
).
134.
G.
Cui
,
Y.
Ai
, and
W.
Fang
, “
Conical intersection is responsible for the fluorescence disappearance below 365 nm in cyclopropanone
,”
J. Phys. Chem. A
114
,
730
734
(
2010
).
135.
G.
Cui
and
W.
Fang
, “
Ab initio based surface-hopping dynamics study on ultrafast internal conversion in cyclopropanone
,”
J. Phys. Chem. A
115
,
1547
1555
(
2011
).
136.
T. H.
Cromartie
, “
Irreversible inactivation of the flavoenzyme alcohol oxidase by cyclopropanone
,”
Biochem. Biophys. Res. Commun.
105
,
785
790
(
1982
).
137.
J. S.
Wiseman
and
R. H.
Abeles
, “
Mechanism of inhibition of aldehyde dehydrogenase by cyclopropanone hydrate and the mushroom toxin coprine
,”
Biochemistry
18
,
427
435
(
1979
).
138.
J. S.
Wiseman
,
J. S.
Nichols
, and
M. X.
Kolpak
, “
Mechanism of inhibition of horseradish peroxidase by cyclopropanone hydrate
,”
J. Biol. Chem.
257
,
6328
6332
(
1982
).
139.
K. W.
Schmid
,
A.
Hittmair
,
H.
Schmidhammer
, and
B.
Jasani
, “
Non-deleterious inhibition of endogenous peroxidase activity (EPA) by cyclopropanone hydrate: A definitive approach
,”
J. Histochem. Cytochem.
37
,
473
477
(
1989
).
140.
N.
Agmon
, “
The Grotthuss mechanism
,”
Chem. Phys. Lett.
244
,
456
462
(
1995
).
141.
M.
Mališ
and
S.
Luber
, “
ΔSCF with subsystem density embedding for efficient nonadiabatic molecular dynamics in condensed-phase systems
,”
J. Chem. Theory Comput.
17
,
1653
1661
(
2021
).
142.
M.
Iannuzzi
,
B.
Kirchner
, and
J.
Hutter
, “
Density functional embedding for molecular systems
,”
Chem. Phys. Lett.
421
,
16
20
(
2006
).
143.
S.
Luber
, “
Local electric dipole moments for periodic systems via density functional theory embedding
,”
J. Chem. Phys.
141
,
234110
(
2014
).
144.
C. R.
Jacob
and
J.
Neugebauer
, “
Subsystem density-functional theory
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
325
(
2014
).
145.
T. A.
Wesolowski
,
S.
Shedge
, and
X.
Zhou
, “
Frozen-density embedding strategy for multilevel simulations of electronic structure
,”
Chem. Rev.
115
,
5891
5928
(
2015
).
146.
W.-K.
Chen
,
W.-H.
Fang
, and
G.
Cui
, “
Integrating machine learning with the multilayer energy-based fragment method for excited states of large systems
,”
J. Phys. Chem. Lett.
10
,
7836
7841
(
2019
).
147.
S.
Andermatt
,
J.
Cha
,
F.
Schiffmann
, and
J.
VandeVondele
, “
Combining linear-scaling DFT with subsystem DFT in Born–Oppenheimer and Ehrenfest molecular dynamics simulations: From molecules to a virus in solution
,”
J. Chem. Theory Comput.
12
,
3214
3227
(
2016
).
148.
P.
Ramos
and
M.
Pavanello
, “
Nonadiabatic couplings from a variational excited state method based on constrained DFT
,”
J. Chem. Phys.
154
,
014110
(
2021
).
149.
H. G. A.
Burton
and
D. J.
Wales
, “
Energy landscapes for electronic structure
,”
J. Chem. Theory Comput.
17
,
151
169
(
2021
).
150.
P.
Zawadzki
,
K. W.
Jacobsen
, and
J.
Rossmeisl
, “
Electronic hole localization in rutile and anatase TiO2—Self-interaction correction in Δ-SCF DFT
,”
Chem. Phys. Lett.
506
,
42
45
(
2011
).
151.
A.
Görling
, “
Orbital- and state-dependent functionals in density-functional theory
,”
J. Chem. Phys.
123
,
062203
(
2005
).
152.
S.
Kümmel
and
L.
Kronik
, “
Orbital-dependent density functionals: Theory and applications
,”
Rev. Mod. Phys.
80
,
3
60
(
2008
).
153.
A.
Wildman
,
G.
Donati
,
F.
Lipparini
,
B.
Mennucci
, and
X.
Li
, “
Nonequilibrium environment dynamics in a frequency-dependent polarizable embedding model
,”
J. Chem. Theory Comput.
15
,
43
51
(
2019
).
154.
B.
Mignolet
,
B. F.
Curchod
, and
T. J.
Martínez
, “
Communication: XFAIMS–external field ab initio multiple spawning for electron-nuclear dynamics triggered by short laser pulses
,”
J. Chem. Phys.
145
,
191104
(
2016
).
155.
C.
Park
,
K.
Lee
,
M.
Koo
, and
C.
Park
, “
Soft ferroelectrics enabling high-performance intelligent photo electronics
,”
Adv. Mater.
33
,
2004999
(
2021
).
156.
G. A.
Worth
and
G. W.
Richings
, “
Optimal control by computer
,”
Annu. Rep. Prog. Chem., Sect. C: Phys. Chem.
109
,
113
139
(
2013
).
157.
R.
Mitrić
,
J.
Petersen
, and
V.
Bonačić-Koutecký
, “
Laser-field-induced surface-hopping method for the simulation and control of ultrafast photodynamics
,”
Phys. Rev. A
79
,
053416
(
2009
).
158.
J. J.
Bajo
,
J.
González-Vázquez
,
I. R.
Sola
,
J.
Santamaria
,
M.
Richter
,
P.
Marquetand
, and
L.
González
, “
Mixed quantum-classical dynamics in the adiabatic representation to simulate molecules driven by strong laser pulses
,”
J. Phys. Chem. A
116
,
2800
2807
(
2012
).
159.
Z.
Zhou
,
H.-T.
Chen
,
A.
Nitzan
, and
J. E.
Subotnik
, “
Nonadiabatic dynamics in a laser field: Using Floquet fewest switches surface hopping to calculate electronic populations for slow nuclear velocities
,”
J. Chem. Theory Comput.
16
,
821
834
(
2020
).
160.
D.
Keefer
,
S.
Thallmair
,
J. P.
Zauleck
, and
R. D.
Vivie-Riedle
, “
A multi target approach to control chemical reactions in their inhomogeneous solvent environment
,”
J. Phys. B: At., Mol. Opt. Phys.
48
,
234003
(
2015
).
Published open access through an agreement with Universitat Zurich