Combining SchNet and SHARC: The SchNarc Machine Learning Approach for Excited-State Dynamics.

In recent years, deep learning has become a part of our everyday life and is revolutionizing quantum chemistry as well. In this work, we show how deep learning can be used to advance the research field of photochemistry by learning all important properties-multiple energies, forces, and different couplings-for photodynamics simulations. We simplify such simulations substantially by (i) a phase-free training skipping costly preprocessing of raw quantum chemistry data; (ii) rotationally covariant nonadiabatic couplings, which can either be trained or (iii) alternatively be approximated from only ML potentials, their gradients, and Hessians; and (iv) incorporating spin-orbit couplings. As the deep-learning method, we employ SchNet with its automatically determined representation of molecular structures and extend it for multiple electronic states. In combination with the molecular dynamics program SHARC, our approach termed SchNarc is tested on two polyatomic molecules and paves the way toward efficient photodynamics simulations of complex systems.

Excited-state dynamics simulations are powerful tools to predict, understand, and explain photoinduced processes, especially in combination with experimental studies. Examples of photoinduced processes range from photosynthesis, DNA photodamage as the starting point of skin cancer, to processes that enable our vision.[1−5] As they are part of our everyday lives, their understanding can help to unravel fundamental processes of nature and to advance several research fields, such as photovoltaics,[6,7] photocatalysis,[8] or photosensitive drug design.[9]

Because the full quantum mechanical treatment of molecules remains challenging, exact quantum dynamics simulations are limited to systems containing only a couple of atoms, even if fitted potential energy surfaces (PESs) are used.[10−26] In order to treat larger systems in full dimensions, i.e., systems with up to hundreds of atoms, and on long time scales, i.e., in the range of several 100 ps, excited-state machine learning (ML) molecular dynamics (MD), where the ML model is trained on quantum chemistry data, has evolved as a promising tool in the last couple of years.[27−33]

Such nonadiabatic MLMD simulations are in many senses analogous to excited-state ab initio molecular dynamics simulations. The only difference is that the costly electronic structure calculations are mostly replaced by an ML model, providing quantum properties like the PESs and the corresponding forces. The nuclei are assumed to move classically on those PESs. This mixed quantum–classical dynamics approach allows for a very fast on-the-fly evaluation of the necessary properties at the geometries visited during the dynamics simulations.

In order to account for nonadiabatic effects, i.e., transitions from one state to another, further approximations have to be introduced.[34] One method, which is frequently used to account for such transitions, is the surface-hopping method originally developed by Tully.[35] A popular extension for this method including not only nonadiabatic couplings (NACs) but also other couplings, e.g., spin–orbit couplings (SOCs), is the SHARC (surface hopping including arbitrary couplings) approach.[36−38] Importantly, NACs, also called derivative couplings, are used to determine the hopping directions and probabilities between states of the same spin multiplicity.[36,37,39−41] The NAC vector (denoted as CNAC) between two states, i and j, can be computed as[39,42,43]where the second-order derivatives are neglected. As a further difficulty, NACs are often missing from quantum chemistry implementations and hopping probabilities and directions are thus often approximated.[44−52] SOCs (denoted as CSOC) are present between states of different spin multiplicityand determine the rate of intersystem crossing. They are obtained as off-diagonal elements of the Hamiltonian matrix in standard electronic-structure calculations.[37,53]

Most of the recent studies involving ML dynamics deal with ground-state MD simulations, see e.g. refs (54−78), where one of the most promising ML models is SchNet,[79,80] a deep continuous-filter convolutional-layer neural network. In contrast to popular ML models like RuNNer,[64] n2p2,[81] TensorMol,[82] ANI,[63] or the DeePMD model[83] that require hand-crafted molecular descriptors, SchNet belongs to the class of message-passing neural networks.[84] Other examples of such networks are the DTNN,[85] PhysNet,[86] or HIP-NN.[87] The advantage is that the descriptors of molecules are automatically designed by a deep neural network and are based on the provided data set. In this way, the descriptors are tailored to the encountered chemical environments. Thus, we choose SchNet as a convenient platform for our developments.

