Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

DOI: 10.1039/C9SC01742A
(Edge Article)
Chem. Sci., 2019, Advance Article

Julia Westermayr^{a},
Michael Gastegger^{b},
Maximilian F. S. J. Menger^{ac},
Sebastian Mai^{a},
Leticia González^{a} and
Philipp Marquetand*^{a}
^{a}Institute of Theoretical Chemistry, Faculty of Chemistry, University of Vienna, 1090 Vienna, Austria. E-mail: philipp.marquetand@univie.ac.at
^{b}Machine Learning Group, Technical University of Berlin, 10587 Berlin, Germany
^{c}Dipartimento di Chimica e Chimica Industriale, University of Pisa, Via G. Moruzzi 13, 56124 Pisa, Italy

Received
9th April 2019
, Accepted 2nd August 2019

First published on 5th August 2019

Photo-induced processes are fundamental in nature but accurate simulations of their dynamics are seriously limited by the cost of the underlying quantum chemical calculations, hampering their application for long time scales. Here we introduce a method based on machine learning to overcome this bottleneck and enable accurate photodynamics on nanosecond time scales, which are otherwise out of reach with contemporary approaches. Instead of expensive quantum chemistry during molecular dynamics simulations, we use deep neural networks to learn the relationship between a molecular geometry and its high-dimensional electronic properties. As an example, the time evolution of the methylenimmonium cation for one nanosecond is used to demonstrate that machine learning algorithms can outperform standard excited-state molecular dynamics approaches in their computational efficiency while delivering the same accuracy.

Computer simulations of photodynamics typically rely on molecular dynamics simulations of coupled nuclei and electrons. These simulations require the computation of high-dimensional potential energy surfaces (PESs), i.e., the electronic energy levels of the molecule for all possible molecular configurations, using quantum chemistry. The calculation of these PESs is usually the most expensive part of the dynamics simulations^{15} and therefore, different approximations are necessary and ubiquitous. For the electronic ground state, the time-consuming quantum chemical calculations are often replaced with force fields^{16} but no standard force fields are available to describe electronically excited states. Another drawback of most conventional force fields is their inability to describe the breaking and formation of chemical bonds. Recently, increasing effort has been devoted to ML potentials,^{17,18} where an accurate representation of the ground state PES including bond breaking^{19} and formation is promised.^{16,20–32} Similarly, modified Shepard interpolation is used to construct PESs in low-dimensional systems and adapt them in out-of-confidence regions.^{33,34} However, the problem of obtaining accurate full-dimensional PESs for excited states in order to simulate long time photodynamics has not been solved yet. A few studies have focused on the prediction of excited state dynamics as well as on excited-state properties such as spectral densities with ML.^{35–44} The breakdown of the Born–Oppenheimer approximation, leading to critical regions in the coupled excited state PESs,^{45} poses yet another obstacle to quantum chemistry (QC) and consequently also ML.^{39–41} Among those critical regions are conical intersections (or state crossings), where two PESs get into close proximity. The underlying elements that become important in such areas are nonadiabatic couplings (or spin–orbit couplings). They induce non-radiative transitions between two electronic states of the same (or different) spin-multiplicities involving ultrafast rearrangements of both nuclei and electrons. These challenges led to the need for intermittent quantum chemistry calculations^{39,40} or omittance of couplings between different PESs^{41} in ML driven photodynamics. Hence long time photodynamics are still lacking and the possibility to additionally represent the aforementioned nonadiabatic derivative couplings between PESs fundamental to model photodynamics has not been demonstrated yet. Here we overcome all these different bottlenecks using deep neural networks (NNs) and achieve the simulation of photodynamics for long time scales. We expand on the idea of using ML to obtain potentials for electronic excited states, as well as arbitrary couplings within a framework that combines ML with trajectory surface hopping molecular dynamics (Fig. 1). Our ML approach is fully capable of describing all necessary properties for executing nonadiabatic excited-state molecular dynamics on the order of nanoseconds. These properties include electronic energies, gradients, spin–orbit couplings, nonadiabatic couplings, and dipole moments of molecules. Additionally, the underlying potentials and couplings can be used to optimize critical points of the configurational space, such as potential minima or crossing points, which are important for interpreting photochemical mechanisms.

