Excited state, non-adiabatic dynamics of large photoswitchable molecules using a chemically transferable machine learning potential

Simon Axelrod

Department of Chemistry and Chemical Biology,
Harvard University, Cambridge, MA, 02138
Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA, 02139

Eugene Shakhnovich

Department of Chemistry and Chemical Biology,
Harvard University, Cambridge, MA, 02138

Rafael Gómez-Bombarelli
[Corresponding author<EMAIL_ADDRESS>Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA, 02139

Light-induced chemical processes are ubiquitous in nature and have widespread technological applications. For example, photoisomerization can allow a drug with a photo-switchable scaffold such as azobenzene to be activated with light. In principle, photoswitches with desired photophysical properties like high isomerization quantum yields can be identified through virtual screening with reactive simulations. In practice, these simulations are rarely used for screening, since they require hundreds of trajectories and expensive quantum chemical methods to account for non-adiabatic excited state effects. Here we introduce a diabatic artificial neural network (DANN) based on diabatic states to accelerate such simulations for azobenzene derivatives. The network is six orders of magnitude faster than the quantum chemistry method used for training. DANN is transferable to azobenzene molecules outside the training set, predicting quantum yields for unseen species that are correlated with experiment. We use the model to virtually screen 3,100 hypothetical molecules, and identify novel species with extremely high predicted quantum yields. The model predictions are confirmed using high-accuracy non-adiabatic dynamics. Our results pave the way for fast and accurate virtual screening of photoactive compounds.

[Evans et al., 2013]
Rachel C Evans, Peter Douglas, and Hugh D Burrow.
Applied photochemistry.
Springer, 2013.

[Kolpak and Grossman, 2011]
Alexie M Kolpak and Jeffrey C Grossman.
Azobenzene-functionalized carbon nanotubes as high-energy density
solar thermal fuels.
Nano letters, 110 (8):0 3156–3162, 2011.

[Mai and González, 2020]
Sebastian Mai and Leticia González.
Molecular photochemistry: Recent developments in theory.
Angewandte Chemie International Edition, 590
(39):0 16832–16846, 2020.

[Broichhagen et al., 2015]
Johannes Broichhagen, James Allen Frank, and Dirk Trauner.
A roadmap to success in photopharmacology.
Accounts of chemical research, 480 (7):0
1947–1960, 2015.

[Lerch et al., 2016]
Michael M Lerch, Mickel J Hansen, Gooitzen M van Dam, Wiktor Szymanski, and
Ben L Feringa.
Emerging targets in photopharmacology.
Angewandte Chemie International Edition, 550
(37):0 10978–10999, 2016.

[Mai and González, 2019]
Sebastian Mai and Leticia González.
Unconventional two-step spin relaxation dynamics of [Re (CO) 3
(im)(phen)]+ in aqueous solution.
Chemical science, 100 (44):0 10405–10411,

[Yu et al., 2020]
Jimmy K Yu, Christoph Bannwarth, Ruibin Liang, Edward G Hohenstein, and Todd J
Nonadiabatic dynamics simulation of the wavelength-dependent
photochemistry of azobenzene excited to the n$\pi^*$ and $\pi\pi^*$ excited
Journal of the American Chemical Society, 1420
(49):0 20680–20690, 2020a.

[Bannwarth et al., 2020]
Christoph Bannwarth, Jimmy K Yu, Edward G Hohenstein, and Todd J Martínez.
Hole–hole tamm–dancoff-approximated density functional theory: A
highly efficient electronic structure method incorporating dynamic and static
The Journal of Chemical Physics, 1530 (2):0
024110, 2020.

[Tully, 1998]
JohnáC Tully.
Mixed quantum–classical dynamics.
Faraday Discussions, 110:0 407–419, 1998.

[Shalashilin, 2009]
Dmitrii V Shalashilin.
Quantum mechanics with the basis set guided by ehrenfest
trajectories: Theory and application to spin-boson model.
The Journal of chemical physics, 1300 (24):0
244101, 2009.

[Ben-Nun et al., 2000]
Michal Ben-Nun, Jason Quenneville, and Todd J Martínez.
Ab initio multiple spawning: Photochemistry from first
principles quantum molecular dynamics.
The Journal of Physical Chemistry A, 1040
(22):0 5161–5175, 2000.

[Shiozaki et al., 2011]
Toru Shiozaki, Werner Győrffy, Paolo Celani, and Hans-Joachim Werner.
Communication: Extended multi-state complete active space
second-order perturbation theory: Energy and nuclear gradients.
Journal of Chemical Physics, 1350 (8):0
081106–081106, 2011.

[Olsen et al., 1988]
Jeppe Olsen, Björn O Roos, Poul Jørgensen, and Hans Jørgen Aa.
Determinant based configuration interaction algorithms for complete
and restricted configuration interaction spaces.
The Journal of chemical physics, 890 (4):0
2185–2192, 1988.

[Ma et al., 2011]
Dongxia Ma, Giovanni Li Manni, and Laura Gagliardi.
The generalized active space concept in multiconfigurational
self-consistent field methods.
The Journal of chemical physics, 1350 (4):0
044128, 2011.

[Li Manni et al., 2011]
Giovanni Li Manni, Francesco Aquilante, and Laura Gagliardi.
Strong correlation treated via effective hamiltonians and
perturbation theory.
The Journal of chemical physics, 1340 (3):0
034114, 2011.

[Ivanic, 2003]
Joseph Ivanic.
Direct configuration interaction and multiconfigurational
self-consistent-field method for multiple active spaces with variable
occupations. I. Method.
The Journal of chemical physics, 1190 (18):0
9364–9376, 2003.

[Song and Martínez, 2020]
Chenchen Song and Todd J Martínez.
Reduced scaling extended multi-state caspt2 (xms-caspt2) using
supporting subspaces and tensor hyper-contraction.
The Journal of chemical physics, 1520 (23):0
234113, 2020.

[Song et al., 2021]
Chenchen Song, Jeffrey B Neaton, and Todd J Martínez.
Reduced scaling formulation of caspt2 analytical gradients using the
supporting subspace method.
The Journal of Chemical Physics, 1540 (1):0
014103, 2021.

[Seritan et al., 2020]
Stefan Seritan, Christoph Bannwarth, B. Scott Fales, Edward G. Hohenstein, Sara
I. L. Kokkila-Schumacher, Nathan Luehr, James W. Snyder Jr., Chenchen Song,
Alexey V. Titov, Ivan S. Ufimtsev, and Todd J. Martínez.
TeraChem: Accelerating electronic structure and ab initio molecular
dynamics with graphical processing units.
The Journal of chemical physics, 1520 (22):0
224110, 2020.

[Marti and Reiher, 2011]
Konrad H Marti and Markus Reiher.
New electron correlation theories for transition metal chemistry.
Physical Chemistry Chemical Physics, 130
(15):0 6750–6759, 2011.

[Sharma and Chan, 2012]
Sandeep Sharma and Garnet Kin-Lic Chan.
Spin-adapted density matrix renormalization group algorithms for
quantum chemistry.
The Journal of chemical physics, 1360 (12):0
124121, 2012.

[Hermes and Gagliardi, 2019]
Matthew R Hermes and Laura Gagliardi.
Multiconfigurational self-consistent field theory with density
matrix embedding: The localized active space self-consistent field method.
Journal of chemical theory and computation, 150
(2):0 972–986, 2019.

[Marian et al., 2019]
Christel M Marian, Adrian Heil, and Martin Kleinschmidt.
The DFT/MRCI Method.
Wiley Interdisciplinary Reviews: Computational Molecular
Science, 90 (2):0 e1394, 2019.

[Li Manni et al., 2014]
Giovanni Li Manni, Rebecca K Carlson, Sijie Luo, Dongxia Ma, Jeppe Olsen,
Donald G Truhlar, and Laura Gagliardi.
Multiconfiguration pair-density functional theory.
Journal of chemical theory and computation, 100
(9):0 3669–3680, 2014.

[Gagliardi et al., 2017]
Laura Gagliardi, Donald G Truhlar, Giovanni Li Manni, Rebecca K Carlson, Chad E
Hoyer, and Junwei Lucas Bao.
Multiconfiguration pair-density functional theory: A new way to
treat strongly correlated systems.
Accounts of chemical research, 500 (1):0
66–73, 2017.

[Stein and Reiher, 2016]
Christopher J. Stein and Markus Reiher.
Automated selection of active orbital spaces.
Journal of chemical theory and computation, 120
(4):0 1760, 2016.

[Stein and Reiher, 2019]
Christopher J. Stein and Markus Reiher.
autoCAS: A program for fully automated multiconfigurational
Journal of computational chemistry, 400 (25):0
2216, 2019.

[Yu et al., 2020]
Jimmy K Yu, Christoph Bannwarth, Edward G Hohenstein, and Todd J Martínez.
Ab initio nonadiabatic molecular dynamics with hole–hole
Tamm–Dancoff approximated density functional theory.
Journal of Chemical Theory and Computation, 160
(9):0 5499–5511, 2020b.

[Filatov, 2015]
Michael Filatov.
Spin-restricted ensemble-referenced kohn–sham method: basic
principles and application to strongly correlated ground and excited states
of molecules.
Wiley Interdisciplinary Reviews: Computational Molecular
Science, 50 (1):0 146–167, 2015.

[Lee et al., 2018]
Seunghoon Lee, Michael Filatov, Sangyoub Lee, and Cheol Ho Choi.
Eliminating spin-contamination of spin-flip time dependent density
functional theory within linear response formalism by the use of zeroth-order
mixed-reference (mr) reduced density matrix.
The Journal of chemical physics, 1490 (10):0
104101, 2018.

[Cusati et al., 2012]
Teresa Cusati, Giovanni Granucci, Emilio Martínez-Núnez, Francesca
Martini, Maurizio Persico, and Saulo Vázquez.
Semiempirical Hamiltonian for simulation of azobenzene
The Journal of Physical Chemistry A, 1160
(1):0 98–110, 2012.

[Inamori et al., 2020]
Mayu Inamori, Takeshi Yoshikawa, Yasuhiro Ikabata, Yoshifumi Nishimura, and
Hiromi Nakai.
Spin-flip approach within time-dependent density functional
tight-binding method: Theory and applications.
Journal of computational chemistry, 410 (16):0
1538–1548, 2020.

[de Wergifosse et al., 2019]
Marc de Wergifosse, Christoph Bannwarth, and Stefan Grimme.
A simplified spin-flip time-dependent density functional theory
approach for the electronic excitation spectra of very large diradicals.
The Journal of Physical Chemistry A, 1230
(27):0 5815–5825, 2019.

[Qiao et al., 2020]
Zhuoran Qiao, Matthew Welborn, Animashree Anandkumar, Frederick R Manby, and
Thomas F Miller III.
OrbNet: Deep learning for quantum chemistry using symmetry-adapted
atomic-orbital features.
The Journal of Chemical Physics, 1530 (12):0
124111, 2020a.

[Schütt et al., 2021]
Kristof T. Schütt, Oliver T. Unke, and Michael Gastegger.
Equivariant message passing for the prediction of tensorial
properties and molecular spectra.
arXiv preprint arXiv:2102.03150, 2021.

[Unke et al., 2021]
Oliver T Unke, Stefan Chmiela, Michael Gastegger, Kristof T Schütt,
Huziel E Sauceda, and Klaus-Robert Müller.
SpookyNet: Learning force fields with electronic degrees of freedom
and nonlocal effects.
arXiv preprint arXiv:2105.00304, 2021.

[Ang et al., 2021]
Shi Jun Ang, Wujie Wang, Daniel Schwalbe-Koda, Simon Axelrod, and Rafael
Active learning accelerates ab initio molecular dynamics on reactive
energy surfaces.
Chem, 70 (3):0 738, 2021.

[Wang et al., 2020]
Wujie Wang, Tzuhsiung Yang, William H. Harris, and Rafael Gómez-Bombarelli.
Active learning and neural network potentials accelerate molecular
screening of ether-based solvate ionic liquids.
Chemical Communications, 560 (63):0 8920,

[Chen et al., 2018]
Wen-Kai Chen, Xiang-Yang Liu, Wei-Hai Fang, Pavlo O Dral, and Ganglong Cui.
Deep learning for nonadiabatic excited-state dynamics.
The journal of physical chemistry letters, 90
(23):0 6702–6708, 2018.

[Dral et al., 2018]
Pavlo O Dral, Mario Barbatti, and Walter Thiel.
Nonadiabatic excited-state dynamics with machine learning.
The journal of physical chemistry letters, 90
(19):0 5660–5663, 2018.

[Hu et al., 2018]
Deping Hu, Yu Xie, Xusong Li, Lingyue Li, and Zhenggang Lan.
Inclusion of machine learning kernel ridge regression potential
energy surfaces in on-the-fly nonadiabatic molecular dynamics simulation.
The journal of physical chemistry letters, 90
(11):0 2725–2732, 2018.

[Li et al., 2021]
Jingbai Li, Patrick Reiser, Benjamin R Boswell, André Eberhard, Noah Z
Burns, Pascal Friederich, and Steven A Lopez.
Automatic discovery of photoisomerization mechanisms with nanosecond
machine learning photodynamics simulations.
Chemical Science, 2021.

[Westermayr et al., 2019]
Julia Westermayr, Michael Gastegger, Maximilian FSJ Menger, Sebastian Mai,
Leticia González, and Philipp Marquetand.
Machine learning enables long time scale molecular photodynamics
Chemical science, 100 (35):0 8100–8107, 2019.

[Westermayr et al., 2020]
Julia Westermayr, Michael Gastegger, and Philipp Marquetand.
Combining SchNet and SHARC: The SchNarc machine learning approach
for excited-state dynamics.
The journal of physical chemistry letters, 110
(10):0 3828–3834, 2020.

[Westermayr and Marquetand, 2020]
Julia Westermayr and Philipp Marquetand.
Machine learning for electronically excited states of molecules.
Chemical Reviews, 2020a.

[Makhov et al., 2014]
Dmitry V Makhov, William J Glover, Todd J Martinez, and Dmitrii V Shalashilin.
Ab initio multiple cloning algorithm for quantum nonadiabatic
molecular dynamics.
The Journal of chemical physics, 1410 (5):0
054110, 2014.

[Zhu et al., 2004]
Chaoyuan Zhu, Shikha Nangia, Ahren W Jasper, and Donald G Truhlar.
Coherent switching with decay of mixing: an improved treatment of
electronic coherence for non-born–oppenheimer trajectories.
The Journal of chemical physics, 1210 (16):0
7658–7670, 2004.

[Richings et al., 2015]
GW Richings, Iakov Polyak, KE Spinlove, GA Worth, Irene Burghardt, and Benjamin
Quantum dynamics simulations using gaussian wavepackets: the vmcg
International Reviews in Physical Chemistry, 340
(2):0 269–308, 2015.

[Abedi et al., 2010]
Ali Abedi, Neepa T Maitra, and Eberhard KU Gross.
Exact factorization of the time-dependent electron-nuclear wave
Physical review letters, 1050 (12):0 123002,

[Abedi et al., 2014]
Ali Abedi, Federica Agostini, and EKU Gross.
Mixed quantum-classical dynamics from the exact decomposition of
electron-nuclear motion.
EPL (Europhysics Letters), 1060 (3):0 33001,

[Min et al., 2017]
Seung Kyu Min, Federica Agostini, Ivano Tavernelli, and Eberhard KU Gross.
Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous
approach to quantum (de) coherence.
The journal of physical chemistry letters, 80
(13):0 3048–3055, 2017.

[Curchod and Agostini, 2017]
Basile FE Curchod and Federica Agostini.
On the dynamics through a conical intersection.
The journal of physical chemistry letters, 80
(4):0 831–837, 2017.

[Ha et al., 2018]
Jong-Kwon Ha, In Seong Lee, and Seung Kyu Min.
Surface hopping dynamics beyond nonadiabatic couplings for quantum
The journal of physical chemistry letters, 90
(5):0 1097–1104, 2018.

[Beck et al., 2000]
Michael H Beck, Andreas Jäckle, Graham A Worth, and H-D Meyer.
The multiconfiguration time-dependent Hartree (MCTDH) method: a
highly efficient algorithm for propagating wavepackets.
Physics reports, 3240 (1):0 1–105, 2000.

[Wang and Thoss, 2003]
Haobin Wang and Michael Thoss.
Multilayer formulation of the multiconfiguration time-dependent
Hartree theory.
The Journal of chemical physics, 1190 (3):0
1289–1299, 2003.

[Burghardt et al., 1999]
Irene Burghardt, H-D Meyer, and LS Cederbaum.
Approaches to the approximate treatment of complex molecular systems
by the multiconfiguration time-dependent hartree method.
The Journal of chemical physics, 1110 (7):0
2927–2939, 1999.

[Tully, 1990]
John C Tully.
Molecular dynamics with electronic transitions.
The Journal of Chemical Physics, 930 (2):0
1061–1071, 1990.

[Zhu and Nakamura, 1992]
Chaoyuan Zhu and Hiroki Nakamura.
The two-state linear curve crossing problems revisited. ii.
analytical approximations for the stokes constant and scattering matrix: The
landau–zener case.
The Journal of chemical physics, 970 (11):0
8497–8514, 1992.

[Zhu and Nakamura, 1993]
Chaoyuan Zhu and Hiroki Nakamura.
The two-state linear curve crossing problems revisited. iii.
analytical approximations for stokes constant and scattering matrix:
Nonadiabatic tunneling case.
The Journal of chemical physics, 980 (8):0
6208–6222, 1993.

[Shu and Truhlar, 2020]
Yinan Shu and Donald G Truhlar.
Diabatization by machine intelligence.
Journal of Chemical Theory and Computation, 160
(10):0 6456–6464, 2020.

[Williams and Eisfeld, 2018]
David MG Williams and Wolfgang Eisfeld.
Neural network diabatization: A new ansatz for accurate
high-dimensional coupled potential energy surfaces.
The Journal of chemical physics, 1490 (20):0
204106, 2018.

[Guan et al., 2019]
Yafu Guan, Dong H Zhang, Hua Guo, and David R Yarkony.
Representation of coupled adiabatic potential energy surfaces using
neural network based quasi-diabatic hamiltonians: 1, 2 2 a' states of lifh.
Physical Chemistry Chemical Physics, 210
(26):0 14205–14213, 2019.

[Shao et al., 2003]
Yihan Shao, Martin Head-Gordon, and Anna I Krylov.
The spin–flip approach within time-dependent density functional
theory: Theory and applications to diradicals.
The Journal of chemical physics, 1180 (11):0
4807–4818, 2003.

[Francl et al., 1982]
Michelle M Francl, William J Pietro, Warren J Hehre, J Stephen Binkley, Mark S
Gordon, Douglas J DeFrees, and John A Pople.
Self-consistent molecular orbital methods. XXIII. A
polarization-type basis set for second-row elements.
The Journal of Chemical Physics, 770 (7):0
3654–3665, 1982.

[Becke, 1993]
Axel D Becke.
A new mixing of Hartree–Fock and local density-functional
The Journal of chemical physics, 980 (2):0
1372–1377, 1993.

[Levine et al., 2006]
Benjamin G Levine, Chaehyuk Ko, Jason Quenneville, and Todd J MartÍnez.
Conical intersections and double excitations in time-dependent
density functional theory.
Molecular Physics, 1040 (5-7):0 1039–1051,

[Filatov, 2013]
Michael Filatov.
Assessment of density functional methods for obtaining geometries at
conical intersections in organic molecules.
Journal of chemical theory and computation, 90
(10):0 4526–4541, 2013.

[Mead and Truhlar, 1982]
C Alden Mead and Donald G Truhlar.
Conditions for the definition of a strictly diabatic electronic
basis for molecular systems.
The Journal of Chemical Physics, 770 (12):0
6090–6098, 1982.

[Baer and Englman, 1992]
Michael Baer and Robert Englman.
A study of the diabatic electronic representation within the
Born-Oppenheimer approximation.
Molecular Physics, 750 (2):0 293–303, 1992.

[Toniolo et al., 2005]
A Toniolo, C Ciminelli, Maurizio Persico, and TJ Martínez.
Simulation of the photodynamics of azobenzene on its first excited
state: Comparison of full multiple spawning and surface hopping treatments.
The Journal of chemical physics, 1230 (23):0
234308, 2005.

[Köppel et al., 2001]
Horst Köppel, Joachim Gronki, and Susanta Mahapatra.
Construction scheme for regularized diabatic states.
The Journal of Chemical Physics, 1150 (6):0
2377, 2001.

[Yue et al., 2018]
Ling Yue, Yajun Liu, and Chaoyuan Zhu.
Performance of TDDFT with and without spin-flip in trajectory
surface hopping dynamics: cis$\rightleftharpoons$trans azobenzene
Physical Chemistry Chemical Physics, 200
(37):0 24123–24139, 2018.

[Bandara et al., 2010]
HM Dhammika Bandara, Tracey R Friss, Miriam M Enriquez, William Isley,
Christopher Incarvito, Harry A Frank, Jose Gascon, and Shawn C Burdette.
Proof for the concerted inversion mechanism in the trans$\to$cis
isomerization of azobenzene using hydrogen bonding to induce isomer locking.
The Journal of organic chemistry, 750 (14):0
4817–4827, 2010.

[Bandara et al., 2011]
HM Dhammika Bandara, Shannon Cawley, José A Gascón, and Shawn C
Short-circuiting azobenzene photoisomerization with
electron-donating substituents and reactivating the photochemistry with
chemical modification, 2011.

[Gómez-Bombarelli et al., 2016]
Rafael Gómez-Bombarelli, Jorge Aguilera-Iparraguirre, Timothy D Hirzel,
David Duvenaud, Dougal Maclaurin, Martin A Blood-Forsythe, Hyun Sik Chae,
Markus Einzinger, Dong-Gwang Ha, Tony Wu, et al.
Design of efficient molecular organic light-emitting diodes by a
high-throughput virtual screening and experimental approach.
Nature materials, 150 (10):0 1120–1127, 2016.

[Velema et al., 2014]
Willem A Velema, Wiktor Szymanski, and Ben L Feringa.
Photopharmacology: beyond proof of principle.
Journal of the American Chemical Society, 1360
(6):0 2178–2191, 2014.

[Bandara and Burdette, 2012]
HM Dhammika Bandara and Shawn C Burdette.
Photoisomerization in different classes of azobenzene.
Chemical Society Reviews, 410 (5):0
1809–1825, 2012.

[Dias et al., 1992]
AR Dias, ME Minas Da Piedade, JA Martinho Simoes, JA Simoni, C Teixeira,
HP Diogo, Yang Meng-Yan, and G Pilcher.
Enthalpies of formation of cis-azobenzene and trans-azobenzene.
The Journal of Chemical Thermodynamics, 240
(4):0 439–447, 1992.

[Yu et al., 2014]
Le Yu, Chao Xu, Yibo Lei, Chaoyuan Zhu, and Zhenyi Wen.
Trajectory-based nonadiabatic molecular dynamics without calculating
nonadiabatic coupling in the avoided crossing case:
Trans$\rightleftharpoons$cis photoisomerization in azobenzene.
Physical Chemistry Chemical Physics, 160
(47):0 25883–25895, 2014.

[Qiao et al., 2020]
Zhuoran Qiao, Feizhi Ding, Matthew Welborn, Peter J. Bygrave, Animashree
Anandkumar, Frederick R. Manby, and Thomas F. Miller III.
Multi-task learning for electronic structure to predict and explore
molecular potential energy surfaces.
arXiv preprint arXiv:2011.02680, 2020b.

[Schwalbe-Koda et al., 2021]
Daniel Schwalbe-Koda, Aik Rui Tan, and Rafael Gómez-Bombarelli.
Differentiable sampling of molecular geometries with
uncertainty-based adversarial attacks.
Nature Communications, 120 (5104), 2021.

[Van Voorhis et al., 2010]
Troy Van Voorhis, Tim Kowalczyk, Benjamin Kaduk, Lee-Ping Wang, Chiao-Lun
Cheng, and Qin Wu.
The diabatic picture of electron transfer, reaction barriers, and
molecular dynamics.
Annual review of physical chemistry, 61:0 149–170,

[Schuurman and Yarkony, 2007]
Michael S Schuurman and David R Yarkony.
On the vibronic coupling approximation: A generally applicable
approach for determining fully quadratic quasidiabatic coupled electronic
state Hamiltonians.
The Journal of chemical physics, 1270 (9):0
094104, 2007.

[Shao et al., 2015]
Yihan Shao, Zhengting Gan, Evgeny Epifanovsky, Andrew TB Gilbert, Michael
Wormit, Joerg Kussmann, Adrian W Lange, Andrew Behn, Jia Deng, Xintian Feng,
et al.
Advances in molecular quantum chemistry contained in the Q-Chem 4
program package.
Molecular Physics, 1130 (2):0 184–215, 2015.

[Nosé, 1984]
Shuichi Nosé.
A unified formulation of the constant temperature molecular dynamics
The Journal of chemical physics, 810 (1):0
511–519, 1984.

[Hoover, 1985]
William G Hoover.
Canonical dynamics: Equilibrium phase-space distributions.
Physical review A, 310 (3):0 1695, 1985.

[Hudock et al., 2007]
Hanneli R Hudock, Benjamin G Levine, Alexis L Thompson, Helmut Satzger, David
Townsend, N Gador, Susanne Ullrich, Albert Stolow, and Todd J Martínez.
Ab initio molecular dynamics and time-resolved photoelectron
spectroscopy of electronically excited uracil and thymine.
The Journal of Physical Chemistry A, 1110
(34):0 8500–8508, 2007.

[Mai et al., 2018]
Sebastian Mai, Philipp Marquetand, and Leticia González.
Nonadiabatic dynamics: The SHARC approach.
Wiley Interdisciplinary Reviews: Computational Molecular
Science, 80 (6):0 e1370, 2018.

[Herbert et al., 2016]
John M Herbert, Xing Zhang, Adrian F Morrison, and Jie Liu.
Beyond time-dependent density functional theory using only single
excitations: Methods for computational studies of excited states in complex
Accounts of chemical research, 490 (5):0
931–941, 2016.

[Martínez et al., 1997]
TJ Martínez, M Ben-Nun, and RD Levine.
Molecular collision dynamics on several electronic states.
The Journal of Physical Chemistry A, 1010
(36):0 6389–6402, 1997.

[Kryachko and Yarkony, 2000]
Eugene S Kryachko and David R Yarkony.
Diabatic bases and molecular properties.
International Journal of Quantum Chemistry, 760
(2):0 235–243, 2000.

[Hoyer et al., 2016]
Chad E Hoyer, Kelsey Parker, Laura Gagliardi, and Donald G Truhlar.
The DQ and DQ$\Phi$ electronic structure diabatization methods:
Validation for general applications.
The Journal of chemical physics, 1440 (19):0
194101, 2016.

[Nakamura and Truhlar, 2001]
Hisao Nakamura and Donald G Truhlar.
The direct calculation of diabatic states based on configurational
The Journal of Chemical Physics, 1150 (22):0
10353–10372, 2001.

[Plasser et al., 2019]
Felix Plasser, Sandra Gómez, Maximilian FSJ Menger, Sebastian Mai, and
Leticia González.
Highly efficient surface hopping dynamics using a linear vibronic
coupling model.
Physical Chemistry Chemical Physics, 210 (1):0
57–69, 2019.

[Guan et al., 2017]
Yafu Guan, Bina Fu, and Dong H Zhang.
Construction of diabatic energy surfaces for LiFH with artificial
neural networks.
The Journal of chemical physics, 1470 (22):0
224307, 2017.

[Shu et al., 2021]
Yinan Shu, Zoltan Varga, Antonio Gustavo Sampaio de Oliveira-Filho, and
Donald G Truhlar.
Permutationally restrained diabatization by machine intelligence.
Journal of Chemical Theory and Computation, 170
(2):0 1106–1116, 2021.

[Richter et al., 2011]
Martin Richter, Philipp Marquetand, Jesús González-Vázquez, Ignacio
Sola, and Leticia González.
SHARC: ab initio molecular dynamics with surface hopping in
the adiabatic representation including arbitrary couplings.
Journal of chemical theory and computation, 70
(5):0 1253–1258, 2011.

[Runge and Gross, 1984]
Erich Runge and Eberhard KU Gross.
Density-functional theory for time-dependent systems.
Physical Review Letters, 520 (12):0 997, 1984.

[Gozem et al., 2014]
Samer Gozem, Federico Melaccio, Alessio Valentini, Michael Filatov, Miquel
Huix-Rotllant, Nicolas Ferré, Luis Manuel Frutos, Celestino Angeli,
Anna I Krylov, Alexander A Granovsky, et al.
Shape of multireference, equation-of-motion coupled-cluster, and
density functional theory potential energy surfaces at a conical
Journal of chemical theory and computation, 100
(8):0 3074–3084, 2014.

[Nikiforov et al., 2014]
Alexander Nikiforov, Jose A Gamez, Walter Thiel, Miquel Huix-Rotllant, and
Michael Filatov.
Assessment of approximate computational methods for conical
intersections and branching plane vectors in organic molecules.
The Journal of chemical physics, 1410 (12):0
124122, 2014.

[Malkin and Fischer, 1962]
Shmuel Malkin and Ernst Fischer.
Temperature dependence of photoisomerization. Part II. Quantum
yields of cis-trans isomerizations in azo-compounds.
The Journal of Physical Chemistry, 660 (12):0
2482–2486, 1962.

[Rau and Yu-Quan, 1988]
Hermann Rau and Shen Yu-Quan.
Photoisomerization of sterically hindered azobenzenes.
Journal of Photochemistry and Photobiology A: Chemistry,
420 (2-3):0 321–327, 1988.

[Mohr and Grummt, 2004]
Gerhard J Mohr and Ulrich-W Grummt.
Photochemistry of the amine-sensor dye 4-N,
Journal of Photochemistry and Photobiology A: Chemistry,
1630 (3):0 341–345, 2004.

[Knie et al., 2014]
Christopher Knie, Manuel Utecht, Fangli Zhao, Hannes Kulla, Sergey Kovalenko,
Albert M Brouwer, Peter Saalfrank, Stefan Hecht, and David Bléger.
ortho‐Fluoroazobenzenes: Visible light switches with very
long‐lived Z isomers.
Chemistry–A European Journal, 200 (50):0
16492–16501, 2014.

[Moreno et al., 2015]
J Moreno, M Gerecke, AL Dobryakov, IN Ioffe, AA Granovsky, D Bléger,
S Hecht, and SA Kovalenko.
Two-photon-induced versus one-photon-induced isomerization dynamics
of a bistable azobenzene derivative in solution.
The Journal of Physical Chemistry B, 1190
(37):0 12281–12288, 2015.

[Rau, 1984]
Hermann Rau.
Further evidence for rotation in the $\pi$, $\pi^*$ and inversion in
the n, $\pi^*$ photoisomerization of azobenzenes.
Journal of photochemistry, 260 (2-3):0
221–225, 1984.

[Fischer, 1967]
Ernst Fischer.
Calculation of photostationary states in systems A
$\rightleftharpoons$ B when only A is known.
The Journal of Physical Chemistry, 710 (11):0
3704–3706, 1967.

[Zimmerman et al., 1958]
George Zimmerman, Lue-Yung Chow, and Un-Jin Paik.
The photochemical isomerization of azobenzene.
Journal of the American Chemical Society, 800
(14):0 3528–3531, 1958.

[Marquetand et al., 2011]
Philipp Marquetand, Martin Richter, Jesús González-Vázquez, Ignacio
Sola, and Leticia González.
Nonadiabatic ab initio molecular dynamics including spin–orbit
coupling and laser fields.
Faraday discussions, 153:0 261–273, 2011.

[Bortolus and Monti, 1979]
Pietro Bortolus and Sandra Monti.
Cis-trans photoisomerization of azobenzene. Solvent and triplet
donors effects.
Journal of Physical Chemistry, 830 (6):0
648–652, 1979.

[Ronayette et al., 1974]
Jacques Ronayette, René Arnaud, Philippe Lebourgeois, and Jacques Lemaire.
Isomérisation photochimique de l'azobenzène en solution. i.
Canadian Journal of Chemistry, 520 (10):0
1848–1857, 1974.

[Bortolus and Monti, 1987]
Pietro Bortolus and Sandra Monti.
Cis $rightleftharpoons$ trans photoisomerization of
azobenzene-cyclodextrin inclusion complexes.
Journal of Physical Chemistry, 910 (19):0
5046–5050, 1987.

[Tiberio et al., 2010]
Giustiniano Tiberio, Luca Muccioli, Roberto Berardi, and Claudio Zannoni.
How does the trans-cis photoisomerization of azobenzene take place in
organic solvents?
ChemPhysChem, 110 (5):0 1018, 2010.

[Granucci and Persico, 2007]
Giovanni Granucci and Maurizio Persico.
Excited state dynamics with the direct trajectory surface hopping
method: azobenzene and its derivatives as a case study.
Theoretical Chemistry Accounts, 1170 (5):0
1131–1143, 2007a.

[Landry and Subotnik, 2012]
Brian R Landry and Joseph E Subotnik.
How to recover Marcus theory with fewest switches surface hopping:
Add just a touch of decoherence.
The Journal of chemical physics, 1370 (22):0
22A513, 2012.

[Chen and Reichman, 2016]
Hsing-Ta Chen and David R Reichman.
On the accuracy of surface hopping dynamics in condensed phase
non-adiabatic problems.
The Journal of Chemical Physics, 1440 (9):0
094104, 2016.

[Zobel et al., 2021]
J Patrick Zobel, Moritz Heindl, Felix Plasser, Sebastian Mai, and Leticia
Surface hopping dynamics on vibronic coupling models.
Accounts of chemical research, 540 (20):0
3760–3771, 2021.

[Ibele et al., 2021]
Lea M Ibele, Yorick Lassmann, Todd J Martínez, and Basile FE Curchod.
Comparing (stochastic-selection) ab initio multiple spawning with
trajectory surface hopping for the photodynamics of cyclopropanone, fulvene,
and dithiane.
The Journal of Chemical Physics, 1540 (10):0
104110, 2021.

[Paszke et al., 2019]
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory
Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban
Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan
Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith
PyTorch: An imperative style, high-performance deep learning
In Advances in Neural Information Processing Systems 32, pages
8024–8035. 2019.

[Klicpera et al., 2020]
Johannes Klicpera, Janek Groß, and Stephan Günnemann.
Directional message passing for molecular graphs.
arXiv preprint arXiv:2003.03123, 2020.

[Westermayr and Marquetand, 2020]
Julia Westermayr and Philipp Marquetand.
Deep learning for UV absorption spectra with SchNarc: First steps
toward transferability in chemical compound space.
The Journal of Chemical Physics, 1530 (15):0
154112, 2020b.

[Barbatti and Sen, 2016]
Mario Barbatti and Kakali Sen.
Effects of different initial condition samplings on photodynamics
and spectrum of pyrrole.
International Journal of Quantum Chemistry, 1160
(10):0 762–771, 2016.

[Larsen et al., 2017]
Ask Hjorth Larsen, Jens Jørgen Mortensen, Jakob Blomqvist, Ivano E Castelli,
Rune Christensen, Marcin Dułak, Jesper Friis, Michael N Groves, Bjørk
Hammer, Cory Hargus, Eric D Hermes, Paul C Jennings, Peter Bjerre Jensen,
James Kermode, John R Kitchin, Esben Leonhard Kolsbjerg, Joseph Kubal,
Kristen Kaasbjerg, Steen Lysgaard, Jón Bergmann Maronsson, Tristan Maxson,
Thomas Olsen, Lars Pastewka, Andrew Peterson, Carsten Rostgaard, Jakob
Schiøtz, Ole Schütt, Mikkel Strange, Kristian S Thygesen, Tejs Vegge, Lasse
Vilhelmsen, Michael Walter, Zhenhua Zeng, and Karsten W Jacobsen.
The atomic simulation environment—a Python library for working
with atoms.
Journal of Physics: Condensed Matter, 290
(27):0 273002, 2017.

[Verlet, 1967]
Loup Verlet.
Computer “experiments” on classical fluids. I. Thermodynamical
properties of Lennard-Jones molecules.
Physical review, 1590 (1):0 98, 1967.

[Yue et al., 2017]
Ling Yue, Le Yu, Chao Xu, Yibo Lei, Yajun Liu, and Chaoyuan Zhu.
Benchmark performance of global switching versus local switching for
trajectory surface hopping molecular dynamics simulation:
Cis$\rightleftharpoons$trans azobenzene photoisomerization.
ChemPhysChem, 180 (10):0 1274–1287, 2017.

[Mai et al., 2015]
Sebastian Mai, Philipp Marquetand, and Leticia González.
A general method to describe intersystem crossing dynamics in
trajectory surface hopping.
International Journal of Quantum Chemistry, 1150
(18):0 1215–1231, 2015.

[Plasser et al., 2012]
Felix Plasser, Giovanni Granucci, Jiri Pittner, Mario Barbatti, Maurizio
Persico, and Hans Lischka.
Surface hopping dynamics using a locally diabatic formalism: Charge
transfer in the ethylene dimer cation and excited state dynamics in the
2-pyridone dimer.
The Journal of chemical physics, 1370 (22):0
22A514, 2012.

[128]
Neural Force Field.
< https://github.com/learningmatter-mit/NeuralForceField>.

[Granucci and Persico, 2007]
Giovanni Granucci and Maurizio Persico.
Critical appraisal of the fewest switches algorithm for surface
The Journal of chemical physics, 1260 (13):0
134114, 2007b.

[Becke, 1988]
Axel D Becke.
Density-functional exchange-energy approximation with correct
asymptotic behavior.
Physical review A, 380 (6):0 3098, 1988.

[Perdew, 1986]
John P Perdew.
Density-functional approximation for the correlation energy of the
inhomogeneous electron gas.
Physical Review B, 330 (12):0 8822, 1986.

[Van Voorhis and Head-Gordon, 2002]
Troy Van Voorhis and Martin Head-Gordon.
A geometric approach to direct minimization.
Molecular Physics, 1000 (11):0 1713–1721,

[Smith et al., 2019]
Justin S Smith, Benjamin T Nebgen, Roman Zubatyuk, Nicholas Lubbers, Christian
Devereux, Kipton Barros, Sergei Tretiak, Olexandr Isayev, and Adrian
Outsmarting quantum chemistry through transfer learning.

[Bahdanau et al., 2014]
Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio.
Neural machine translation by jointly learning to align and
arXiv preprint arXiv:1409.0473, 2014.

[Kim et al., 2017]
Yoon Kim, Carl Denton, Luong Hoang, and Alexander M Rush.
Structured attention networks.
arXiv preprint arXiv:1702.00887, 2017.

[Vaswani et al., 2017]
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones,
Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.
Attention is all you need.
In Advances in neural information processing systems, pages
5998–6008, 2017.

[Veličković et al., 2017]
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero,
Pietro Lio, and Yoshua Bengio.
Graph attention networks.
arXiv preprint arXiv:1710.10903, 2017.

[Kurita et al., 2014]
Motoki Kurita, Miku Makihara, and Hideyuki Nakano.
Photochromic reactions of 4-[bis (9, 9-dimethylfluoren-2-yl) amino]
azobenzene in low molecular-mass organogels.
Soft materials, 120 (1):0 42–46, 2014.

[Gegiou et al., 1968]
Dina Gegiou, KA Muszkat, and Ernst Fischer.
Temperature dependence of photoisomerization. V. Effect of
substituents on the photoisomerization of stilbenes and azobenzenes.
Journal of the American Chemical Society, 900
(15):0 3907–3918, 1968.

[Sierocki et al., 2006]
Paul Sierocki, Huub Maas, Patrick Dragut, Gabriele Richardt, Fritz Vögtle,
Luisa De Cola, Fred Brouwer, and Jeffrey I Zink.
Photoisomerization of azobenzene derivatives in nanostructured
The Journal of Physical Chemistry B, 1100
(48):0 24390–24398, 2006.

[Siewertsen et al., 2009]
Ron Siewertsen, Hendrikje Neumann, Bengt Buchheim-Stehn, Rainer Herges,
Christian Nather, Falk Renth, and Friedrich Temps.
Highly efficient reversible Z-E photoisomerization of a bridged
azobenzene with visible light through resolved $S_1$(n$\pi^*$) absorption
Journal of the American Chemical Society, 1310
(43):0 15594–15595, 2009.

[Lentes et al., 2019]
Pascal Lentes, Eduard Stadler, Fynn Röhricht, Arne Brahms, Jens
Gröbner, Frank D Sönnichsen, Georg Gescheidt, and Rainer Herges.
Nitrogen bridged diazocines: Photochromes switching within the
near-infrared region with high quantum yields in organic solvents and in
Journal of the American Chemical Society, 1410
(34):0 13592–13600, 2019.

[Birnbaum and Style, 1954]
PP Birnbaum and DWG Style.
The photo-isomerization of some azobenzene derivatives.
Transactions of the Faraday Society, 50:0 1192–1196,

[Nagamani et al., 2005]
S Anitha Nagamani, Yasuo Norikane, and Nobuyuki Tamaoki.
Photoinduced hinge-like molecular motion: studies on xanthene-based
cyclic azobenzene dimers.
The Journal of organic chemistry, 700 (23):0
9304–9313, 2005.

[Olmsted III et al., 1983]
John Olmsted III, Jerry Lawrence, and Geary G Yee.
Photochemical storage potential of azobenzenes.
Solar Energy, 300 (3):0 271–274, 1983.

[Chen et al., 2012]
Tao Chen, Atsushi Yamaguchi, Kazumasa Igarashi, Naoya Nakagawa, Hidenori
Nishioka, Hiroyuki Asanuma, and Mikio Yamashita.
Ultrafast photoisomerization and its single-shot pump pulse
efficiency of trans-azobenzene derivative: Compound for photosensitive DNA.
Optics Communications, 2850 (6):0 1206–1211,

[Goulet-Hanssens et al., 2014]
Alexis Goulet-Hanssens, T Christopher Corkery, Arri Priimagi, and Christopher J
Effect of head group size on the photoswitching applications of
azobenzene Disperse Red 1 analogues.
Journal of Materials Chemistry C, 20 (36):0
7505–7512, 2014.

[Mourot et al., 2011]
Alexandre Mourot, Michael A Kienzler, Matthew R Banghart, Timm Fehrentz,
Florian ME Huber, Marco Stein, Richard H Kramer, and Dirk Trauner.
Tuning photochromic ion channel blockers.
ACS chemical neuroscience, 20 (9):0 536–543,

[Sadovski et al., 2009]
Oleg Sadovski, Andrew A Beharry, Fuzhong Zhang, and G Andrew Woolley.
Spectral tuning of azobenzene photoswitches for biological
Angewandte Chemie International Edition, 480
(8):0 1484–1486, 2009.

[Samanta et al., 2012]
Subhas Samanta, Chuanguang Qin, Alan J Lough, and G Andrew Woolley.
Bidirectional photocontrol of peptide conformation with a bridged
azobenzene derivative.
Angewandte Chemie International Edition, 510
(26):0 6452–6455, 2012.

[Samanta et al., 2013]
Subhas Samanta, Andrew A Beharry, Oleg Sadovski, Theresa M McCormick,
Amirhossein Babalhavaeji, Vince Tropepe, and G Andrew Woolley.
Photoswitching azo compounds in vivo with red light.
Journal of the American Chemical Society, 1350
(26):0 9777–9784, 2013a.

[Gascón-Moya et al., 2015]
Marta Gascón-Moya, Arnau Pejoan, Mercè Izquierdo-Serra, Silvia Pittolo,
Gisela Cabré, Jordi Hernando, Ramon Alibés, Pau Gorostiza, and
Félix Busqué.
An optimized glutamate receptor photoswitch with sensitized
azobenzene isomerization.
The Journal of organic chemistry, 800 (20):0
9915–9925, 2015.

[Rau et al., 1990]
Hermann Rau, Gerhard Greiner, Guenter Gauglitz, and Herbert Meier.
Photochemical quantum yields in the A (+$h\nu$)
$\rightleftharpoons$ B (+$h\nu$, $\Delta$) system when only the spectrum
of A is known.
Journal of physical chemistry, 940 (17):0
6523–6524, 1990.

[Albini et al., 1983]
Angelo Albini, Elisa Fasani, and Silvio Pietra.
The photochemistry of azo dyes. Photoisomerisation versus
photoreduction from 4-diethylaminoazobenzene and
Journal of the Chemical Society, Perkin Transactions 2,
0 (7):0 1021–1024, 1983.

[Lee et al., 2009]
Yu Jin Lee, Sung Ik Yang, Dae Seung Kang, and Sang-Woo Joo.
Solvent dependent photo-isomerization of 4-dimethylaminoazobenzene
carboxylic acid.
Chemical Physics, 3610 (3):0 176–179, 2009.

[Koller et al., 2007]
FO Koller, C Sobotta, TE Schrader, T Cordes, WJ Schreier, A Sieg, and P Gilch.
Slower processes of the ultrafast photo-isomerization of an
azobenzene observed by IR spectroscopy.
Chemical Physics, 3410 (1-3):0 258–266, 2007.

[Umemoto et al., 2014]
Tetsuro Umemoto, Yuta Ohtani, Takamasa Tsukamoto, Tetsuya Shimada, and Shinsuke
Pinning effect for photoisomerization of a dicationic azobenzene
derivative by anionic sites of the clay surface.
Chemical Communications, 500 (3):0 314–316,

[Tahara et al., 1997]
Ruriko Tahara, Tatsuya Morozumi, Hiroshi Nakamura, and Masatsugu Shimomura.
Photoisomerization of azobenzocrown ethers. Effect of complexation
of alkaline earth metal ions.
The Journal of Physical Chemistry B, 1010
(39):0 7736–7743, 1997.

[Dong et al., 2015]
M Dong, A Babalhavaeji, MJ Hansen, L Kalman, and GA Woolley.
Red, far-red, and near infrared photoswitches based on azonium
Chemical Communications, 510 (65):0
12981–12984, 2015.

[Samanta et al., 2013]
Subhas Samanta, Amirhossein Babalhavaeji, Ming-xin Dong, and G Andrew Woolley.
Photoswitching of ortho-substituted azonium ions by red light in
whole blood.
Angewandte Chemie, 1250 (52):0 14377–14380,

[Carroll et al., 2015]
Elizabeth C Carroll, Shai Berlin, Joshua Levitz, Michael A Kienzler, Zhe Yuan,
Dorte Madsen, Delmar S Larsen, and Ehud Y Isacoff.
Two-photon brightness of azobenzene photoswitches designed for
glutamate receptor optogenetics.
Proceedings of the National Academy of Sciences, 1120
(7):0 E776–E785, 2015.

[Archut et al., 1998]
Andreas Archut, Fritz Vögtle, Luisa De Cola, Gianluca Camillo Azzellini,
Vincenzo Balzani, PS Ramanujam, and Rolf H Berg.
Azobenzene-functionalized cascade molecules: photoswitchable
supramolecular systems.
Chemistry–A European Journal, 40 (4):0
699–706, 1998.

[Lednev et al., 1998]
Igor K Lednev, Tian-Qing Ye, Laurence C Abbott, Ronald E Hester, and John N
Photoisomerization of a capped azobenzene in solution probed by
ultrafast time-resolved electronic absorption spectroscopy.
The Journal of Physical Chemistry A, 1020
(46):0 9161–9166, 1998.

[Beharry et al., 2008]
Andrew A Beharry, Oleg Sadovski, and G Andrew Woolley.
Photo-control of peptide conformation on a timescale of seconds with
a conformationally constrained, blue-absorbing, photo-switchable linker.
Organic & biomolecular chemistry, 60 (23):0
4323–4332, 2008.

Supplementary Information

§ EXTENDED METHODS

§.§ Relevance of adiabatic gap

Accurate prediction of the energy gap is crucial for generating reliable NAMD results. The gap controls the hopping rate, and in turn observables like the photoisomerization quantum yield, excited state lifetime and time-resolved photoelectron spectrum [Hudock et al., 2007]. To understand the importance of the gap, consider that in most approaches, the hopping rate is determined by the derivative form of the NACV [Tully, 1990, Mai et al., 2018]. The derivative coupling between states $n$ and $m$ can be written as $\vec{g}_{nm} / \Delta E_{nm}$, where $\Delta E_{nm}$ is the energy gap and $\vec{g}_{nm}$ is the force coupling [Schuurman and Yarkony, 2007, Herbert et al., 2016]. The gap in the denominator leads to singular derivative coupling at CIs, and therefore to a guaranteed hop. The coupling can alternatively be obtained from the gap alone from through its first and second derivatives [Westermayr et al., 2020]. The energy difference also features prominently in the Zhu-Nakamura method, in which the hopping rate is approximately exponential in the square of the gap [Yu et al., 2014]. Therefore, it is crucial to accurately predict the energy gap in any NAMD simulation.

§.§ Diabatization

The adiabatic energies of a system, $\{ E_n(\vec{R}) \}$, are the energies produced by a QC calculation. They depend on the nuclear coordinates $\vec{R}$, and form the usual Born-Oppenheimer PESs. In the adiabatic basis the nuclear kinetic energy operator, related to $\nabla_{R}^2$, is non-diagonal. Its off-diagonal elements are related to the NACVs, which generate transitions between adiabatic PESs [Martínez et al., 1997]. The derivative form of the NACV also becomes singular at CIs, which is an undesirable numerical property. The diabatic basis is designed to remove this singularity [Kryachko and Yarkony, 2000]. The diabatic Hamiltonian $\mathbf{H}_d(\vec{R})$ is a rotation of the adiabatic energies into a new basis, given by $\mathbf{H}_d(\vec{R}) = \mathbf{U}(\vec{R}) \hspace{0.05cm} \mathrm{diag}(\{ E_n(\vec{R}) \} ) \hspace{0.05cm} \mathbf{U}^{\dagger}(\vec{R})$. Here $\mathrm{diag}$ denotes a diagonal matrix, and $\mathbf{U}(\vec{R})$ is a rotation matrix that depends on the nuclear coordinates (see Eq. (<ref>eq:u_def)). $\mathbf{H}_d$ can also be viewed as the clamped nucleus Hamiltonian expressed in the basis of diabatic wave functions. These wave functions are given by $\psi_{\mathrm{d}, n}(\vec{r}; \vec{R})= \sum_m U_{nm} (\vec{R}) \hspace*{0.05cm} \psi_{\mathrm{ad}, m}(\vec{r}; \vec{R})$, where $\psi_{\mathrm{ad}, m}$ is the $m^{\mathrm{th}}$ adiabatic wave function and $\vec{r}$ denotes the electronic coordinates.

The diabatic states are defined such that the nuclear kinetic energy is approximately diagonal [Kryachko and Yarkony, 2000]. The states are instead coupled through off-diagonal elements in the potential energy matrix, known as diabatic couplings. The diabatic couplings are smooth functions of the nuclear coordinates. The diagonal elements are also smooth, maintaining their orbital character and switching energy ordering through a CI (Fig. <ref>fig:diabat(b)).
In many applications, diabatic states are preferred over adiabatic ones because the singular coupling is removed. Here we prefer them because diabatic energies are easier to fit than adiabatic energies. That is, even if one is only interested in adiabatic energies and not NACVs, it is easier to learn the diabatic energies and then diagonalize $\mathbf{H}_d$ than to learn the adiabatic energies directly. This is because the diabatic states possess no CI cusps or avoided crossings (Fig. <ref>fig:diabat), and are thus more easily fit by interpolating functions such as neural networks. As discussed below, diabatic states improve the model accuracy for species outside the training set.

While many diabatization methods exist, the most common ones cannot be straightforwardly applied to the current problem. Property-based methods were developed for charge-transfer type problems [Hoyer et al., 2016], while orbital-based approaches [Nakamura and Truhlar, 2001] are not implemented for TDDFT.
Diabatic models that are parameterized by electronic structure data [Plasser et al., 2019, Schuurman and Yarkony, 2007, Williams and Eisfeld, 2018] are not designed for systems undergoing large distortions. Approaches that solve for the adiabatic-to-diabatic transformation matrix [Guan et al., 2017] are difficult to implement in practice, because the matrix varies rapidly near a CI.

Recent work introduced neural network diabatization based on reference geometries [Shu and Truhlar, 2020, Shu et al., 2021]. In this procedure one assumes that the diabatic Hamiltonian is diagonal at a set of known reference geometries. At such geometries the elements of $\mathbf{H}_d$ are equal to the adiabatic energies, but possibly reordered. For example, for two states and two reference geometries, one would have $\mathbf{H}_d = \mathrm{diag}(E_0, E_1)$ at the first geometry and $\mathbf{H}_d = \mathrm{diag}(E_1, E_0)$ at the second (Fig. <ref>fig:diabat(b)). A neural network is then trained to produce $\mathbf{H}_d$, such that its eigenvalues are always equal to the adiabatic energies, and its elements are as above at the reference geometries. These two constraints generate $\mathbf{H}_d$ at all intermediate geometries.

This method was successfully applied to thiophenol dissociation, yielding results consistent with the fourfold way [Shu and Truhlar, 2020]. However, as shown by the NACV error in Table <ref>, this method alone does not generate true diabatic states for azobenzene. To understand why, consider that near a CI, the true diabatic coupling $d_{01}$ must closely resemble the coupling shown in Fig. <ref>fig:diabat(b) ($d_{nm}$ is shorthand for $(\mathbf{H}_d)_{nm}$). This is because $d_{01}$ must be linear in displacements about a CI, and quadratic only for Renner–Teller type intersections [Köppel et al., 2001]. However, for two diabatic states, only the square of the diabatic coupling $\vert d_{01} \vert^2$ enters into the expression for the adiabatic energies. The model might then generate an off-diagonal element similar to $|d_{01}|$. This error would incur no penalty on the network, because the diagonal elements would properly switch ordering and the adiabatic energies would be correct. Hence the model would not generate smooth diabatic states.

To remedy this issue we used the combined loss described in the Methods section. The NACV component of this loss was also used in Ref. [Guan et al., 2019]. We note that the number of diabatic states and choice of orderings in Eq. (<ref>eq:ref_geom_azo) may depend on the system. For azobenzene we were interested in fitting $M_{\mathrm{ad}}=2$ adiabatic states, and in general one should use $M_{\mathrm{d}} \geq M_{\mathrm{ad}}$ diabatic states. In this work chose $M_{\mathrm{d}}=3$ because it significantly improved on the results of $M_{\mathrm{d}} = 2$. In fact, with $M_{\mathrm{d}}=2$, the $R^2$ correlation of the force NACV was negative. The orderings in Eq. (<ref>eq:ref_geom_azo) were chosen by taking a small sample of azobenzene configurations, training several small models with different diabatic orderings, and picking the one with the best results. Thus no system knowledge is required to choose the orderings. The only system knowledge required is the set of reference geometries[Typically the reference geometries are reactants and products. For reactions in which the product is not known a priori, one might use the following approach. First, train an adiabatic network on a single molecule. We have found that adiabatic models match or outperform diabatic ones for single species. Then use the model to simulate dynamics. By using active learning to improve the model's coverage of configurational space, the simulations can eventually discover the product. If derivatives of the molecule are expected to have similar reactions, then their products are also now known. With knowledge of the reactants and products, and hence the reference geometries, one can now build a diabatic model. If this is not possible, then one can train the model without $\mathcal{L}_{\mathrm{ref}}$. Table <ref> shows that this will likely decrease the model performance near CIs by a small amount.], but the reference loss can be omitted with negligible impact on model performance.

The diabatic model leads to true CIs, where the ground and excited state energies are exactly equal. To see why, consider Fig. <ref>fig:results(b), which shows the two diagonal diabatic elements (red and blue) and the off-diagonal coupling (light gray). The off-diagonal element passes linearly through zero in the $\vec{h}$ direction. The diagonal elements cross in the $\vec{g}$ direction, and so $\Delta \equiv d_{11} - d_{00}$ passes linearly through zero. Therefore, one can begin at a geometry for which $\Delta = 0$, and move in the $\vec{h}$ direction until $d_{01}=0$ without changing $\Delta$. The final geometry will therefore be a CI. The fact that both $d_{01}$ and $\Delta$ are locally linear, and hence must pass through zero, is known theoretically [Köppel et al., 2001] and properly predicted by the model.

§.§ Network loss

The neural network loss terms are defined as:

ℒ_core =        ∑_n  ρ_E_n  ·mse (E_n ) + ρ_f_n ·mse (f⃗_n)

+ ∑_n > m ρ_ΔE_nm  ·mse( ΔE_nm )

ℒ_ref =         ∑_n ρ_ref ·mse ( d_nn(R⃗) ) |_R⃗ ∈{cis, trans}

ℒ_nacv  =         ∑_n > m  ρ_nacv ·mse( A_n( R⃗ ) A_m ( R⃗ )  g⃗_nm(R⃗) ),

where $\Delta E_{nm} = E_n - E_m$ is the energy gap. Each loss function is a sum over mean squared error (MSE) terms scaled by different weights $\rho$. For scalar molecular quantities $X$ the loss is given by $\mathrm{MSE}(X) = \sum_{j=1}^{M} (X_j - \hat{X}_j)^2 / M$, where $\hat{X}$ is the predicted quantity and the sum is over $M$ geometries in a batch. For atomic vectorial quantities the mean is additionally over the $3 L_j$ vector components, where $L_j$ is the number of atoms in the $j^{\mathrm{th}}$ geometry.

The loss term $\mathcal{L}_{\mathrm{core}}$ contains the usual energy and force losses, plus an additional term for gap errors. While the gap term was not used in previous work, we found it to be crucial for the systems studied here. Indeed, without this loss term, one would expect the gap MAE in adiabatic models to be the geometric sum of the energy MAEs. This is what we found when using $\mathcal{L}_{\mathrm{core}}$ without the gap loss. Table <ref> shows that, when using the full core loss, the gap MAE is actually lower than each individual energy MAE.

The reference loss, $\mathcal{L}_{\mathrm{ref}}$, is a sum over geometries $\vec{R}$ which are considered to be cis or trans. At these geometries the target $d_{nn}$ are given by Eq. (<ref>eq:ref_geom_azo). A geometry is considered to be cis or trans if its central CNNC atoms deviate from those of the equilibrium structure by <0.15 Å. The distance is computed as the root-mean-square deviation (RMSD) among the atoms after alignment.

The NACV loss imposes diabaticity. It involves the force NACV $\vec{g}_{nm}$, and a phase correction $A_{n}(\vec{R}) = \pm 1$. The phase correction is chosen separately for each geometry to minimize the error between the predicted and target force coupling. This factor is necessary because each adiabatic wavefunction can have an arbitrary sign [Richter et al., 2011]. The signs cancel for diagonal terms like $E_n$ and $\vec{f}_n$, but not for off-diagonal terms like $\vec{g}_{nm}$. The $A_{n}$ account for these arbitrary sign changes.

§.§ Data generation

To train a useful model one must generate reliable QC data. TDDFT [Runge and Gross, 1984] typically offers a good compromise between speed and accuracy for modeling excited states. However, it has known instabilities near CIs [Levine et al., 2006], and, as a result of the Brillouin theorem, produces the wrong branching space dimensionality for $S_0/S_1$ intersections [Gozem et al., 2014]. These issues, which can be traced back to the single-reference description of the excitation, can be partially alleviated with SF-TDDFT [Shao et al., 2003]. In this approach the excitation is performed with respect to a high-spin reference state. This yields some transitions that have double-excitation character with respect to the singlet ground state. Here we used SF-TDDFT with the BHHLYP functional [Becke, 1993] and the 6-31G* basis [Francl et al., 1982]. SF-TDDFT/BHHLYP is well-benchmarked for CIs in a number of molecules [Filatov, 2013, Nikiforov et al., 2014]. Because SF is not spin complete, the singlet states must be identified based on their square spins $\langle S^2 \rangle$. We used the approach of Ref. [Yue et al., 2018], which identifies singlets as the two states with the lowest $\langle S^2 \rangle$ from the three excitations of lowest energy. Calculations were performed with the Q-Chem package [Shao et al., 2015].

Near-CI regions of the combinatorial species were sampled by setting the central CNNC dihedral angle of relaxed geometries to $\pm$90 degrees and/or the CNN/NNC angles to $180$ degrees. The other internal coordinates were not changed. The former corresponds to the rotation pathway and the latter to inversion. This led to 11 possible combinations of rotation and inversion, including pure rotation, pure inversion, concerted inversion, and inversion-assisted rotation [Bandara and Burdette, 2012].

§ MODEL ACCURACY AND ABLATION STUDIES

$E_0$     $E_1$     $\Delta E_{01}$     $(\Delta E)_{\mathrm{small}}$     sgn$(\Delta E)_{\mathrm{small}} $     $\vec{F}_0$     $\vec{F}_1$     $\vec{g}_{01}$

DANN     3.06     3.77     1.89     0.97     -0.29     1.72      2.31     1.36

$-\mathcal{L}_{\mathrm{nacv}}$     2.31     2.49     1.65     1.22     0.52     1.72     2.39     2.21

$-\mathcal{L}_{\mathrm{ref}}$     2.21     2.58     1.68     1.06     0.16     1.72     2.32     1.37

adiabat      2.11     3.24     2.43     1.24     0.60     1.69     2.31     —

adiabat (st. 1)      3.28     2.99     1.88     1.49     0.95     1.67     2.30     —

median       16.86     19.04     22.80     36.04     -36.04     17.27     17.18     2.05

Performance of different ablated models. Each column, apart from $\mathrm{sgn} (\Delta E)_{\mathrm{small}}$, shows the MAE of a different quantity. $\mathrm{sgn} (\Delta E)_{\mathrm{small}}$ is the mean signed error of $(\Delta E)_{\mathrm{small}}$, given by mean{$( E)_{\mathrm{small}}^{\mathrm{pred}} - (\Delta E)_{\mathrm{small}}^{\mathrm{target}}$}. “st. 1” indicates the first stage of training of the adiabatic model, as described in Supplementary Sec. <ref>. The best score in each category is shown in bold.

$E_0$     $E_1$     $\Delta E_{01}$     $(\Delta E_{01})_{\mathrm{small}}$     $\vec{F}_0$     $\vec{F}_1$     $\vec{g}_{01}$

MAE ($\downarrow$)      0.75     0.96     0.72     0.40     0.94     1.13     0.87

$R^2$ ($\uparrow$)      1.00     0.99     0.99     0.97     1.00     0.99     0.88

Test set MAE for the DANN model trained on all species.

Here we compare the DANN model to a model trained without a NACV loss (“$-\mathcal{L}_{\mathrm{nacv}}$”), a model trained without a reference geometry loss (“$-\mathcal{L}_{\mathrm{ref}}$”), an adiabatic model, and a median predictor. The latter predicts the median ground- and excited-state energy for each species, and the median value among all species for other quantities. The MAEs for unseen species are compared in Table <ref>. Half the geometries were sampled randomly and half by proximity to a CI, as described in Supplementary Sec. <ref>. Of particular interest are $(\Delta E)_{\mathrm{small}}$, the MAE of the gap when it is under 0.2 eV, and $\mathrm{sgn}(\Delta E)_{\mathrm{small}}$, the average overestimation of small gaps by the model. DANN actually underestimates the small gaps, indicating that it does not “miss” conical intersections. Using an adiabatic model or removing the NACV loss increases both the error and overestimation of $(\Delta E)_{\mathrm{small}}$. The same is true of removing $\mathcal{L}_{\mathrm{ref}}$, but the effect is smaller. While the changes in $(\Delta E)_{\mathrm{small}}$ appear minor, we show in Supplementary Sec. <ref> that they lead to noticeable quantum yield differences for a number of species.

Also of note is the difference in NACV error among the different models. Removing $\mathcal{L}_{\mathrm{nacv}}$ substantially increases the NACV error, even when $\mathcal{L}_{\mathrm{ref}}$ is used. Since $\mathcal{L}_{\mathrm{nacv}}$ imposes diabaticity, this shows that the reference loss alone does not give accurate diabatic states. We also see that the adiabatic model is much worse at predicting $\Delta E$ away from the CI than the diabatic model. This is because of the large loss weight used for $(\Delta E)_{\mathrm{small}}$ in the second stage of adiabatic training (see Supplementary Sec. <ref>). Indeed, the results after the first stage are much better for $\Delta E$, but far worse for $(\Delta E)_{\mathrm{small}}$. Finally, the models without reference or NACV losses have far lower energy errors than DANN. It may be possible to decrease DANN's energy errors by increasing the weight of the energy loss.

A key benefit of the diabatic model is that it accurately predicts the gap near CIs. It also produces the NACVs, adiabatic energies, and forces all with one model. Without the diabatic model, one might learn the force coupling as a separate property. This has been done in other work by taking the gradient of a learnable scalar [Westermayr et al., 2020, Li et al., 2021]. However, one would still have to divide by the adiabatic gap to obtain the derivative coupling. Thus gap prediction errors from an adiabatic model would still be problematic. A separate option would be to compute the NACV through the Hessian of the gap [Westermayr et al., 2020]. This is quite slow, and would have similar problems with the gap.

Table <ref> shows the test set statistics for the DANN model trained on all species. This model was used for virtual screening, as described in the main text. Here the test set simply consists of 5,000 held-out geometries. Unlike in Table <ref>, the test set contains mostly seen species, and the geometries were not generated from DANN-NAMD with the trained model. The results are quite similar to the “seen species” rows in Table <ref>tab:accuracy.

§ EXPERIMENTAL DATA

We performed an extensive literature search to find quantum yields of azobenzene derivatives. To best compare to the vacuum conditions simulated here, we chose quantum yields measured in non-polar solvents. When multiple results were available in different non-polar solvents, we computed the prediction error relative to the mean experimental value. To compare to the simulated $S_1$ dynamics, we selected yields measured at the peak of the highest-wavelength absorption band, which is a dipole-forbidden $n-\pi^*$ transition in unsubstituted azobenzene. All sources used here specifically reported yields at wavelengths close to the $n-\pi^*$ and $\pi-\pi^*$ absorption peaks, so it was straightforward to choose the appropriate wavelength. Lastly, we only selected yields at or around room temperature, since yields can have a strong temperature dependence [Malkin and Fischer, 1962].

§ SOURCES OF ERROR

§.§ Experimental error

There are two main sources of uncertainty in the experimental quantum yield results. The first is the error in the absorption coefficients $\varepsilon$ of the two isomers. The absorption strengths are used in standard rate equations to compute the yield; see, for example, Refs. [Rau and Yu-Quan, 1988, Mohr and Grummt, 2004] and the supplementary of Refs. [Knie et al., 2014, Moreno et al., 2015]. The error can be traced back to the overlapping absorption spectra of the two isomers, which must be disentangled. The magnitude of the error depends on the method used to compute $\varepsilon$, which itself is related to the year of publication. For example, Ref. [Rau and Yu-Quan, 1988], published in 1988, indicated that the absorption coefficients were the primary source of error, but did not give an estimate for their uncertainty. They isolated the different isomers using chromatography and computed their absorption coefficients separately. When the isomers could not be isolated in sufficient amounts, the authors used the time-dependent absorption approach of Ref. [Rau, 1984]. This method had relative errors of of $\pm$25% for $n-\pi^*$ quantum yields [Rau, 1984] (i.e, yield $\to \mathrm{yield} \times (1 \pm 0.25)$). Ref. [Mohr and Grummt, 2004], published in 2004, reported a yield of $0.7$ to $1.0$ for compound 20, a large range stemming from the overlap of the isomers' spectra and the method used to separate them  [Fischer, 1967]. More modern methods have lower errors; for example, Ref. [Moreno et al., 2015], published in 2015, reported a relative error of only $\pm$10% in the supplementary.

The second important source of error is the overlap of the $n-\pi^*$ and $\pi-\pi^*$ absorption bands. The more strongly the two bands overlap, the higher the error of using only the $S_1$ state in the dynamics.
Moreover, many experiments do not irradiate at the precise maximum of the $n-\pi^*$ band. Often a representative $n-\pi^*$ or $\pi-\pi^*$ wavelength is chosen, and then used for each derivative [Knie et al., 2014, Malkin and Fischer, 1962], even though the derivatives have different absorption peaks.

To get an idea of the range of errors this could introduce, consider the results of Ref. [Zimmerman et al., 1958] from 1958. This work computed the quantum yield of unsubstituted azobenzene over several wavelengths. The trans to cis yield was measured as 21% at 405 nm and 27% at 436 nm ($10^{-3}$ M solution). The difference is because 405 nm light excites more of the $\pi-\pi^*$ band (313 nm) than 436 nm light. 405 nm excitation leads to a lower yield because the $\pi-\pi^*$ yield is only 11%. The difference is even more pronounced for the cis quantum yield, which drops from 55% to 40% moving from 436 nm to 546 nm. This result does not have a simple explanation, as there is no lower energy band beyond 436 nm. In principle, the error due to overlapping bands could be mitigated with an explicit dipole-electric field coupling term [Marquetand et al., 2011], which would excite multiple states in different proportions.
For derivatives with completely overlapping $n-\pi*$ and $\pi-\pi*$ bands, the $S_1$ approximation would incur a maximum error of the $S_2$ yield minus the $S_1$ yield. This is because the $S_2$ state is much brighter than $S_1$, and would therefore dominate the experimental yield. The yield difference is about 10 percentage points for both cis and trans unsubstituted azobenzene [Bortolus and Monti, 1979, Ronayette et al., 1974]. Therefore, in the worst case scenario, we would expect an error of about 10 percentage points from the $S_1$ approximation. The error would likely be closer to the 6% reported in Ref. [Zimmerman et al., 1958] at wavelengths of moderate overlap.

We can further quantify these errors by computing the range of results from different studies.  Measurements of the $t \to c$ azobenzene yield range from 20% to 28% between 1979 and 1987 [Bortolus and Monti, 1979, Rau, 1984, Bortolus and Monti, 1987], though this could also be related to solvent effects (Sec. <ref>).The yields of compounds 9 and 10 range from 16% to 24% and 44% to 50%, respectively, with references from 1962 and 1988 [Malkin and Fischer, 1962, Rau and Yu-Quan, 1988]. Aggregating the results of all three compounds gives an average yield range of 7.3 percentage points, and an average relative error of $\pm 14\%$.

A final source of uncertainty is specific to Refs. [Bandara et al., 2010, Bandara et al., 2011]. Several species were reported to have zero yield over a large range of irradiation wavelengths. However, as noted in the footnote to Table <ref>, the dipole-allowed $S_2$ transition was highly redshifted and thus overlapped strongly with the $S_1$ transition. The absorbance changes after irradiation were quite small, indicating a near-zero yield, but could have been due to the $S_1$ transition. Hence it is possible that the $S_1$ yield is not precisely zero.

§.§ Computational error

There are several sources of computational error. The first is error in the PES, which can be decomposed into error from SF-TDDFT and error from the model. While it is intractable to perform SF-TDDFT for all species with experimental data, our results for unsubstituted azobenzene suggest that it is rather accurate. As reported in the main text, we computed yields of 60 ± 4% and 26 ± 3% for $c \to t$ and $t \to c$, respectively. Experimental measurements between 1979 to 1987 gave $t\to c$ yields between 20% and 28% [Bortolus and Monti, 1979, Rau, 1984, Bortolus and Monti, 1987, Ronayette et al., 1974]. Experimental measurements in 1974 and 1979 gave $c\to t$ yields of 55% and 56% [Bortolus and Monti, 1979, Ronayette et al., 1974], while a measurement in 1958 gave 48% $\pm$ 5% [Zimmerman et al., 1958]. All measurements were performed in non-polar solvents. The yields computed with the original model were 59% and 37% for $c \to t$ and $t \to c$, respectively. Hence the main source of error seems to be the model, rather than SF-TDDFT.

Solvent effects may also affect the yield. These effects can be decomposed into a systematic term and a non-systematic term. It has been argued that non-polar solvents systematically reduce the quantum yield relative to vacuum [Tiberio et al., 2010]. However, careful analysis of the experimental data reported in [Granucci and Persico, 2007] and cited in [Tiberio et al., 2010] does not support this conclusion. Indeed, given the good agreement between non-polar experimental results and both SF-TDDFT  and hh-TDA DFT [Yu et al., 2020], systematic effects of non-polar solvents are likely to be small. Solvent-specific effects are part of the range of reported experimental values, since different works often use different non-polar solvents. Typical ranges were discussed in Sec. <ref>.

A final source of error is the approximations in surface hopping. A large body of literature has examined surface hopping's accuracy; see, for example, Refs. [Landry and Subotnik, 2012, Chen and Reichman, 2016, Zobel et al., 2021, Ibele et al., 2021]. In this work, we can roughly evaluate its performance by comparing its results to those of ab initio multiple spawning (AIMS) [Yu et al., 2020]. AIMS is a fully quantum mechanical method and can thus be treated as a benchmark. We computed the $t \to c$ yield as 26 $\pm$ 3% using surface hopping with SF-TDDFT, and Ref. [Yu et al., 2020] computed the yield as 24% $\pm$ 6% using AIMS with hh-TDA DFT. While there may be some error cancellation from the different quantum chemistry methods, the good agreement is still encouraging.

§ INFLUENCE OF DIFFERENT FUNCTIONAL GROUPS

Here we discuss how the model transferability depends on the functional groups in a compound. To answer this question, we analyzed each of the 40 unseen species in the test set, and computed the model error for each geometry. As described in the main text, the geometries were taken from DANN-NAMD with the trained model; half were selected by proximity to a CI and half selected randomly. For each functional group, we aggregated the errors of all geometries that contained the group, and then computed the mean.

The functional groups with the lowest gap errors are shown in Fig. <ref>(a), and those with the highest are shown in panel (b). Table <ref> shows the number of times that each functional group appears in the training set, together with the errors for all properties. Of the groups with the lowest errors, we see that both B and C are well-represented in the training set, which may explain their accurate results. The other groups each have about 1,000 samples in the training set. This is smaller than that of B and C, but not negligible. The compounds are also rather simple, which may partly explain their low error.

Of the groups with the highest gap error, only the nitro group G is well-represented in the training set. There are almost 10,000 training geometries with nitro groups, yet the error is still quite high. This may be because $\mathrm{NO_2}$ is a strong electron-withdrawing group, and can thus have a significant effect on the electronic structure of azobenzene. Of the remaining substituents, both H and J are rather complicated, and each has only about 200 samples in the training set. Groups F and I both have about 1,000 samples in the training set, and are thus moderately represented. I may have large errors because of the electron-withdrawing effects of the three fluorines, or because it is found together with J in compound 20. F, the tert-butyl group, likely has large errors simply because it is bulky and thus leads to distorted geometries.

We see that the transferability of the model depends on how well-represented a functional group is in the training set, and how complicated its electronic effects are. The former is supported by Fig. <ref>, which shows that the functional group error is anti-correlated with its prevalence in the training set. However, the effect is stronger for the gap than the excited state forces ($\rho = -0.35$ and $\rho=-0.15$, respectively), and neither fully explains the error. The remaining error can best be explained by analyzing the groups themselves.

Functional groups in the test set of unseen species. (a) The five functional groups with the lowest gap error. Error increases from left to right. C is fused to a benzene ring in azobenzene. (b) The five groups with the highest gap error. Error decreases from left to right.

Test set error by functional group, plotted against the number of times the group appears in the training set. The Spearman rank correlation $\rho$ is also shown. (a) Gap error. (b) Excited state force error.

Num. train     $E_0$     $E_1$     $\Delta E_{01}$     $(\Delta E)_{\mathrm{small}}$     sgn$(\Delta E)_{\mathrm{small}} $     $\vec{F}_0$     $\vec{F}_1$     $\vec{g}_{01}$

A       1,219     1.05     0.94     0.63     0.61     -0.04      1.24     1.32     1.19

B      71,940     0.95     1.02     0.81     0.52     -0.01     1.33     1.66     1.62

C      17,624     1.08     1.14     0.85     0.37     0.12      1.24     1.52     1.62

D      994     0.73     1.00     0.91     0.44     -0.07     1.18     1.70     1.54

E      1,154     1.09     1.21     0.95     0.72     -0.04     1.61     2.42     1.41

F      979     2.28     3.36     3.39     1.14     -0.31     1.53     1.53     1.96

G      9,578     0.94     3.49     2.90     0.68     -0.16     1.59      4.11     1.83

H      200     1.62     3.17     2.18     0.44     -0.13     0.99     1.67      0.68

I      1,225     1.33     2.90     2.10     0.78     -0.17     1.63     2.90     0.91

J      236     1.33     2.90     2.10     0.78     -0.17     1.63     2.90     0.91

Test set error by functional group. The groups are divided into those with the lowest gap errors (top) and those with the highest (bottom). “Num. train” refers to the number of geometries in the training set that contain the functional group. While each of the functional groups is contained in the training set, the precise combinations and arrangements of groups (and hence the molecular graphs) are unique to the test set.

§ SPIN CONTAMINATION

Analysis of spin contamination in the training set. (a) Distribution of $\langle S^2 \rangle$ in the ground state. (b) As in (a), but for the excited state. (c) Cumulative probability as a function of $\langle S^2 \rangle$ in the ground state. (d) As in (c), but for the excited state.

Here we discuss the amount of spin contamination in the training set. Figure <ref> shows the square spin, $\langle S^2 \rangle$, for both the ground and excited states. The spin contamination is quite low for the ground state, with $\langle S^2 \rangle$ under 1.0 for 96% of geometries. It is much higher for the excited state, with 76%, 82%, and 93% of geometries having $\langle S^2 \rangle$ under 1.0, 1.2, and 1.4, respectively. Most of the geometries with excited-state $\langle S^2 \rangle$ $\gtrapprox 1.2$ were near the non-reactive $S_1/S_2$ CI, characterized by near-planarity and CNN angles near 108$^\circ$ [Yu et al., 2020]. The $S_0/S_1$ CIs, on the other hand, did not have severe spin contamination. The spin-contamination near the $S_1/S_2$ CI led to difficulty in fine-tuning the model for unsubstituted azobenzene. We also visually inspected a sample of geometries with extreme spin contamination, $\langle S^2 \rangle > 1.8$, and found that nearly all had broken apart. We kept these geometries in the training set so that the model could learn the high energy of bond breaking.

§ SURFACE HOPPING RESULTS WITH DIFFERENT MODELS

Figure <ref> compares ZN hopping statistics of adiabatic and diabatic models for unseen species. Panel (a) shows the distribution of hopping percentages among species. The hopping percentage is defined as the percent of trajectories for a given species that end in the ground state. The $y$-axis is the percent of all species that correspond to each bin. Panel (b) has an analogous plot, but with the hopping percentage replaced by the lifetime. The lifetime was estimated from an exponential fit of $S_1$ population, $p =  \mathrm{exp}(-(t-t_{\mathrm{on}})/\tau) \ \Theta(t-t_{\mathrm{on}}) + \Theta(t_{\mathrm{on}} - t )$. Here $p  \in  [0, 1]$ is the $S_1$ population, $t$ is the time, $t_{\mathrm{on}}$ is the fitted turn-on time, $\tau$ is the fitted lifetime, and $\Theta$ is the Heaviside step function. Trajectories that did not contain any hops were assigned the maximum lifetime of all the other trajectories.

The cis derivatives have lifetimes around 50 fs, consistent with computational results for cis azobenzene [Yue et al., 2018]. Nearly 100% of all cis trajectories ended in the ground state. The trans derivatives have a wide distribution of lifetimes. Some are between 1 and 2 ps, which is similar to trans azobenzene [Yue et al., 2018]. Others are in the range of tens to hundreds of ps, which reflects the high proportion of trajectories that never returned to the ground state. These are very likely incorrect. They may be because of barriers between $S_1$ minima and $S_0/S_1$ CIs, which are known to exist for trans azobenzene [Yu et al., 2020]. Relatively small errors in the barrier may lead to large over-estimations of the lifetime.

Excited state barriers make it even more important to have accurate PESs near CIs. Since trajectories spend little time near crossing regions, trapping becomes even more severe when the gap near CIs is overestimated. The diabatic model helps to address this problem: in each plot we can see that the diabatic model leads to more hopping. For example, using the adiabatic model, 21% of species have hopping percentages under 10%. This number is reduced to 13% for the diabatic model.

The diabatic model also improves the quantum yield. Figure <ref> shows the ZN quantum yields with the diabatic model, and Fig. <ref> shows the diabatic FS yields. The Spearman rank correlation for the trans species is fairly high with both the ZN and FS methods, with the model accurately predicting low yields for a number of species. Figure <ref> shows the ZN yield with the adiabatic model. The correlation among all species is reasonable, but for trans species is rather low. This is because of three molecules with zero predicted yield, all of which became trapped in the excited state. The diabatic model only predicts zero yield for one of these. Further, the diabatic model properly predicts zero yield for species 35, while the adiabatic model does not. For these reasons the diabatic model has a fairly high correlation with experiment. Still, it is clear that excited state trapping of trans species is not a fully solved problem. Preferential sampling of excited state barriers may help to further address this issue in the future.

Comparison of ZN hopping statistics in diabatic and adiabatic models. (a) Proportion of trajectories in each species that hopped to the ground state. (b) Excited state lifetime for each species.

Experimental vs. predicted ZN yields using the diabatic model. (a) Cis isomers. (b) Trans isomers. (c) All species.

Experimental vs. predicted ZN yields using the adiabatic model. (a) Cis isomers. (b) Trans isomers. (c) All species.

Experimental vs. predicted FS yields using the diabatic model. (a) Cis isomers. (b) Trans isomers. (c) All species.

§ ARCHITECTURE

Hyperparameter     Meaning     Value or name

$F$     dimension of hidden atomic features     128

$n_{\mathrm{conv}}$     number of convolutions     5

$n_{\mathrm{RBF}}$     number of radial basis functions (RBF)     20

$R_{\mathrm{cut}}$     cutoff distance for convolutions     5.0 Å

activation function     activation used in message-passing and readout     Swish

learnable $k$     whether $k$ parameters in RBF are learnable     true

skip     output is sum of per-convolution outputs     true only for adiabatic

$M_{\mathrm{d}}$     number of diabatic states     3

$M_{\mathrm{ad}}$     number of ground truth adiabatic states     2

Model hyperparameters

All models were implemented in PyTorch [Paszke et al., 2019]. The model hyperparameters are given in Table <ref>, and an in-depth explanation of the PaiNN parameters can be found in Ref. [Schütt et al., 2021]. Note that we used five convolutions instead of the three used originally, as this substantially improved model performance.
Following DimeNet [Klicpera et al., 2020], we allowed the $k$ values in the radial basis functions to be updated during training. For the adiabatic model, we predicted each property as a sum over per-convolution properties, which was also used in DimeNet. In particular, each convolution had a readout network to convert the atomic features to an output. The final property was obtained by summing each of these outputs.

Several variations on the architecture were tested. For example, we trained both two- and three-state diabatic models on 5,000 azobenzene configurations. We found that adding a third diabatic state significantly decreased the error for all properties. We then trained models using all possible three-state reference orderings and chose the best one. We also experimented with intensive pooling for off-diagonal energies (see Supplementary Sec. <ref>). In particular, we generated a molecular fingerprint through an attention-weighted average of atomic fingerprints. We then mapped this fingerprint to the $d_{nm}$ for $n \neq m$. Even though the $d_{nm}$ are intensive for $n\neq m$, this approach did not improve model predictions.

We also experimented with several adiabatic models. The model in the main text predicted $E_0$ and $E_1$ directly. This approach was used in previous work for single-molecule non-adiabatic dynamics [Chen et al., 2018, Dral et al., 2018, Hu et al., 2018, Li et al., 2021, Westermayr et al., 2019, Westermayr et al., 2020, Westermayr and Marquetand, 2020] and for the prediction of absorption spectra across chemical space [Westermayr and Marquetand, 2020]. We also examined three models that predicted $E_0$ directly and $E_1$ as the sum of $E_0$ and a learned gap. We tried three different pooling methods for the learned gap $\Delta E$: taking the mean over atomwise gaps, taking an attention-weighted average over atomwise gaps, and applying a dense readout network to an attention-weighted sum of atomic fingerprints. All approaches performed quite poorly compared to learning $E_0$ and $E_1$ separately as summed atomwise energies. This is problematic, because only adiabatic models that predict $E_1$ as $E_0 + \Delta E$ can guarantee the positivity of the gap (e.g. by squaring $\Delta E$ or applying a softplus function). Hence the adiabatic model in the main text sometimes generated predictions with $E_1 < E_0$. This is impossible in the diabatic model by construction.

§ TRAINING

$\rho_{E_0}$      $\rho_{E_1}$     $\rho_{\Delta  E_{01}}$     $\rho_{f_0}$     $\rho_{f_1}$     $\rho_{\mathrm{ref}}$     $\rho_{\mathrm{nacv}}$

0.2     0.1     0.5      1     1     0.01     1

Loss parameters used to train all diabatic models.

lr     $\mathrm{lr}_{\mathrm{min}}$     $\rho_{\Delta  E_{01}}$      $\rho_{\mathrm{small}}$

Stage 1      $10^{-4}$     $10^{-5}$     0.5     0

Stage 2     $10^{-5}$     $10^{-6}$      1.0     100

Variable loss parameters used for training the adiabatic model. lr is the starting learning rate and $\mathrm{lr}_{\mathrm{min}}$ is the minimum learning rate, at which point training is stopped.

The training set contained 562,037 geometries from 8,197 species. 5,000 geometries from 308 species were used for validation. The remaining 74,322 geometries from 40 species were held out for testing, so that the predicted yields of unseen species could be compared with experiment.  A different random seed was used to determine the training and validation splits for each committee model, and also to initialize the different models. After training, we ran FS DANN-NAMD on 40 holdout species using the trained diabatic model. For each species we selected 50 geometries randomly and 50 by CI proximity (Eq. (<ref>)), for a total of 4,000 geometries. These geometries were used as the test set, giving the “unseen” statistics in Table <ref>tab:accuracy. The DANN-NAMD geometries from the diabatic network were used for all models, including the adiabatic and ablated ones in Table <ref>. The “seen” statistics were generated using the validation set. For all statistics, a phase correction was applied to $\vec{g}_{01}$ to minimize the prediction error, as in Eq. (<ref>).

The model for screening new azobenzene derivatives was trained on all species. Inclusion of the holdout species provided an additional 74,322 geometries. We used 631,367 geometries for training, 5,000 for validation, and 5,000 for testing. This corresponded to 8,215 training species, 332 validation species, and 303 test species.

Training was performed over energies and forces/force couplings in units of kcal/mol and kcal/mol/Å, respectively. Per-species reference energies were subtracted from each energy. These were obtained by summing atomic reference energies, computed using multi-variable linear regression from (atom type, count) to relaxed geometry energy. Configurations with 10-$\sigma$ energy and force outliers were removed prior to training. Those with forces or energies $\geq$450 kcal/mol(/Å) from the mean were also removed. We found that more stringent removal of outliers led to less stable trajectories. For example, removing 3-$\sigma$ outliers and maximal energies/forces of $\geq$300 kcal/mol(/Å) led to unstable ground state trajectories for six of the 40 unseen species. The 10-$\sigma$ and 450 kcal/mol(/Å) criteria led to no diverging ground state trajectories. A small proportion of excited state trajectories were still unstable. We discarded all excited state trajectories that produced NaN geometries or energies, which were at most a few percent of the trajectories for a given species.

Models were trained with the Adam algorithm using a batch size of 60. Geometries were sampled for each batch using Eq. (<ref>), as described below. The loss was given by Eq. (<ref>eq:overall_loss), using parameters in Table <ref>. Note that the range of NACVs is approximately 10 times smaller than the range of forces. This means that $\rho_{\mathrm{nacv}}=1$, $\rho_{f_i}=1$ gives much higher weight to the forces. We experimented with a range of NACV coefficients between 0.1 and 10, and found that $\rho_{\mathrm{nacv}} = 1$ gave the best performance.

The learning rate was initialized to $10^{-4}$ and reduced by a factor of two if the validation loss had not improved in 10 epochs.
The final model was selected as the one with the lowest validation loss. Training was performed on a single 32 GB Nvidia Volta V100 GPU, and took 13 days to complete.

For diabatic models, the training was stopped once the learning rate reached $10^{-5}$. Adiabatic models were trained in a two-step process, using different loss functions in each stage. Each step ended once the learning rate fell below a certain value. The following loss function was used with different loss trade-offs at different stages:

\begin{align}
& \mathcal{L}  = \mathcal{L}_{\mathrm{core}} + \mathcal{L}_{\mathrm{small}}, \nonumber \\
& \mathcal{L}_{\text{small}} =  \sum_{n > m} \rho_{\mathrm{small}} \cdot  \mathrm{mse}\big( \Delta E^{^{\mathrm{small}}}_{nm} \big),
\end{align}

where $\mathcal{L}_{\text{small}}$ penalizes errors in gaps under 0.2 eV. The parameters for each stage are given in Table <ref>, and the $\rho_{E}$ and $\rho_f$ coefficients are the same as in Table <ref>. The first stage emphasized energy gaps and gradients, while the second stage emphasized small gaps to fine-tune the model near conical intersections. We also experimented with scheduled training and using $\mathcal{L}_{\mathrm{small}}$ for diabatic models, but did not find any improvements.

§ BALANCED SAMPLING

A custom data sampler was used during training because the dataset was imbalanced in the following ways. First, the combinatorial species only had a few geometries each, and so would be very rarely sampled during training. Second, there were more equilibrium trans geometries than cis geometries. Third, there were more equilibrium geometries in general than near-CI configurations. Our sampler addressed these imbalances by giving higher sampling probability to underrepresented species and/or configurations.

Sampling probabilities for different geometries in the training set, with and without the custom sampler. (a) Cumulative probability $p(i)$ of sampling any species $\leq i$. The species are ordered from fewest to most geometries. (b) Probability of sampling a geometry in each of the three clusters. (c) Probability density of sampling a geometry with a given energy gap.

To explain the sampling procedure, let us define two types of sampling weights. Sampling weights that are balanced by species are denoted by $w$, and those that are not are denoted by $v$. For example, the weights $w_{\mathrm{cluster}}(g_{i, A})$ and $v_{\mathrm{cluster}}(g_{i, A})$ are both assigned to the $i^{\mathrm{th}}$ geometry in species $A$, denoted $g_{i, A}$, based on the cluster of configurations that it belongs to. This cluster is denoted $k$, where $k \in$ (cis, trans, other). For a geometry $g_{i, A}$ in cluster $k$, the weights are given as follows:

w_cluster(g_i,A) =  ( 1 /n_spec  ·n_clusters  ) ( 1/n_k, A  )

v_cluster(g_i, A) = ( 1 /n_geoms ·n_clusters  ) ( n_A/n_k, A )

Here $n_{k, A}$ is the number of geometries of species $A$ in cluster $k$, and $n_A$ is the number of geometries in species $A$. $n_{\mathrm{spec}}$ is the total number of species, $n_{\mathrm{clusters}}$ is the total number of clusters, and $n_{\mathrm{geoms}}$ is the total number of geometries. The above definitions ensure that there is an equal probability $P_k$ of sampling any cluster $k$:

\begin{align}
& \text{Using }w: \ P_k = \sum_A \sum_{g_{i, A} \in \ k} w_{\mathrm{cluster}}(g_{i, A}) = \sum_A \frac{1}{n_{\mathrm{spec}} \cdot n_{\mathrm{clusters}}  } = \frac{1}{n_{\mathrm{clusters}} } \\
& \text{Using }v: \ P_k = \sum_A \sum_{g_{i, A} \in \ k} v_{\mathrm{cluster}}(g_{i, A}) = \sum_A \frac{n_A}{n_{\mathrm{geoms}} \cdot n_{\mathrm{clusters}} } = \frac{1}{n_{\mathrm{clusters}} }
% =  \left( \frac{n_{k, A} \cdot  n_{\mathrm{spec}} \cdot  n_{k, A}^{-1}}{n_{\mathrm{spec}} \cdot \sum_j n_{j, A}^{-1} }\right)   =  \frac{1}{ \sum_j n_{j, A}^{-1} } = \mathrm{const} \\
% & \text{Using }v: \ P_k = \sum_A \sum_{g_i \in \ k, A} v_{\mathrm{cluster}}(g_i) =  \left( \frac{n_{k, A} \cdot  n_{\mathrm{spec}} \cdot  n_{k, A}^{-1}}{n_{\mathrm{spec}} \cdot \sum_{j, B} n_{j, B}^{-1} }\right)
\end{align}

Here we have used the fact that $\sum_{g_i \in \ k, A} \ 1 = n_{k, A}$, $\sum_{A} \ 1 = n_{\mathrm{spec}}$, and $\sum_A n_A = n_{\mathrm{geoms}}$. The difference between the two weights is that $w$ leads to balanced sampling among species, while $v$ does not:

\begin{align}
& \text{Using }w: \ P_A = \sum_k \sum_{g_{i, A} \in k} w_{\mathrm{cluster}}(g_{i, A}) = \sum_k \ \frac{1}{n_{\mathrm{spec}} \cdot n_{\mathrm{clusters}}  } = \frac{1}{n_{\mathrm{spec}}} \\
& \text{Using }v: \ P_A = \sum_k \sum_{g_{i, A} \in k} w_{\mathrm{cluster}}(g_{i, A}) = \sum_k \ \frac{n_A}{n_{\mathrm{geoms}} \cdot n_{\mathrm{clusters}}  } = \frac{n_A}{n_{\mathrm{geoms}}}.
\end{align}

Here we have used the fact that $\sum_k \ 1 = n_{\mathrm{cluster}}$. We see that $w$ gives equal weights for each species, whereas $v$ gives weights proportional to the number of geometries in that species.

Similarly to the cluster weights, we define Zhu-Nakamura weights $w_{\mathrm{ZN}}$ and $v_{\mathrm{ZN}}$, such that geometries with a higher hopping probability $p_{\mathrm{ZN}}$ are sampled more often. The corresponding expressions are:

w_ZN(g_i, A) = ( 1/n_spec ·∑_j p_ZN(g_j, A) ) p_ZN(g_i, A)

v_ZN(g_i, A) = ( 1 /∑_j, B p_ZN(g_j, B) ) p_ZN(g_i, A).

The overall sampling weight for a given geometry is determined by each $w$ and $v$, together with two user-defined weights: $P_{\mathrm{spec}} \in [0, 1]$, the importance of species balance, and $P_{\mathrm{ZN}} \in [0, 1]$, the importance of sampling geometries with high hopping rates. $P_{\mathrm{cluster}} = 1 - P_{\mathrm{ZN}}$ is the importance of sampling different clusters in a balanced way. The final sampling weight is then given by

p(g_i, A)     =  P_ZN [ P_spec ·w_ZN(g_i, A) + (1 - P_spec) ·v_ZN(g_i, A) ]

+ P_cluster [ P_spec ·w_cluster(g_i, A) + (1 - P_spec) ·v_cluster(g_i, A) ].

In applying Eq. (<ref>) we used $P_{\mathrm{spec}} = 0.6$ and $P_{\mathrm{cluster}} = 0.5$. We defined a geometry as cis or trans if its RMSD with respect to the corresponding reference structure was $\leq 0.25$ Å, and “other” otherwise. In contrast to diabatic reference geometries, the RMSD was computed using all atoms, not just the central CNNC atoms. Note that while our clustering approach was specific to cis-trans isomerization, many black-box clustering methods exist, such as hierarchical clustering [Schwalbe-Koda et al., 2021]. Any of these methods could have been used in place of our user-defined clusters.

In principle $p_{\mathrm{ZN}}$ should depend on trajectory-specific factors such as velocity, and is therefore not a function of geometry alone. However, during simulations, hops almost always occurred below 0.5 eV. Further, since $p_{\mathrm{ZN}}$ is approximately exponential in the square of the gap, we approximated it with

\begin{align}
p_{\mathrm{ZN}}(\Delta E) \approx \mathrm{exp}(-\Delta E ^2 / (2 \Delta E_0^2)), \label{sm_eq:p_zn}
\end{align}

where $\Delta E$ is the energy gap and $\Delta E_0 = 0.15$ eV. This choice of $\Delta E_0$ gave $p(\Delta E = 0.5 \text{ eV}) \approx 3\times 10^{-3}$. This meant that geometries were usually selected only if their gap was under 0.5 eV.

The sampling probabilities for the training set are shown in Fig. <ref>. Panel (a) shows the cumulative probability of sampling different species, with the compounds ordered from fewest to most geometries. Without the custom sampler the probability was quite small for most molecules. This was because the majority were combinatorially generated and had few geometries. 90% of the probability was contained in the last 41 species. This probability was reduced to 35% when using the balanced sampler. Panel (b) shows that most geometries were trans isomers and few were cis. The custom sampler gave approximately equal probability to the two isomers, and the highest probability to the other group. Most geometries in this group had small gaps, and hence were highly weighted by $w_{\mathrm{ZN}}$ and $v_{\mathrm{ZN}}$. This is demonstrated in panel (c), which shows that near-CI geometries had the highest probabilities when using the balanced sampler.

§ DYNAMICS

Simulations were performed in two stages. First the starting geometries and velocities were generated for each NAMD trajectory. For ab initio simulations we optimized the ground state geometry of each species and computed its normal modes. These modes were used to sample geometries and velocities from the harmonic oscillator Wigner distribution at temperature $T$ = 300 K [Barbatti and Sen, 2016]. For neural trajectories we ran classical MD for 15 ps with the Nosé-Hoover thermostat [Nosé, 1984, Hoover, 1985], and selected random samples to start the NAMD trajectories. Classical MD was implemented with ASE [Larsen et al., 2017]. The effective mass $Q$ was set to $(3 N - 6) \cdot \tau^2  k_{\mathrm{B}} T $, where $N$ is the number of atoms, $k_{\mathrm{B}}$ is Boltzmann's constant, $T = 300$ K, and $\tau = 25$ fs is the relaxation time. Initial atomic velocities were sampled from a Maxwell-Boltzmann distribution at 300 K. Optimized geometries were used for the initial atomic positions. The total linear and angular momenta of the system were set to zero, and the target kinetic energy was set to $ (3N -6) \cdot k_{\mathrm{B}} T/ 2$. The time step of the simulation was set to 0.5 fs, and the neighbor list of the system was updated every 10 time steps (5 fs). All pairs of atoms within with 7 Å were considered neighbors. In each step the distance between all pairs of neighbors was computed, and only pairs within 5 Å of each other were used in the model. This procedure meant that atoms entering the 5 Å cutoff between neighbor list updates would not be missed (i.e., an extra 2 Å “skin” was addded). All ground state and excited state trajectories in this work were propagated with the velocity Verlet algorithm [Verlet, 1967]. Frames for NAMD were taken only after the first 1 ps of ground state MD to allow for equilibration.

Next we ran ZN dynamics [Yu et al., 2014, Yue et al., 2017, Yue et al., 2018] using $N_{\mathrm{trj}}$ trajectories, where $N_{\mathrm{trj}}$ was 10 during the active learning cycle and 500 during final inference.
Hops were restricted to gaps $\leq 0.5$ eV to avoid unphysical transitions, the time step was set to 0.5 fs [Yu et al., 2014, Yue et al., 2017, Yue et al., 2018], and the neighbor list was again updated every 10 steps. Excited state trajectories were run for 1.5 ps during active learning and 5 ps for final inference. All other details of the implementation can be found in Refs. [Yu et al., 2014, Yue et al., 2017, Yue et al., 2018]. For each species we performed inference in parallel over 100-250 geometries at a time, depending on the number of atoms in the molecule. This was done by batching together geometries from 100-250 trajectories and evaluating the model on the batch. The batch size depended on the size of the molecule, which determined the GPU memory consumption. FS simulations were also performed for final inference. The ZN parameters above were used for FS DANN-NAMD of molecules in the holdout set. For virtual screening we used 100 trajectories and 2 ps of excited state dynamics. Species with a hopping rate under 10%, or with molecular graphs that changed during ground state dynamics, were deemed unreliable. These molecules were excluded from Fig. <ref>fig:screening, but not from the average hopping percentage reported in the main text.

The quantum yield was calculated as $Y = n_R / n_T$, where $n_R$ is the number of reactive trajectories and $n_T$ is the total number of trajectories. The uncertainty was computed as the standard deviation of 1,000 bootstrapped samples. To identify reactivity, we computed the RMSD between the CNNC atoms in the final frame and the corresponding atoms in the optimized cis/trans geometries. A trajectory was considered reactive if it started as cis and ended closer to trans, or vice-versa. Trajectories that ended in the excited state were excluded from $n_R$ and $n_T$. The yield was reported as $0 \pm 0$ if all trajectories ended in the excited state.

§ FEWEST SWITCHES IMPLEMENTATION

Tully's surface hopping propagates the nuclei on one electronic surface at a time, called the active surface. The expansion coefficients of the electronic wavefunction are also propagated, and are used to make stochastic changes in the active surface. The coefficient vector $\mathbf{c}$, expressed in a basis denoted by “rep”, evolves as [Mai et al., 2015]

\begin{align}
\frac{d}{d t} \mathbf{c}^{\mathrm{rep}} = - \frac{i}{\hbar} [\mathbf{H}^{\mathrm{rep}} - i \mathbf{T}^{\mathrm{rep}}] \mathbf{c}^{\mathrm{rep}}. \label{sm_eq:c_eqn}
\end{align}

The Hamiltonian matrix is $\mathbf{H}^{\mathrm{rep}} = \bra{\psi^{\mathrm{rep}}_{n} } H \ket{\psi^{\mathrm{rep}}_{m} }$, and the coupling term is

\begin{align}
(\mathbf{T}^{\mathrm{rep}})_{nm} = \left\langle  \psi^{\mathrm{rep}}_{n}(\vec{R}(t)) \left\vert \frac{\partial}{\partial t} \right\vert \psi^{\mathrm{rep}}_{m} (\vec{R}(t))  \right\rangle   = \vec{v} \cdot \vec{k}^{\mathrm{rep}}_{nm}, \label{sm_eq:T}
\end{align}

where $\vec{v}$ is the classical velocity of the nuclei. In principle the nuclei can be propagated on a PES in any electronic basis, and hops can be performed between states in that basis. For example, nuclei could be propagated on diabatic PESs and hops could occur between diabatic states. However, it is well-established that surface hopping in the adiabatic basis gives the best results [Richter et al., 2011]. Hence the nuclei should be propagated in the adiabatic basis, and hops should be decided using $\mathbf{c}^{\mathrm{ad}}$.

While the adiabatic basis should be used for hopping, $\mathbf{c}$ can be propagated in the diabatic basis, and then transformed into the adiabatic basis for all subsequent manipulations (e.g. deciding on hops, adding decoherence, etc.). Indeed, using Eq. (<ref>) in the adiabatic basis can require small time steps for so-called trivial crossings, where derivative NACVs are quite narrow and large [Plasser et al., 2012]. Since global diabatic states are almost never available in ab initio simulations, a local diabatization method is often used to propagate $\mathbf{c}$ [Plasser et al., 2012, Mai et al., 2018]. This method is very stable and can be used with a fairly large time step.
In this work we have a global diabatic basis that can replace local diabatization. Following Refs. [Plasser et al., 2012, Mai et al., 2015], we then propagate $\mathbf{c}$ as follows:

\begin{align}
& \mathbf{c}^{\mathrm{ad}}(t + \Delta t) = \mathbf{P}^{\mathrm{ad}}(t, t+ \Delta t) \hspace*{0.05cm} \mathbf{c}^{\mathrm{ad}}(t), \label{sm_eq:p_rep} \\
& \mathbf{P}^{\mathrm{ad}}(t, t+ \Delta t) = \mathbf{U}^{\dagger}(t + \Delta t) \hspace*{0.05cm} \mathbf{P}^{\mathrm{d}}(t, t+ \Delta t)  \label{sm_eq:p_di_to_ad} \hspace*{0.05cm} \mathbf{U}(t) \\
&  \mathbf{P}^{\mathrm{d}}(t, t+ \Delta t)  = \Pi_{k=1}^{K}  \ \mathrm{exp} \left[ -i \mathbf{H}^{\mathrm{d}}_k \Delta t / K \right]  \label{sm_eq:p_d} \\
& \mathbf{H}_k^{\mathrm{d}} = \mathbf{H}^{\mathrm{d}}(t) + \frac{k}{K} \left[ \mathbf{H}^{\mathrm{d}}(t + \Delta t) - \mathbf{H}^{\mathrm{d}}(t) \right]. \label{sm_eq:p_interp}
\end{align}

Here $\mathbf{P}^{\mathrm{ad}}(t, t + \Delta t)$ is the propagator of $\mathbf{c}^{\mathrm{ad}}$ from $t$ to $t + \Delta t$ (Eq. (<ref>)), and $\Delta t$ is the timestep. The adiabatic propagator is a transformation of the diabatic propagator into the adiabatic basis (Eq. (<ref>)), using the transformation matrix $\mathbf{U}$ (Eq. (<ref>eq:u_def)). The diabatic propagator is a matrix product of $K$ sub-propagators (Eq. (<ref>)), where $K$ is the number of electronic substeps for each nuclear step. “exp” denotes a matrix exponential, not element-wise exponentiation. The $k^{\mathrm{th}}$ sub-propagator uses a linear interpolation between $\mathbf{H}^{d}(t)$ and $\mathbf{H}^{d}(t + \Delta t)$ to approximate $\mathbf{H}_d(t + k \Delta t / K)$ (Eq. (<ref>)). The diabatic Hamiltonian is produced by the neural network model. The probability of hopping from state $n$ to $m$ is [Mai et al., 2015]

\begin{align}
p_{n \to m} = \left(1 - \frac{\vert c^{\mathrm{ad}}_n(t + \Delta t) \vert^2}{ \vert c^{\mathrm{ad}}_n(t) \vert^2} \right) \frac{\mathrm{Re} \left[ c_m^{\mathrm{ad}}(t + \Delta t) \left(P^{\mathrm{ad}}_{mn}\right)^* \left( c^{\mathrm{ad}}_n(t)\right)^* \right] }{\vert c^{\mathrm{ad}}_n(t) \vert^2  - \mathrm{Re}\left[ c^{\mathrm{ad}}_n (t + \Delta t) \left(P^{\mathrm{ad}}_{nn}\right)^*  \left( c^{\mathrm{ad}}_n(t)\right)^* \right] }, \label{sm_eq:p_hop}
\end{align}

where $\mathrm{Re}(x)$ is the real part of $x$. Note that while our approach follows that of SHARC [Mai et al., 2015], we have written our own code and made it publicly available [128]. Our repository also contains code for surface-hopping with the ZN method.

In this work we initialized $\mathbf{c}^{\mathrm{ad}}(t=0) = [0, 1, 0]$, used Eqs. (<ref>)-(<ref>) to generate $\mathbf{c}^{\mathrm{ad}}(t + \Delta t)$, and computed the hopping probability with Eq. (<ref>). $\mathbf{c}^{\mathrm{ad}}$ was expressed in a three-state basis because our model used three diabatic states, and hence $\mathbf{U}$ had three dimensions. Since the model was only trained on the first and second adiabatic energies, we set $p=0$ for all transitions to the third state. Hops to state $m$ were performed if [Mai et al., 2015]

\begin{align}
\sum_{i=1}^{m-1} p_{n \to i} < r \leq \sum_{i=1}^{m-1} p_{n \to i} + p_{n \to m},
\end{align}

where $0 \leq r \leq 1$ is a random number. The momentum was rescaled after each hop to conserve the total energy. The momentum was multiplied by a constant factor, rather than being re-scaled in the direction of the NACV, to avoid the overhead of a NACV calculation (see below) [Mai et al., 2018]. If the factor was complex—a so-called “frustrated” hop—then no transition was made [Mai et al., 2018].

We also tested propagation in the adiabatic basis. In this case we used the NACV to evaluate $\mathbf{T}$, and constructed $P^{\mathrm{ad}}$ from the interpolation of $\mathbf{H}^{\mathrm{ad}} -i \mathbf{T}$ [Mai et al., 2015]. Both Eq. (<ref>) and Tully's original hopping expression [Tully, 1990] were used. The momentum was re-scaled in the direction of the NACV [Landry and Subotnik, 2012]. In all cases the results were very similar to those of the diabatic basis. Diabatic propagation was ultimately chosen because of its reported stability [Plasser et al., 2012], and because of its computational efficiency. In particular, the diabatic propagation requires only diabatic energies and one adiabatic gradient. It therefore uses only one forward and one backward pass through the neural network. By contrast, constructing the NACVs requires gradients for each diabatic element (Eq. (<ref>eq:hd_nacv)). For three diabatic states, this corresponds to one forward pass and six backward passes. While shared convolution layers and caching mean that $N$ gradients do not take $N\times$ as long as one gradient, we found that they still added significant time. Hence diabatic propagation was chosen as the most efficient method.

The decoherence correction of Ref. [Granucci and Persico, 2007] was used to counter the over-coherence of Tully's original method. A sign correction was also used to remove random sign changes in the eigenvectors of $\mathbf{H}_d$. The sign correction $A_m$ for the $m^{\mathrm{th}}$ eigenvector $\mathbf{v}_m$ was computed as

\begin{align}
&  A_m = \mathrm{sgn} (S_{n'm}) \\
& S_{nm} = \mathbf{v}_n^{\dagger}(t) \mathbf{v}_{m}(t + \Delta t) = [\mathbf{U}^{\dagger}(t) \mathbf{U}(t + \Delta t)]_{nm} \\
&  n' = \mathrm{argmax}_{n} \{ \vert S_{nm} \vert \}.
\end{align}

The transformation matrix and NACVs were corrected through

\begin{align}
& U_{nm}(t + \Delta t) \to A_m U_{nm}(t + \Delta t) \\
& \vec{k}_{nm}(t + \Delta t) \to A_n A_m \hspace*{0.05cm} \vec{k}_{nm}(t + \Delta t) \\
& \vec{g}_{nm}(t + \Delta t) \to A_n A_m \hspace*{0.05cm} \vec{g}_{nm}(t + \Delta t).
\end{align}

Note also that $S_{nm} \approx \bra{\psi_n^{\mathrm{ad}}(t)} \ket{\psi_m^{\mathrm{ad}}(t + \Delta t)}$:

\begin{align}
\bra{\psi_n^{\mathrm{ad}}(t)} \ket{\psi_m^{\mathrm{ad}}(t + \Delta t)} & = \sum_{ij} U_{in}^*(t) U_{jm}(t + \Delta t) \bra{\psi_i^{\mathrm{d}}(t)} \ket{\psi_j^{\mathrm{d}}(t + \Delta t)} \nonumber \\
& \approx  \sum_{ij} U^{\dagger}_{ni}(t) U_{jm}(t + \Delta t) \bra{\psi_i^{\mathrm{d}}(t)} [(1 -  (\vec{v} \Delta t) \cdot \nabla_R) \ket{\psi_j^{\mathrm{d}}(t)}] \nonumber \\
& = \sum_{ij} U^{\dagger}_{ni}(t) U_{jm}(t + \Delta t) \hspace*{0.05cm} \bra{\psi_n^{\mathrm{d}}(t)} \ket{\psi_m^{\mathrm{d}}(t )}  \nonumber \\
& = [\mathbf{U}^{\dagger}(t) \mathbf{U}(t + \Delta t) ]_{nm} = S_{nm}.
\end{align}

Here we used the fact that there is no derivative coupling between diabatic states, and that the diabatic states at a given time are orthonormal (Eq. (<ref>)).

All simulations were performed with a time step of 0.5 fs and $K=25$ substeps for electronic propagation. Unlike in ZN dynamics, hops were not restricted to gaps under 0.5 eV. We found that restricting hops improved the ZN results but hurt the FS results. For example, with no maximum gap, most neural ZN trajectories starting from trans-azobenzene hopped at $\sim 1.5$ eV. This did not match the results of Ref. [Yue et al., 2018], which used ab initio ZN with the same level of DFT theory. After adding the restriction, most hops occurred at gaps under 0.1 eV, in agreement with Ref. [Yue et al., 2018]. By contrast, most FS hops occurred under 0.1 eV even without a maximum restriction, though some still occurred at large gaps. The lifetime and yield without a maximum also better matched previous calculations.

§ ACTIVE LEARNING

New geometries were selected for QC calculations using two different criteria. In the first three active learning loops, new geometries were chosen based on the prediction variance of two neural networks. Our aim was to select geometries with fairly uncertain predictions. We did not want geometries with extremely high uncertainty, as these usually corresponded to broken graphs that were outside the target learning space for the model. We therefore used a log-normal target distribution for the uncertainty. Target uncertainties were randomly sampled from this distribution, and geometries with the closest variance to the targets were selected. The log-normal probability of obtaining a sample $x$ is given by $P(x) = \mathrm{exp}\left[ -(\mathrm{ln}(x/s) -\mu )^2 / 2\sigma^2 \right]/(x \sigma \sqrt{2 \pi})$, where $s$, $\sigma$ and $\mu$ are positive numbers. Under this distribution, both completely certain and completely uncertain predictions have zero probability, with a peak for predictions with medium uncertainty. The decay of the probability at large uncertainty is slow, such that highly uncertain geometries can still be selected with reasonable probability. We set $\mu=0$, $\sigma=1$, and $s=7$, giving a target distribution with a mode of 2.5 kcal/mol(/Å) and a mean of 11.5 kcal/mol(/Å). Committee variances were computed for both forces and energies for each electronic state, and the largest variance was compared with the target variance from the distribution. For each species we selected 16 geometries from ground state MD and 33 from surface hopping. Only geometries from the 164 literature species were chosen. This led to approximately 8,000 new data points in each active learning cycle. Other, simpler methods of geometry selection are certainly possible; for example, random selection is also known to work quite well [Wang et al., 2020, Ang et al., 2021]. Our approach successfully increased the model quality throughout active learning, but we did not thoroughly compare it to other methods. Such a comparison may be of interest in the future.

In the next two loops we used only azobenzene, with the goal of densely sampling the CI region. Half the geometries were selected randomly from the excited state dynamics, and half were selected based on the gap. Because the model overestimated the gap in uncertain regions, we did not want to simply select geometries with low predicted gap. Rather, in each trajectory we identified avoided crossing geometries, and assigned equal sampling probability to all geometries before and after that crossing, up to a maximum predicted gap of 1.7 eV. An avoided crossing geometry was defined as having a gap that was lower than in the previous and subsequent time step. In this way we aimed to identify regions that could have small gaps, even if the model did not predict them to be so.

§ VALIDATION

§.§ Ab initio NAMD

Ab initio simulations were used for unsubstituted cis and trans azobenzene. Starting configurations and velocities were generated with ab initio MD using Q-Chem to drive the dynamics. Ten MD simulations were performed for each species. Each simulation was initiated with a different geometry produced by ground-state neural network MD. The simulations were run for 10 ps each using a time step of 0.5 fs, with the the BP86 functional [Becke, 1988, Perdew, 1986] and 6-31G* basis [Francl et al., 1982]. Initial velocities were sampled from a thermal distribution at 300 K, with rotation and translation projected out. The nuclei were propagated with the Velocity Verlet algorithm. A Nosé-Hoover chain of length 5, timescale $\tau=25$ fs, and temperature $T=300$ K was used as a thermostat. The Fock extrapolation order was set to six, and twelve Fock matrices from previous steps were used in the extrapolation.

167 and 178 excited state trajectories were generated for trans and cis azobenzene, respectively. Each trajectory was initialized with a different set of coordinates and velocities, randomly sampled from the ten ground state simulations after 1 ps of equilibration. Spin-flip TDDFT was used with the BHHLYP functional [Becke, 1993] and 6-31G* basis. DIIS with geometric direct minimization (GDM) [Van Voorhis and Head-Gordon, 2002] was used to improve SCF convergence. We found this to be critical: with the usual DIIS algorithm, the SCF cycle failed to converge for 36% of trajectories. This usually occurred within the first 200 fs.

The dynamics were run with in-house scripts, which can be found at <https://github.com/learningmatter-mit/NeuralForceField>. We propagated the elements of the electronic wavefunction in the adiabatic basis, and corrected the sign of the force NACV by minimizing its change from the previous step. The momentum was re-scaled in the direction of the NACV after a hop [Landry and Subotnik, 2012]. All other parameters were unchanged from the DANN-NAMD simulations. Cis trajectories were propagated for 200 fs, and trans trajectories for a maximum of 1.5 ps. For the former, we extended trajectories that had hopped within the last 50 fs, or not hopped at all, until at least 50 fs had passed since they hopped. For the latter, we stopped a trajectory if it had lasted at least 500 fs and had hopped more than 300 fs earlier.

§.§ Transfer learning

To validate the yield results for substituted compounds, we performed DANN-NAMD for the top candidates using a set of highly accurate models. Each model was fine-tuned for a single species only, which allowed it to achieve high accuracy on that one molecule. This transfer learning strategy was used in place of ab initio NAMD because the latter would be prohibitively slow for all but the smallest molecules. For example, consider molecule 169, which has only 54 atoms. We found that a single gradient or NACV calculation for this species took approximately 50 minutes with 8 CPU cores. Since trans derivatives must be simulated for at least 1 ps, and since the time step must be no larger than 0.5 fs, we would need to perform 2,000 QC calculations for each trajectory. Assuming parallel calculation of the NACVs and gradients at each step, an ab initio simulation would take 70 days.

By contrast, the fine-tuning approach allows us to perform highly accurate simulations for tens to hundreds of species.

lr     $\mathrm{lr}_{\mathrm{min}}$     patience     factor     batch size     $\rho_{E_0}$      $\rho_{E_1}$     $\rho_{\Delta  E_{01}}$     $\rho_{f_0}$     $\rho_{f_1}$     $\rho_{\mathrm{ref}}$     $\rho_{\mathrm{nacv}}$

$10^{-4}$     $10^{-5}$     50     0.5     1     0.3     0.3     0.5      1     1     0     1

Training parameters used for transfer learning. “Patience” refers to the number of epochs without an improvement in validation loss before the learning rate is reduced. “Factor” is the amount by which the learning rate is lowered.

Each new model was refined from the original network using QC data from a single species. The initial training geometries were sampled from DANN-NAMD simulations with the original DANN model. A committee of three DANN models, each trained on the entire dataset, was used to select the initial geometries (see below). Three fine-tuned models were then trained. Each was initialized with a different random seed and trained with different data splits. 90% of the data was used for training and 10% for validation. The first model was subsequently used for DANN-NAMD. Geometries were selected from these simulations using the newly-trained committee, and each geometry received QC calculations. This data was added to the training set, each committee model was re-trained, and the cycle was repeated as in Fig. <ref>fig:schematics(b). After each cycle we evaluated the model accuracy using the geometries generated by DANN-NAMD.

Initially we selected only 50 geometries in each round of active learning. Once we narrowed our focus to two species, we sampled 500 geometries in each round. The geometries were sampled according to the following four strategies. One third was sampled randomly. One third was chosen by the prediction variance in the excited state forces. Those with the highest variance were selected. We used only the uncertainty in the forces, and not the energies, because the former is a better indicator of trajectory instability [Schwalbe-Koda et al., 2021]. One sixth was chosen by proximity to a CI (Eq. (<ref>) with $\Delta E = 0.2$ eV). The final sixth was chosen to sample excited-state barriers. In particular, we sampled geometries with probability

\begin{align}
p \propto \mathrm{exp} \left[ (E - E_{\mathrm{min}}) / (k_\mathrm{B} T) \right]. \label{sm_eq:p_barrier}
\end{align}

Here $E$ is the excited state energy, $E_{\mathrm{min}}$ is the minimum excited sate energy in the trajectory, and $T=300$ K. We only sampled configurations from 50 fs or later in the simulations. This was done to avoid the initial high-energy geometries encountered before relaxation. Equation (<ref>) is inversely proportional to the room-temperature Boltzmann probability, and thus assigns the highest weight to the highest-energy configurations.

Transfer learning has previously been used to account for solvent effects [Ang et al., 2021] and to reach higher levels of QC theory [Smith et al., 2019, Ang et al., 2021] with only a small portion of the training set. Typically the majority of the network weights are frozen during re-training. This leaves modifiable parameters in only the final few layers. However, we found that the best results were obtained without any parameter freezing. We therefore allowed all parameters to be modified during re-training. Further, we found it best to start with the normal learning rate rather than a reduced one. We also did not use a reference loss, as its effect on the original DANN results was minor, and possibly even harmful. The full set of training parameters can be found in Table <ref>. All trajectory parameters were unchanged from the screening phase. For final inference we performed DANN-NAMD for 5 ps with 2,000 trajectories, using 500 ps of ground state MD.

Sampled by      Metric      $E_0$     $E_1$     $\Delta E_{01}$     $\vec{F}_0$     $\vec{F}_1$     $\vec{g}_{01}$

ZN     MAE ($\downarrow$)      0.48     0.42       0.44     0.83     0.82     0.52

$R^2$ ($\uparrow$)     1.00      1.00     0.96     0.99     0.99      0.86

Barrier     MAE ($\downarrow$)     0.62     0.93     0.86     0.67     0.82     0.66

$R^2$ ($\uparrow$)     0.99     0.99     0.99     1.00     0.99     0.88

Random     MAE ($\downarrow$)     0.74     0.80     0.86     0.71     0.77     0.89

$R^2$ ($\uparrow$)     0.99     0.99     0.98     1.00     1.00     0.88

Uncertainty     MAE ($\downarrow$)     0.74     1.26     1.09     0.83     1.36     1.04

$R^2$ ($\uparrow$)     0.99     0.99     0.99     0.99     0.97     0.73

Test set statistics for the final TL model of species 165. Results are divided by the method used to select samples. “ZN” uses Zhu-Nakamura gap-based sampling [Eq. (<ref>)] and “barrier” uses Eq. (<ref>). 1,244 geometries were used for fine-tuning.

Sampled by      Metric      $E_0$     $E_1$     $\Delta E_{01}$     $\vec{F}_0$     $\vec{F}_1$     $\vec{g}_{01}$

ZN     MAE ($\downarrow$)     0.67     0.54     0.39     0.87     0.84     0.64

$R^2$ ($\uparrow$)     0.98     0.98      0.99     0.99     0.99     0.82

Barrier     MAE ($\downarrow$)     0.89     0.74     0.41     0.61     0.57     0.51

$R^2$ ($\uparrow$)     0.99     0.95     1.00     1.00     1.00     0.98

Random     MAE ($\downarrow$)     0.58     0.66     0.51     0.67     0.70     0.78

$R^2$ ($\uparrow$)     0.99     1.00     1.00     1.00     1.00     0.96

Uncertainty     MAE ($\downarrow$)     0.72     1.66     1.90     0.83     1.53     1.01

$R^2$ ($\uparrow$)     0.99     0.93     0.97     0.99     0.95     0.92

As in Table <ref>, but for species 169. 2,445 geometries were used for fine-tuning.

The accuracy of the fine-tuned networks is shown in Tables <ref> and <ref>. Statistics are shown for 500 geometries that were sampled from DANN-NAMD for each molecule using the final trained networks. The model errors are well below 1 kcal/mol(/Å) for geometries sampled by all methods other than uncertainty. The error on the uncertainty-selected configurations is under 2 kcal/mol(/Å) in all cases. Since these geometries are specifically chosen to have the highest error, and since their errors are still rather small, we can be confident in the accuracy of the models. We note, however, that the average error for the uncertain geometries depends on how many trajectories are run. Running more trajectories means sampling more configurations, and hence finding more geometries with high uncertainty. As mentioned above, we ran 100 trajectories for 5 ps each for both screening and transfer learning (2,000 trajectories were used for the final predictions). We saved frames every 15 fs, leading to 33,333 geometries in total. Since we picked the 167 most uncertain frames, our sample is roughly the same as choosing the top 0.5% of all geometries with the highest error.

One reason that we did not use transfer learning for the unsubstituted compounds was severe spin contamination. This is an artifact of unrestricted SF-TDDFT, and is a well-documented problem in NAMD for trans azobenzene [Yue et al., 2018]. After the first round of active learning, we found that the force error for the unsubstituted models increased to 8 kcal/mol/Å. This error did not drop significantly with more data, even though the new geometries were fairly close to the Franck-Condon region. This issue did not occur with any of the derivatives. We found that the error was correlated with $\langle S^2 \rangle$, and reasoned that our method of singlet selection (Sec. <ref><ref>) was likely choosing triplet excited states. This highlights the issues inherent in SF-TDDFT, and reinforces the need for low-cost, spin-complete alternatives [Filatov, 2015, Lee et al., 2018, Yu et al., 2020].

Lastly, we note that both the ab initio and transfer-learned cis quantum yields were noticeably higher than in other SF-TDDFT studies [Yue et al., 2018]. This is likely because we used MD to initiate the trajectories, rather than normal-mode or Wigner sampling based on the harmonic approximation. Our experiments showed that normal-mode sampling led to decreased cis yields, closer to those reported in Ref. [Yue et al., 2018]. Unlike the trans isomer, cis azobenzene is somewhat flexible, with significant torsions occurring during ground-state MD. This indicates that the harmonic approximation should be avoided when possible. Indeed, using MD ground state sampling together with FS dynamics for cis azobenzene, we obtained a yield of 60 $\pm$ 4%; this is in excellent agreement with experimental values in non-polar solution, which are close to 55% on average [Bandara et al., 2010]. It is in much better agreement than the value of 34% obtained with FS surface hopping in Ref. [Yue et al., 2018]. Since we used the same electronic structure method and the same surface hopping approach, we can be confident that the difference is mainly due to MD sampling.

§.§ Conical intersection pathways

We labeled the CI pathway of each trans trajectory according to its last hopping geometry. Each geometry was compared to two reference planar CIs and two reference rotational CIs. The trajectory was labeled by the reference CI that was the closest to the hopping geometry. That is, if a trajectory hopped at a geometry closer to one of the rotational CIs, it was labeled a rotational trajectory, and similarly for the planar CI.

The reference CIs were chosen from the set of all hopping geometries in the trajectories. The two rotational CIs were those with dihedral angles closest to $90^{\circ}$ and $270^{\circ}$, respectively, and $\mathrm{max}(\alpha_{\mathrm{CNN}}, \alpha_{\mathrm{NNC}})$ closest to $138^{\circ}$. The two planar CIs were those with dihedral angles closest to $186^{\circ}$ and $174^{\circ}$, respectively, and both $\mathrm{max}(\alpha_{\mathrm{CNN}}$ and $\alpha_{\mathrm{NNC}})$ closest to $148^{\circ}$. The angles are those of the optimized CI geometries in Ref. [Yu et al., 2020]. For the derivative 169, we further optimized each of the reference structures to minimize the gap and hence obtain a true CI. For each trajectory we computed the RMSD between the CNNC atoms of the hopping geometry and the CNNC atoms of the reference CIs.

We additionally enforced that a hopping geometry could only be considered planar if $\vert 180 - \theta \vert \leq \theta_0 $, where $\theta$ is the CNNC dihedral angle and $\theta_0$ is a cutoff value. We chose $\theta_0 = 75^{\circ}$ for azobenzene and $65^{\circ}$ for 169. Without this constraint, many of the hopping geometries labeled “planar” in 169 actually led to isomerization. This made a dihedral constraint necessary. On the other hand, when we only labeled all geometries with $\vert \theta - 180 \vert \leq 45^{\circ}$ as planar [Yu et al., 2020], we found through visual inspection that many non-reactive CI geometries were mislabeled as rotational. Using the minimal distance to a reference CI, together with a modest dihedral constraint, led to the most qualitatively reasonable results. We confirmed that none of the geometries labeled as “planar CI” with our metric led to switching, which further reinforced the soundness of our approach.

§ FIGURE DETAILS

§.§ Computational speed-up

Here we describe how the speeds of ML and QC calculations were computed. For QC, we first computed the run time of one gradient calculation on a single geometry, denoted $t_{\mathrm{calc}}$. The node time was then calculated as $t_{\mathrm{node}} = t_{\mathrm{calc}} / n_{\mathrm{calc}}$, with $n_{\mathrm{calc}} = (\text{cores per node}) / (\text{cores per job})$. We assumed 40 cores per node. All QC jobs were performed with 8 cores, and so $n_{\mathrm{calc}}$ was equal to 5.

For ML we performed a batched calculation on ten copies of each geometry, and computed one gradient. The node time was computed as $t_{\mathrm{node}} = t_{\mathrm{calc}} / n_{\mathrm{calc}}$, with $n_{\mathrm{calc}} = 10 \cdot (\text{total memory per gpu}) / (\text{memory used in calculation}) \cdot (\text{GPUs per node})$. We set $(\text{total memory per gpu}) = 32$ GB and $(\text{GPUs per node}) = 2$. We used a script that performed one batched calculation on ten copies of a random geometry, and re-ran it 7,000 times. For each iteration we used the  command with the PID of the current job to access the GPU memory. We re-ran the script many times, rather than running many calculations in one script, because  does not account for all freed memory until a job is finished. That is, after one calculation is finished,  shows that GPU memory is still being occupied by the job. This occurs even after all local variables are deleted. Hence running multiple calculations in one script would yield an overestimated GPU memory.

In the above calculation, we implicitly assumed that multiplying the batch size by $x$ would lead to an $x$-fold increase in memory. In practice we have found that the memory is increased by less than that. This means that the ML speedup in Fig. <ref>fig:results is a conservative estimate.

§.§ Diabatic energies