An arising difficulty compared to ground-state energies and properties is that for the excited states not only one, but several PESs as well as the couplings between them have to be taken into account. Only a small but quickly increasing number of studies deal with the treatment of excited states and their properties using ML.[15−28,31,33,88]

In addition to the higher dimensionality that can be tackled with ML,[89,90] the learning of couplings proves challenging because properties resulting from electronic wave functions of two different states, Ψ and Ψ, have their sign dependent on the phase of the wave functions.[32,33,91,92] Because the wave function phase is not uniquely defined in quantum chemistry calculations, random phase jumps occur, leading to sign jumps of the coupling values along a reaction path. Hence, the couplings can not be learned directly as obtained from a quantum chemistry calculation. An option is to use a phase correction algorithm to preprocess data and remove these random phase jumps. Assuming that the effect of the Berry phase remains minor on the training set, smooth properties are obtained that are learnable by ML models.[32] However, this approach is expensive and many quantum chemistry reference computations are necessary to generate the training set. In cases of large polyatomic molecules with many close-lying energetic states, this approach might even be infeasible.

The aim of this Letter is to provide a framework to carry out efficient excited-state MLMD simulations and to combine two popular methods for this purpose: the SHARC approach for photodynamics with states of different multiplicity and SchNet to efficiently and accurately fit potential energies and other molecular properties. We call this combination the SchNarc approach and adapted SchNet for the treatment of excited-state potentials, their forces, and couplings for this purpose. The SchNarc approach can overcome the current limitations of existing MLMD simulations for excited states by allowing (i) a phase-free training to omit the costly preprocessing of raw quantum chemistry data and, to treat (ii) rotationally covariant NACs, which can either be trained or (iii) alternatively be approximated from only ML potentials, their gradients, and Hessians, and to treat (iv) SOCs. When using the phase-free training in combination with the approximated NACs, the costs required for the training set generation can be reduced substantially. Further, each data point can be computed in parallel, which is not possible when phase correction is needed. With all these methodological advances, the SchNarc ML approach simplifies accelerated nonadiabatic dynamics simulations, broadening the range of possible users and the scope of systems, in order to make long time scales accessible.

To validate our developments, the surface hopping dynamics of two molecules, showing slow and ultrafast excited-state dynamics, are investigated. The first molecule is CH2NH2+, of which we take a phase-corrected training set from ref (32). Using the same level of theory (MR-CISD(6,4)/aug-cc-pVDZ) with the program COLUMBUS,[93] the training set is recomputed without applying phase correction to train ML models also on raw data obtained directly from quantum chemistry programs. This training set should be used to validate our phase-free training approach. The ML models are trained on energies, gradients, and NACs for three singlet states using 3000 data points.

Slow photoinduced processes are present in thioformaldehyde, CSH2.[94] The training set is built up of 4703 data points with two singlet states and two triplet states after initial sampling of normal modes and adaptive sampling with simple multilayer feed-forward neural networks according to the scheme described in ref (32). This scheme applied for the training set generation is based on an uncertainty measure. Two (or more) ML models are trained and dynamics simulations are performed. At each time step, the predictions of the different ML models are compared. Whenever the error between the models exceeds a manually defined threshold, the molecular geometry visited at this time step is recomputed with quantum chemistry and the data added to the training set. Our previously proposed network[32] was used for this purpose because the training set was generated before SchNarc was developed. The sampling procedure is largely independent of the network architecture used and could therefore also be carried out with SchNarc but has not been tested here. The program MOLPRO[95] is used for the reference calculations with CASSCF(6,5)/def2-SVP.