An ensemble of two NNs is used not only during the initial adaptive sampling period, but also for the production dynamics simulations in order to check the accuracy of the NN predictions and to discover undersampled regions of conformational space. After 10 ps, the threshold for the RMSE between NN forecasts is not reduced anymore but kept at the previous value when a new data point is added to the training set and NNs are retrained. More details on criteria for the thresholds and iterations are discussed in the ESI.†

Quantum chemical properties that were learned with NNs are energies, gradients, permanent as well as transition dipole moments, and NACs. Other quantities like spin–orbit couplings can also be trained (see the analytical model in the ESI†). Although the (transition) dipole moments are not needed for the present dynamics simulation, calculating them on-the-fly enables the computation of pump–probe schemes, static-field interactions, or time-resolved spectra (see for example ref. 52 and 53). While energies are directly used for training purposes in a single NN, forces are predicted as analytical derivatives of the NNs,^{54} ensuring energy conservation.^{24,32,39} Similarly, permanent dipole moments are directly used in the training. However, couplings (as well as transition dipole moments) need to be pre-processed as they are computed from the wave functions of two different electronic states and therefore depend on the relative phases of these two wave functions. Phase inconsistencies need to be eliminated in order to avoid ill-behaved photodynamics,^{55} as is described in the following subsection.

Such a global phase convention is not mathematically possible for general polyatomic molecules due to the existence of the so-called Berry (or geometric) phase.^{56} Due to the latter, the phase depends on the path between a given geometry and the reference point.^{57} Still the above phase correction is advantageous because it removes phase jumps from almost all parts of configurational space. This is critically necessary to make the data learnable. Only the non-removable phase jumps from the Berry phase remain, but occupy a small volume of configurational space. Hence, our phase correction is assumed to leave the dynamics mostly unaffected. For instance, successful surface hopping algorithms without phase tracking, such as the Zhu–Nakamura theory,^{58,59} exist and substantiate the validity of this approximation. In the case of the Zhu–Nakamura theory, dynamics are comparable to conventional surface-hopping molecular dynamics simulations propagated from NAC vectors.^{59–63} Note that the approximated phase correction for generation of the training set above cannot be circumvented by learning the absolute value of couplings since the relative sign between nonadiabatic coupling vectors of each atom in the x, y and z directions should be retained.

In order to make off-diagonal elements learnable for ML models, phases are tracked by computing wave function overlaps between adjacent molecular geometries.^{55,64,65} If the geometries are close enough, the overlaps will be sufficiently large and contain values close to +1 or −1, allowing a detection of phase changes. In cases where molecular geometries are too far apart, the overlap will generally be close to zero, offering no information about a phase change. In this case, we resort to interpolation between the two molecular geometries and iterative computation of wave function overlaps. In principle, the interpolation can be carried out between the new geometry and any geometry already inside the training set as long as the wavefunction of this previous geometry is stored. Storing the wavefunctions for at least a few geometries and identifying the most suitable one for interpolation via root mean square deviations of the geometry should be considered for larger and more flexible molecules.

Especially for large molecules, where many states lie close in energy, the so-called “intruder states” might become problematic. Such states are excluded at the reference geometry, but are included at another geometry due to an energy change, thus leading to small overlaps for the phase tracking algorithm. In such situations, different possibilities for adapting the phase correction algorithm should be considered. For instance, additional electronic states could be computed with QC. Those should not be included in the training data, but only used to continuously track the phase of all relevant states. This process then still stays affordable, since the additional states do not require a computation of gradients or couplings and do not have to be considered further. Additional details on the phase correction algorithm are given in Section S2 in the ESI.†

The reference quantum chemical computations were carried out with COLUMBUS^{66} using the accurate multi-reference configuration interaction method including single and double excitations and a double-zeta basis set (abbreviated to MR-CISD/aug-cc-pVDZ and in the following sections labelled as QC1). For comparison, we carried out quantum chemical computations with another basis set, 6-31++G**, using the same MR-CISD method (abbreviated to QC2 in the following sections). NNs were implemented in Python using the numpy^{67} and theano^{68} packages. They were trained on energies, forces, dipole moments and nonadiabatic couplings, obtained with the QC1 method using the adaptive sampling scheme described above, resulting in about 4000 data points (mean absolute error (MAE) energies among all states: ; MAE forces among all states: ; see also Tables S2, S4 and S5 in the ESI as well as Fig. S2† for analysis of different states). Using each method, QC1, QC2, and NNs trained on QC1, we simulated the dynamics of the methylenimmonium cation after excitation from the electronic ground state (S_{0}) to the second excited electronic state (S_{2}) over 100 fs using a time step of 0.5 fs.

Optimizations of minima were carried out with the SHARC tools that utilize an external optimizer, ORCA,^{69} where the computed energies and gradients^{70,71} from the NNs were fed in or those from COLUMBUS for comparison.

Fig. 3 Population dynamics of CH_{2}NH_{2}^{+} based on deep NNs and traditional quantum chemistry: comparison between results obtained from (A) QC1 (90 trajectories) and NN (3846 trajectories) and (B) QC1 (90 trajectories) and QC2 (88 trajectories). For completeness, the populations from 90 trajectories propagated with NNs are given in Fig. S5A in the ESI along with geometrical analysis along the trajectories in Fig. S5B.† |

One of the first advantages of the NN driven dynamics simulations is that due to their very low computational cost, a much larger number of trajectories (3846) was simulated than what is typically possible with standard quantum chemistry (90). This enlarged statistics provides smooth population curves for the NN simulations (a comparison of the curves with an identical number of trajectories for NNs and QC1 can be found in Fig. S5A in the ESI along with analysis of energy conservation in Table S11†).

In order to estimate the magnitude of the error obtained with the NNs, we carried out a second ab initio molecular dynamics study with an additional, very similar, quantum chemistry method where only the double-zeta basis set is changed from aug-cc-pVDZ to 6-31++G**. As Fig. 3B shows, the differences between the two levels of theory are of the same order of magnitude as those encountered between NNs and quantum chemistry, indicating that the agreement between the methods is very good. The MAE in population between QC1 and NNs is 0.057 and between QC1 and QC2 it is 0.099. Time constants derived from dynamics with each method also agree well. The time constant from S_{2} to S_{1} is 18.3 fs according to the QC1 method, which is comparable to the QC2 method with 25.0 fs and to NN driven dynamics with 25.2 fs. The time constant obtained for transitions from S_{1} to S_{0} is 51.0 fs for the QC1 method, which is very similar to the value obtained with NNs (52.6 fs), whereas the QC2 method yields a time constant of 73.2 fs.

After nonadiabatic dynamics using deep NNs has been validated for short time scales, we show the major advantage of the method, i.e. that it is able to overcome the problem of limited simulation time and predict long excited-state dynamics. Fig. 4 shows the population dynamics of the methylenimmonium cation on a logarithmic scale up to 1 nanosecond (ns), i.e., 10^{4} times longer than they were simulated using our quantum chemical reference method. Up to 10 ps, we simulated an ensemble of 200 trajectories with 2 NNs using the adaptive sampling scheme described above in order to correctly predict events not yet learned by the NNs. After that, 2 trajectories are propagated up to 1 ns for demonstration purposes using 2 NNs. The populations are thus averaged over 200 trajectories up to 10 ps and over 2 trajectories from 10 ps up to 1 ns, respectively. As can be seen, the molecule relaxes to the ground state after around 300 fs. Due to the remaining kinetic energy a few hops between different states are recorded and can be regarded as noise. A movie of one trajectory over 10 ps is part of the ESI Movie S2.†

The propagation of a CH_{2}NH_{2}^{+} trajectory for 10 ps can be executed in less than 6 hours on one core, which is 300 times faster than the calculation with the quantum chemical reference method. The propagation of 1 ns took 59 days employing two deep NNs serially, whereas an estimated 19 years of computation would have been required with the quantum chemical reference.