The main novelty of the phase-free training is that it removes the influence of the arbitrary phase during the learning process of an ML model. It can be applied to any existing ML model capable of treating excited states and any existing data set for excited states eliminating an expensive and time-consuming phase correction preprocessing. The chosen ML model here is SchNet, for which the new loss function termed phase-less loss function is implemented and tested using the methylenimmonium cation, CH2NH2+. The phase-less loss is based on the standard L2 loss, but here, the squared error of the predicted properties is computed 2 times, with NS being the total number of states. The value of each property, LP (i.e., LSOC and LNAC), that enters the loss function is the minimum function of all possible squared errors εP:withfor vectorial and nonvectorial properties, respectively. The error εP for a specific phase is computed as the mean squared error of a property P from quantum chemistry (index QC) and machine learning (index ML). The property P couples different states, indicated by i and j. Because the wave function of each of the states can have an arbitrary phase, the property P that couples state i and j has to be multiplied with a product of the phases for these states, p·p. The phases for all states together form a vector p with entries of either +1 and −1. Which of the 2 possible combinations for p is chosen is indicated by the index k, also defined in eq . The possible combination that gives the lowest error enters the loss function. This is done for all samples inside of the training set and can be seen as an internal ML-based phase correction. Because of the continuity of ML fitted functions, the sign of properties is consistent for predictions. Note that the relative signs within one vector remain and must be predicted correctly for successful training.

The overall loss function in this work is a combination of such phase-less loss functions and mean squared errors for all properties with a trade-off factor to account for their relative magnitude. The relative magnitude of each property is defined by a manually set trade-off factor (a detailed description of the implementation is given in the Supporting Information).

Results are given in Figure , which shows the population schemes of CH2NH2+ obtained after excitation to the second excited singlet state, S2. The populations obtained from quantum chemistry are shown by dotted lines. SchNarc models are illustrated using solid lines. Panels A and C are obtained from SchNarc models trained on a data set that is not phase-corrected; that is, it contains couplings that can randomly switch their sign. Those are compared to populations obtained from models trained on phase-corrected data (panels B and D). As can be seen, the L2 loss function, as used in the upper plots, leads to an accurate ML model to reproduce ultrafast transitions only in the case of phase-corrected data (panel B), whereas this loss can not be used when trained on raw quantum chemistry data (panel A). In comparison, a SchNarc simulation with an ML model that applies the phase-less loss function is successful in reproducing the populations for both training sets (panels C and D).

Figure 1

Populations obtained from 90 QC (MR-CISD(6,4)/aug-cc-pVDZ) trajectories are shown by dotted lines and are compared to populations resulting from 1000 trajectories initially excited to the S2 obtained from SchNet (solid lines) that is trained on (A) not phase-corrected data and takes the L2 norm as loss function; (B) a similar SchNet model, but trained on phase-corrected data; (C) a SchNet model trained on not phase-corrected data, but using the new phase-less loss function; and (D) a SchNet model trained on phase-corrected data using the new phase-less loss function.

In those simulations, the NACs are multiplied with the corresponding energy gaps, i.e., C̃NAC = CNAC ·ΔE, to get rid of singularities.[21,33] These smooth couplings C̃NAC are not directly learned, but rather constructed as the derivative of a virtual property, analogously to forces that are predicted as derivatives of an energy-ML model. The virtual property is the multidimensional antiderivative of the rightmost expression in eq , (a derivation is given in the Supporting Information). Compared to previous ML models for NACs,[29,30,32,33] where NACs are learned and predicted as direct outputs or even single values, this approach provides rotational and translational covariance, which has recently been achieved in a similar way for the electronic friction tensor.[96]

However, even without the need of preprocessing the training set, the costly computations of NAC vectors for the training set generation remain. Approximations of NACs exist and often involve the computation of the squared energy-gap Hessian.[11,97−101] Their use in dynamics simulations is rather impracticable with quantum chemistry methods, especially in the case of complex systems, because of the expenses of computing second-order derivatives.

Here, we take advantage of the efficiency for second-order derivative computation from ML models with respect to atomic coordinates to obtain the Hessians of the fitted PESs:with R being the atomic coordinates of a molecular system. Note that Hessians are also employed in quantum dynamics simulations,[102,103] which might open further applications for our implementation.

The squared energy-gap Hessian can be further obtained as the sum of two symmetric dyads, which define the branching space.[101] Hence, this Hessian can be employed to obtain the symmetric dyad of the smooth NACs via[11,104]After singular value decomposition, the hopping direction can be computed as the eigenvector, v, of the largest nonzero eigenvalue[101,104,105] with the corresponding eigenvalue, λ, as the squared magnitude of the ML smooth coupling, C̃NAC. The final approximated NAC vectors, CaNAC, between two states are then

The approximated NAC vectors can be employed in the vicinity of a conical intersection; otherwise, the output becomes too noisy. For the latter reason, we define thresholds of 0.5 and 1.0 eV for the energy gaps to compute approximated NACs between coupled singlet–singlet states and triplet–triplet states, respectively. It is worth mentioning that the ML models slightly overestimate the energy gaps at a conical intersection,[32] because ML PESs, in contrast to quantum chemical PESs, are smooth everywhere and can reproduce the cones present in such critical regions only to a certain extent. In ref (105), approximated NACs were applied for a 1D system, and their usefulness in combination with ML was already anticipated.

We turn this idea into reality and show ML excited-state dynamics with approximated NACs for the methylenimmonium cation, CH2NH2+, as presented before, and thioformaldehyde, CSH2. A detailed analysis on the reference computations, the ab initio methods applied, as well as information on the timing of the Hessian evaluation are given in the Supporting Information in sections S2 and S3.1, respectively. The quality of the approximated NACs that is further compared to learned NACs is assessed using a linear vibronic coupling model of sulfur dioxide. The results are given in the Supporting Information in section S4.1 in Figures S1–S3 and support the validity of this approximation. Scatter plots and scans along a reaction coordinate of CH2NH2+ and CSH2 are further computed for energies, gradients, and couplings (and approximated NACs for CH2NH2+) in sections S4.2 (Figures S4 and S5) and S4.3 (Figures S6 and S7), respectively. Computed normal modes from ML Hessians are compared to reference values in Tables S3–S5 and show their accuracy. None of the data points from the ab initio MD simulations, to which we compare our SchNarc models, are included in the training sets, and thus, the dynamics simulations can be seen as an additional test.

The populations of CH2NH2+ are given in Figure . The reference population (panel A) is compared to SchNarc simulations trained on energies and gradients (panel B). As is visible, CH2NH2+ serves as a test system for ultrafast population transfer after photoexcitation. Transitions from the second excited singlet state back to the ground state take place within 100 fs,[32,106] which can be reproduced only with accurate NACs.[33] Those transitions can be reproduced with SchNarc using also the approximated NACs.

Figure 2

Quantum populations of the methylenimmonium cation obtained from (A) 90 trajectories using MR-CISD/aug-cc-pVDZ (QC) and (B) 1000 trajectories using SchNarc (ML). ML models are trained only on energies and gradients, and NACs are approximated from energies, gradients, and Hessians of ML models; this is in contrast to results from Figure , where ML models are also trained on NACs.

The application of the NAC approximation is further tested on CSH2, showing slow population transfer. This model also includes triplet states; thus, SOCs are additionally treated with SchNarc. To the best of our knowledge, for the first time, SOCs are trained with ML as directly obtained from quantum chemistry. The population curves are given in Figure , where panel A gives the reference population for 3000 fs and panel B the SchNarc populations. In contrast to the methylenimmonium cation, the CSH2 molecule serves as a test system for slow populations transfer and shows intersystem crossing strongly dependent on the accuracy of the underlying potentials.[94] Inaccurate ML models would thus be unable to reproduce the reference dynamics. Also the slow population can be reproduced accurately, which proves the validity of the ML approach.

Figure 3

Quantum populations of populations up to 3000 fs of the thioformaldehyde molecule obtained from (A) 100 trajectories using CASSCF(6,5)/def2-SVP (QC) and (B) 9590 trajectories using SchNarc (ML).

In summary, the SchNarc framework combines the SHARC[37] approach for surface hopping and the SchNet[80] approach for ML and introduces several methodological developments, which simplify the use of nonadiabatic ML dynamics substantially. SchNarc takes advantage of SchNet’s automatic generation of representations for the molecular structure and extends it to excited states. The training of ML models is facilitated by using the phase-less loss and the NAC approximation, avoiding quantum chemical NAC calculations at all. Thus, photodynamics simulations are possible based on solely ML PESs, their derivatives, and SOCs. Furthermore, this method allows for an efficient computation of the Hessians of all the excited states at each time step. Hence, SchNarc allows for efficient nonadiabatic dynamics simulations of excited states and light-induced processes including internal conversion and intersystem crossing.