To this aim we optimize two minimum energy conical intersections in CH_{2}NH_{2}^{+}, one between the S_{2} state and the S_{1} state and another one between the S_{1} state and the S_{0} state. We use the QC1 method and NNs to perform potential energy scans around the minimum energy conical intersections optimized at the QC1 level of theory. As can be seen from Fig. 5A–D, typical curved seams of conical intersections between the S_{2} and S_{1} states (Fig. 5A (QC1) and 5B (NN)) and the S_{1} and S_{0} states (Fig. 5C (QC1) and 5D (NN)) are obtained around the minimum energy conical intersections.^{74} The NNs get the shape of this seam correct with slightly larger energy gaps between the crossing surfaces due to the fact that NN potentials need to be differentiable at any point. Analysis of 408 (for the S_{1}/S_{0} CI) and 302 (for the S_{2}/S_{1} CI) configurations around the minimum energy conical intersections – identified by an energy gap smaller than 0.8 eV according to the QC1 method – showed that on average, the gaps are overestimated by 0.068 eV for S_{1}/S_{0} and by 0.014 eV for S_{2}/S_{1} by our NNs. As can be seen from Fig. 5, the potentials around the S_{1}/S_{0} CI are flatter than the potentials around the S_{2}/S_{1} CI, indicating that hopping geometries are closer to the CI in the latter case and that the molecules can also hop farther from the CI in the former case.

Fig. 5 Potential energy scans around the minimum energy conical intersections obtained with QC1 of the S_{2} and S_{1} states (A and B) and S_{1} and S_{0} states (C and D). Panels A and C show the PESs calculated with QC1, while panels B and D illustrate NN potentials. See the caption of Fig. S7 in the ESI† for clarification of the dihedral angle. |

Fig. 6 shows the scatter plots of the optimized geometries of the minimum energy conical intersections projected along two important coordinates together with the hopping geometries and the geometries contained in the training set. As can be seen, the hopping geometries between the S_{2} and S_{1} states are mainly located close to the optimized geometry of the minimum energy conical intersection, while the hopping geometries in the case of the S_{1}/S_{0} crossing are more widely distributed around the optimized geometry. As a consequence, the S_{2}/S_{1} crossing is sampled more comprehensively, since more trajectories pass by near the minimum energy conical intersection. This observation also explains the larger NN energy gap obtained for the second crossing, the S_{1}/S_{0} CI, in Fig. 5.

Fig. 6 Scatter plots showing the distribution of hopping geometries obtained with QC1, QC2, and NN as well as optimized S_{1}/S_{0} (A) and S_{2}/S_{1} (B) minimum energy conical intersections (CIs) along with the geometries that make up the training set with 4000 data points. The actual geometry is depicted on top (geometrical parameters are given in Fig. S7B†). A zoom-in of the regions near the optimized points is shown in Fig. S7A in the ESI† together with a definition of the dihedral and pyramidalization angles. |

The optimizations of the minimum energy conical intersections were independently performed with the trained NN, as well as with the QC1 and QC2 methods for comparison. The optimized molecular geometries (shown in Fig. S6 along with Cartesian coordinates in the ESI†) agree well. As can be seen, the driving force for the transition from the S_{2} state to the S_{1} state is an elongation of the C–N bond in combination with a bipyramidalization. The torsion of the molecule further leads to internal conversion to the ground state, S_{0}. Additionally, each method results in a comparable distribution of hopping geometries around the optimized points, which in practice is of utmost importance^{75} for describing the population transfer in the simulations correctly. There are very few NN hopping geometries at either large pyramidalization angles (S_{1}/S_{0} CI) or long C–N bonds (S_{2}/S_{1} CI), compared to the QC trajectories. This finding correlates with the distribution of training set geometries, which are also absent in these regions of the PES (see the grey circles in Fig. 6). Configurations obtained via sampling of normal modes are clearly visible by a dense alignment of data points. However, the configurations obtained via adaptive sampling are mostly centered in the middle of the plot for the S_{1}/S_{0} CI and close to the optimized CI for the S_{2}/S_{1} crossing, explaining the smaller distribution of NN hopping geometries. Further analysis showed that geometries at large bond lengths are approximately 4 eV higher in energy than the geometries close to the optimized minimum energy conical intersection in the case of the S_{2}/S_{1} crossing. Therefore, trajectories simulated during adaptive sampling probably did not visit those regions of the PES. In the case of the S_{1}/S_{0} crossing, this effect is less pronounced and the geometries with a large pyramidalization angle are approximately 1–1.5 eV larger in energy than the configurations close to the optimized CI, indicating again the much flatter potential.

- I. Goodfellow, Y. Bengio and A. Courville, Deep Learning, MIT Press, 2016 Search PubMed.
- D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. van den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S. Dieleman, D. Grewe, J. Nham, N. Kalchbrenner, I. Sutskever, T. Lillicrap, M. Leach, K. Kavukcuoglu, T. Graepel and D. Hassabis, Nature, 2016, 529, 484–489 CrossRef CAS PubMed.
- K. Bansak, J. Ferwerda, J. Hainmueller, A. Dillon, D. Hangartner, D. Lawrence and J. Weinstein, Science, 2018, 359, 325–329 CrossRef CAS PubMed.
- B. Sanchez-Lengeling and A. Aspuru-Guzik, Science, 2018, 361, 360–365 CrossRef CAS PubMed.
- K. T. Butler, D. W. Davies, H. Cartwright, O. Isayev and A. Walsh, Nature, 2018, 559, 547–555 CrossRef CAS PubMed.
- B. R. Goldsmith, J. Esterhuizen, J.-X. Liu, C. J. Bartel and C. Sutton, AIChE J., 2018, 64, 2311–2323 CrossRef CAS.
- G. Cerullo, D. Polli, G. Lanzani, S. De Silvestri, H. Hashimoto and R. J. Cogdell, Science, 2002, 298, 2395–2398 CrossRef CAS PubMed.
- T. Schultz, E. Samoylova, W. Radloff, I. V. Hertel, A. L. Sobolewski and W. Domcke, Science, 2004, 306, 1765–1768 CrossRef CAS PubMed.
- W. J. Schreier, T. E. Schrader, F. O. Koller, P. Gilch, C. E. Crespo-Hernández, V. N. Swaminathan, T. Charell, W. Zinth and B. Kohler, Science, 2007, 315, 625–629 CrossRef CAS PubMed.
- C. Rauer, J. J. Nogueira, P. Marquetand and L. González, J. Am. Chem. Soc., 2016, 138, 15911–15916 CrossRef CAS PubMed.
- E. Romero, V. I. Novoderezhkin and R. v. Grondelle, Nature, 2017, 543, 355 CrossRef CAS PubMed.
- I. Ahmad, S. Ahmed, Z. Anwar, M. A. Sheraz and M. Sikorski, Int. J. Photoenergy, 2016, 2016, 1–19 CrossRef.
- S. Mathew, A. Yella, P. Gao, R. Humphry-Baker, B. F. E. Curchod, N. Ashari-Astani, I. Tavernelli, U. Rothlisberger, M. K. Nazeeruddin and M. Grätzel, Nat. Chem., 2014, 6, 242 CrossRef CAS PubMed.
- A. P. Bartók, S. De, C. Poelking, N. Bernstein, J. R. Kermode, G. Csányi and M. Ceriotti, Sci. Adv., 2017, 3, e1701816 CrossRef PubMed.
- S. Mai, P. Marquetand and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1370 Search PubMed.
- S. Chmiela, H. E. Sauceda, K.-R. Müller and A. Tkatchenko, Nat. Commun., 2018, 9, 3887 CrossRef PubMed.
- M. Rupp, Int. J. Quantum Chem., 2015, 115, 1058–1073 CrossRef CAS.
- O. A. von Lilienfeld, Angew. Chem., Int. Ed., 2018, 57, 4164–4169 CrossRef CAS PubMed.
- F. Häse, I. F. Galván, A. Aspuru-Guzik, R. Lindh and M. Vacher, Chem. Sci., 2019, 10, 2298–2307 RSC.
- A. P. Bartók, M. C. Payne, R. Kondor and G. Csányi, Phys. Rev. Lett., 2010, 104, 136403 CrossRef PubMed.
- Z. Li, J. R. Kermode and A. De Vita, Phys. Rev. Lett., 2015, 114, 096405 CrossRef PubMed.
- M. Rupp, R. Ramakrishnan and O. A. von Lilienfeld, J. Phys. Chem. Lett., 2015, 6, 3309–3313 CrossRef CAS.
- J. Behler, J. Chem. Phys., 2016, 145, 170901 CrossRef PubMed.
- M. Gastegger, J. Behler and P. Marquetand, Chem. Sci., 2017, 8, 6924–6935 RSC.
- V. Botu, R. Batra, J. Chapman and R. Ramprasad, J. Phys. Chem. C, 2017, 121, 511–522 CrossRef CAS.
- J. S. Smith, O. Isayev and A. E. Roitberg, Chem. Sci., 2017, 8, 3192–3203 RSC.
- J. Behler, Angew. Chem., Int. Ed., 2017, 56, 12828–12840 CrossRef CAS PubMed.
- H. Zong, G. Pilania, X. Ding, G. J. Ackland and T. Lookman, npj Comput. Mater., 2018, 4, 48 CrossRef.
- A. P. Bartók, J. Kermode, N. Bernstein and G. Csányi, Phys. Rev. X, 2018, 8, 041048 Search PubMed.
- R. Xia and S. Kais, Nat. Commun., 2018, 9, 4195 CrossRef PubMed.
- H. Chan, B. Narayanan, M. J. Cherukara, F. G. Sen, K. Sasikumar, S. K. Gray, M. K. Y. Chan and S. K. R. S. Sankaranarayanan, J. Phys. Chem. C, 2019, 123, 6941 CrossRef CAS.
- A. S. Christensen, F. A. Faber and O. A. von Lilienfeld, J. Chem. Phys., 2019, 150, 064105 CrossRef PubMed.
- H. M. Netzloff, M. A. Collins and M. S. Gordon, J. Chem. Phys., 2006, 124, 154104 CrossRef PubMed.
- R. P. A. Bettens and M. A. Collins, J. Chem. Phys., 1999, 111, 816–826 CrossRef CAS.
- J. Behler, K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 115421 CrossRef.
- C. Carbogno, J. Behler, K. Reuter and A. Groß, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 035410 CrossRef.
- F. Häse, S. Valleau, E. Pyzer-Knapp and A. Aspuru-Guzik, Chem. Sci., 2016, 7, 5139–5147 RSC.
- F. Liu, L. Du, D. Zhang and J. Gao, Sci. Rep., 2017, 7, 1–12 CrossRef PubMed.
- D. Hu, Y. Xie, X. Li, L. Li and Z. Lan, J. Phys. Chem. Lett., 2018, 9, 2725–2732 CrossRef CAS PubMed.
- P. O. Dral, M. Barbatti and W. Thiel, J. Phys. Chem. Lett., 2018, 9, 5660–5663 CrossRef CAS PubMed.
- W.-K. Chen, X.-Y. Liu, W.-H. Fang, P. O. Dral and G. Cui, J. Phys. Chem. Lett., 2018, 9, 6702–6708 CrossRef CAS PubMed.
- D. M. G. Williams and W. Eisfeld, J. Chem. Phys., 2018, 149, 204106 CrossRef PubMed.
- C. Xie, X. Zhu, D. R. Yarkony and H. Guo, J. Chem. Phys., 2018, 149, 144107 CrossRef PubMed.
- Y. Guan, D. H. Zhang, H. Guo and D. R. Yarkony, Phys. Chem. Chem. Phys., 2019, 21, 14205–14213 RSC.
- W. Domcke, D. R. Yarkony and H. Köppel, Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, WORLD SCIENTIFIC, 2004 Search PubMed.
- M. Richter, P. Marquetand, J. González-Vázquez, I. Sola and L. González, J. Chem. Theory Comput., 2011, 7, 1253–1258 CrossRef CAS PubMed.
- J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071 CrossRef CAS.
- S. Mai, M. Richter, M. Ruckenbauer, M. Oppel, P. Marquetand and L. González, SHARC2.0: Surface Hopping Including Arbitrary Couplings – Program Package for Non-Adiabatic Dynamics, sharc-md.org, 2018 Search PubMed.
- J. Behler, Int. J. Quantum Chem., 2015, 115, 1032–1050 CrossRef CAS.
- V. Botu and R. Ramprasad, Int. J. Quantum Chem., 2015, 115, 1074–1083 CrossRef CAS.
- J. S. Smith, B. Nebgen, N. Lubbers, O. Isayev and A. E. Roitberg, J. Chem. Phys., 2018, 148, 241733 CrossRef PubMed.
- P. Marquetand, A. Materny, N. E. Henriksen and V. Engel, J. Chem. Phys., 2004, 120, 5871–5874 CrossRef CAS PubMed.
- F. P. Bonafé, F. J. Hernández, B. Aradi, T. Frauenheim and C. G. Sánchez, J. Phys. Chem. Lett., 2018, 9, 4355–4359 CrossRef PubMed.
- M. Gastegger and P. Marquetand, J. Chem. Theory Comput., 2015, 11, 2187–2198 CrossRef CAS PubMed.
- A. V. Akimov, J. Phys. Chem. Lett., 2018, 9, 6096–6102 CrossRef CAS PubMed.
- I. G. Ryabinkin, L. Joubert-Doriol and A. F. Izmaylov, Acc. Chem. Res., 2017, 50, 1785–1793 CrossRef CAS PubMed.
- S. Matsika and P. Krause, Annu. Rev. Phys. Chem., 2011, 62, 621–643 CrossRef CAS PubMed.
- P. Oloyede, G. Mil'nikov and H. Nakamura, J. Chem. Phys., 2006, 124, 144110 CrossRef PubMed.
- T. Ishida, S. Nanbu and H. Nakamura, Int. Rev. Phys. Chem., 2017, 36, 229–286 Search PubMed.
- C. Zhu, H. Kamisaka and H. Nakamura, J. Chem. Phys., 2002, 116, 3234–3247 CrossRef CAS.
- T. Ishida, S. Nanbu and H. Nakamura, J. Phys. Chem. A, 2009, 113, 4356–4366 CrossRef CAS PubMed.
- A.-H. Gao, B. Li, P.-Y. Zhang and K.-L. Han, J. Chem. Phys., 2012, 137, 204305 CrossRef PubMed.
- L. Yu, C. Xu, Y. Lei, C. Zhu and Z. Wen, Phys. Chem. Chem. Phys., 2014, 16, 25883–25895 RSC.
- S. Mai, P. Marquetand and L. González, Int. J. Quantum Chem., 2015, 115, 1215–1231 CrossRef CAS.
- F. Plasser, M. Ruckenbauer, S. Mai, M. Oppel, P. Marquetand and L. González, J. Chem. Theory Comput., 2016, 12, 1207 CrossRef CAS PubMed.
- H. Lischka, R. Shepard, R. M. Pitzer, I. Shavitt, M. Dallos, T. Müller, P. G. Szalay, M. Seth, G. S. Kedziora, S. Yabushita and Z. Zhang, Phys. Chem. Chem. Phys., 2001, 3, 664–673 RSC.
- S. van der Walt, S. C. Colbert and G. Varoquaux, Comput. Sci. Eng., 2011, 13, 22–30 Search PubMed.
- Theano Development Team, arXiv 1605.02688 [cs.SC], 2016.
- F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78 CAS.
- B. G. Levine, J. D. Coe and T. J. Martinez, J. Phys. Chem. B, 2007, 112, 405–413 CrossRef PubMed.
- M. J. Bearpark, M. A. Robb and H. B. Schlegel, Chem. Phys. Lett., 1994, 223, 269 CrossRef CAS.
- M. Barbatti, A. J. A. Aquino and H. Lischka, Mol. Phys., 2006, 104, 1053–1060 CrossRef CAS.
- J. Herbst, K. Heyne and R. Diller, Science, 2002, 297, 822–825 CrossRef CAS PubMed.
- D. R. Yarkony, J. Chem. Phys., 2005, 123, 204101 CrossRef PubMed.
- H. R. Hudock, B. G. Levine, A. L. Thompson, H. Satzger, D. Townsend, N. Gador, S. Ullrich, A. Stolow and T. J. Martínez, J. Phys. Chem. A, 2007, 111, 8500–8508 CrossRef CAS PubMed.

## Footnote |

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9sc01742a |

This journal is © The Royal Society of Chemistry 2019 |