Department of Materials Science and Engineering, Massachusetts Institute of
Technology, Cambridge, MA, 02139

# Thermal half-lives of azobenzene derivatives: virtual screening based on
intersystem crossing using a machine learning potential

Simon Axelrod  Department of Chemistry and Chemical Biology, Harvard
University, Cambridge, MA, 02138  Eugene Shakhnovich  Department of Chemistry
and Chemical Biology, Harvard University, Cambridge, MA, 02138  Rafael Gómez-
Bombarelli<EMAIL_ADDRESS>Department of Materials Science and Engineering,
Massachusetts Institute of Technology, Cambridge, MA, 02139

###### Abstract

Molecular photoswitches are the foundation of light-activated drugs. A key
photoswitch is azobenzene, which exhibits trans-cis isomerism in response to
light. The thermal half-life of the cis isomer is of crucial importance, since
it controls the duration of the light-induced biological effect. Here we
introduce a computational tool for predicting the thermal half-lives of
azobenzene derivatives. Our automated approach uses a fast and accurate
machine learning potential trained on quantum chemistry data. Building on
well-established earlier evidence, we argue that thermal isomerization
proceeds through rotation mediated by intersystem crossing, and incorporate
this mechanism into our automated workflow. We use our approach to predict the
thermal half-lives of 19,000 azobenzene derivatives. We explore trends and
tradeoffs between barriers and absorption wavelengths, and open-source our
data and software to accelerate research in photopharmacology.

## Introduction

Photoswitches are compounds whose properties can be modified by light. They
have applications in many developing technologies, such as organic electronics
1, energy storage 2, and targeted medicine 3. The latter includes
photopharmacology, the field of light-activated drugs. Such drugs are built
around photoswitchable scaffolds, which allows their medicinal activity to be
controlled with light. The most common scaffold is azobenzene, which undergoes
cis $\leftrightarrow$ trans isomerization in response to light. An inactive
drug built around azobenzene can be activated with light at certain times or
in certain regions of the body. This can minimize off-target activity, thereby
minimizing side effects 3.

Figure 1: Possible mechanisms of thermal isomerization. The inversion TS has
$\alpha$ or $\alpha^{\prime}\approx 180^{\circ}$, while the rotational TS has
$\omega\approx 90^{\circ}$. Inversion-assisted rotation combines inversion and
rotation. Concerted inversion combines $\alpha$ inversion and
$\alpha^{\prime}$ inversion.

The development of photoactive drugs is a complex, multi-objective
optimization problem. A formidable number of properties must be optimized for
all drugs. Photoactive compounds must also absorb light at the right
wavelength, typically in the near infrared region; isomerize with high
efficiency in the excited state; display differential bioactivity between the
two isomers; and thermally revert to the stable isomer in a specific time
frame 4. For many applications this time frame should be as long as possible,
and so the isomerization barrier should be as high as possible 4. Yet
substituents that shift absorption from the UV to the visible or near-IR
regions often lower the thermal barrier 3. This highlights the challenge of
the optimization problem.

The design of photoactive drugs can be accelerated with computational
modeling. Property predictors can be applied to large virtual libraries, and
the results can be used to narrow the search space of promising compounds 5.
Quantum chemistry can predict many properties with good accuracy, but the
calculations are quite slow. In our previous work, we showed how a machine
learning (ML) potential trained on quantum chemistry data can be used to
rapidly predict the quantum yield of azobenzene derivatives 6. In this work we
develop an ML-based computational workflow to predict the thermal half-lives
of azobenzene derivatives.

Our contributions are as follows. First, we improve the theory of thermal
azobenzene isomerization. In particular, we provide evidence that thermal
isomerization proceeds through intersystem crossing, not through a typical
singlet transition state (TS). This builds on the theory that was proposed
nearly twenty years ago 7, but which has largely been overlooked (with some
exceptions 8, 9). We also demonstrate the critical importance of
multireference effects in barrier calculations.

Second, we provide a fast and user-friendly computational tool for predicting
the barriers and absorption wavelengths of azobenzene derivatives. Our program
can be run with a single command. The only user input required is the SMILES
strings of the relevant compounds. These can be generated programatically, or
with visual interface programs such as ChemDraw. Further, the program is quite
fast, replicating the results of spin-flip, time-dependent density functional
theory (SF-TDDFT) 10 in milliseconds through use of a transferable ML
potential. This program adds to the growing collection of computational models
for predicting photoswitch properties 6, 11, 12, 13. Our software and pre-
trained models are freely available at https://github.com/learningmatter-
mit/azo_barriers.

Third, we use our tool to perform virtual screening of nearly 19,000
hypothetical azobenzene derivatives. We identify species with high
isomerization barriers and redshifted absorption spectra. The data is freely
available at https://doi.org/10.18126/unc8-336t through the Materials Data
Facility 14, 15. Researchers can use the species with favorable properties as
scaffolds for new photoactive drugs. Further, we explain these results in
terms of substitution patterns and substituent properties. These insights will
accelerate the design of meta-stable, red-shifted azobenzene derivatives in
the future.

## Theory and methods

### Isomerization mechanisms

Four mechanisms have been proposed for azobenzene isomerization 16: rotation,
inversion, inversion-assisted rotation, and concerted inversion (Fig. 1).
Rotation is characterized by $\omega\approx 90^{\circ}$ and
$\alpha\approx\alpha^{\prime}\approx 120^{\circ}$. Both inversion mechanisms
have $\alpha\approx 180^{\circ}$ and $\alpha^{\prime}\approx 120^{\circ}$;
pure inversion has $\omega\approx 180^{\circ}$, while inversion-assisted
rotation has $\omega\approx 90^{\circ}$. (Some works refer to inversion-
assisted rotation simply as rotation 11, but we avoid that terminology here.)
Concerted inversion has $\alpha\approx\alpha^{\prime}\approx 180^{\circ}$, but
can be excluded from possible thermal mechanisms; see Supporting Information
(SI) Sec. S11.1. Here we group together inversion and inversion-assisted
rotation, and refer to both as inversion (SI Sec. S11.2). For asymmetrically
substituted azobenzenes, $\alpha=180^{\circ}$ is distinct from
$\alpha^{\prime}=180^{\circ}$, and $\omega=90^{\circ}$ is distinct from
$\omega=-90^{\circ}$. This gives two inversion TSs and two rotational TSs.

### Standard models

Several works have predicted the thermal lifetimes of azobenzene derivatives
with computational methods 17, 18, 19, 20, 11, 12. However, there are two main
issues with previous calculations. First, all levels of theory overestimate
the experimental enthalpy and entropy of activation. The experimental
activation entropy is $-50.2$ J $\mathrm{mol}^{-1}\mathrm{K}^{-1}$ 21, while
Hartree-Fock, MP2, CC2, and DFT with 12 different functionals predict values
between $+7$ and $+28$ J $\mathrm{mol}^{-1}\mathrm{K}^{-1}$ 19. At room
temperature, the experimental entropic contribution to $\Delta G^{\dagger}$ is
then $-T\Delta S^{\dagger}=+3.6$ kcal/mol, while the computational
contribution is between $-0.5$ and $-2$ kcal/mol. Most of the error persists
even after corrections to the harmonic approximation 19.

Each method also overestimates $\Delta H^{\dagger}$. The experimental
activation enthalpy is $21.1$ kcal/mol 21; B3LYP-D3/6-311++G** predicts $25.2$
kcal/mol, and CC2/aug-cc-pVTZ predicts 29.8 kcal/mol 19. CASPT2(10,8)/6-31G*
gives $\Delta E^{\dagger}=31.0$ kcal/mol 22, which is nearly identical to
CC2/aug-cc-pVTZ 19. These calculations, like most in the literature, were
performed for the inversion TS. In SI Sec. S12, we show that highly accurate
spin-flip coupled cluster methods yield similar overestimates for the
rotational TS.

The errors in $\Delta H^{\dagger}$ and $-T\Delta S^{\dagger}$ partially cancel
for $\Delta G^{\dagger}$. Indeed, recent work has extensively benchmarked
different levels of theory for predicting $\Delta G^{\dagger}$ of various
azoarenes, and found errors near 1 kcal/mol for B3LYP-D3 with some basis sets
11. However, given the significant error cancellation for $\Delta
G^{\dagger}$, and the underestimation of $\Delta H^{\dagger}$ relative to
correlated wavefunction methods, these results should be interpreted with
caution.

The second issue concerns the rotational mechanism. While most DFT
calculations have been applied to inversion 22, 19, 11, CASPT2 calculations
indicate that rotation is in fact preferred for azobenzene 7. Yet Ref. 19
found that the rotational TS cannot be optimized with B3LYP, a finding that we
have also reproduced for various derivatives. It is troubling that rotation
can be energetically favored, yet the rotational TS cannot be reached through
DFT optimization.

### Intersystem crossing

Figure 2: Schematic of triplet-mediated thermal isomerization. Intersystem
crossing from singlet ($S_{0}$) to triplet ($T_{1}$) occurs near
$\omega=70^{\circ}$. The energy is lower at the crossing point than at the
$S_{0}$ TS ($\Delta E^{\mathrm{X}}<\Delta E^{\dagger}$). Re-crossing from
$T_{1}$ to $S_{0}$ then occurs near $\omega=105^{\circ}$, and isomerization is
completed on the $S_{0}$ surface. Figure 3: Approach to active learning and TS
generation used in this work. (a) Active learning loop for training the NN.
(b) Workflow for generating equilibrium and TS geometries. “confgen” stands
for conformer generation, which is described in SI Sec. S3.2. “$S$-$T$ MECP”
denotes the search for the minimum-energy singlet-triplet crossings on each
side of the TS.

To address the overestimation of $\Delta S^{\dagger}$ and $\Delta H^{\dagger}$
common to all levels of theory, we advocate intersystem crossing (ISC) as the
mechanism of thermal isomerization (Fig. 2). This was proposed in Ref. 7
nearly 20 years ago. As discussed below, $\Delta G^{\dagger}$ is replaced with
$\Delta G^{\mathrm{X}}$, the free energy difference at the singlet-triplet
crossing geometry $\mathrm{X}$. The rate prefactor $k_{\mathrm{B}}T/h$ is
replaced with the ISC rate $k_{\mathrm{ISC}}$. Since the crossing geometry has
a lower energy than the TS, this approach corrects the overestimation of
$\Delta H^{\dagger}$. Further, the measured $\Delta S^{\dagger}$ is not a true
activation entropy, but in fact related to both
$\mathrm{log}[k_{\mathrm{ISC}}/(k_{\mathrm{B}}T/h)]$ and $\Delta
S^{\mathrm{X}}$ (SI Eq. (S28)). The negative value of $\Delta S^{\dagger}$
partly reflects the fact that $k_{\mathrm{ISC}}<k_{\mathrm{B}}T/h$. Computing
$k_{\mathrm{ISC}}$ with CASPT2(14,12)/6-31G* yields good agreement with the
effective experimental activation entropy, $\Delta S^{\mathrm{eff}}$, for
unsubstituted azobenzene 7.

We note, however, that isomerization can still proceed through $S_{0}$ for
certain derivatives and environments, and that the $S_{0}$ rate should be
compared to the triplet-mediated rate for any new system. Indeed, Ref. 23
found good agreement between the $S_{0}$ CASPT2 rate and the experimental rate
for phototrexate in DHFR, though no comparison was made to the experimental
activation entropy. For this reason we compute the rate with both Eyring TS
theory and ISC for all molecules in this work.

### Multireference effects

To properly optimize the rotational TS, we show that it is critical to include
multi-reference effects. This is not surprising, since there is a conical
intersection when $\omega$ is close to $90^{\circ}$ 24. In this work we use
SF-TDDFT 10, since it accounts for some double excitations and generally
provides an accurate description of conical intersections 25. We use the
common BHHLYP functional 26 and 6-31G* basis 27. As shown in SI Sec. S12,
single-reference methods produce rotational TS cusps. This is the likely cause
of the failed rotational TS optimizations with standard DFT. Multi-reference
methods, by contrast, produce smooth maxima. We note that the inversion TS is
also close to a conical intersection 7, which further reinforces the need for
a multi-reference treatment of azobenzene TSs.

SF-TDDFT accounts for multi-reference effects while offering a reasonable
balance between cost and accuracy. Empirically it has sub-cubic scaling 6,
which makes it far more affordable than highly accurate spin-flip coupled
cluster methods that scale as $N^{7}$ 28. Further, as shown in SI Sec. S12,
SF-TDDFT has similar errors to CASPT2 for rotational barriers in azobenzene.
The latter scales as $N^{5}$ for a fixed active space, requires manual active
space selection, and does not have analytic gradients in most quantum
chemistry packages, which are essential for training ML potentials. While SF-
TDDFT is not spin-complete, we have found that $S_{0}$ spin contamination is
rather low. As discussed in SI Sec. S8.1, the average square spin in the
$S_{0}$ state is only 0.16, and we excluded all data with square spin
exceeding 1.0. Lastly, the errors in SF-TDDFT are largely systematic, as
demonstrated by the strong correlation between predicted and experimental
activation free energies in Fig. 4. This is discussed further in the Results
section below.

### Computational workflow

Standard quantum chemistry approaches are rather slow. To address this and the
above issues, we develop an ISC workflow based on ML potentials 29 that are
trained on multi-reference SF-TDDFT data.

We generate initial TSs through a relaxed scan, and refine the structures with
a conformer search 30, 31, 32 and eigenvector following (EVF, Fig. 3(b)). We
then use the intrinsic reaction coordinate (IRC) 33 to locate singlet-triplet
crossings on either side of each rotational TS. The geometries are
subsequently refined with a minimum energy crossing point (MECP) optimization
(SI Sec. S3.5). We also apply Eyring TS theory to the TSs and compare the
results to the ISC approach.

Four relaxed scans are performed for each of the four mechanisms (two
inversion and two rotation). Conformers are generated for each TS guess using
a conformational search with fixed CNNC atoms. The five lowest-energy
conformers for each mechanism are then optimized with eigenvector following,
yielding 20 TSs. The TS with the lowest free energy for each rotational
mechanism is used to find the singlet-triplet crossings. ISC does not occur
for the inversion mechanism 7.

Once the workflow is completed, the ISC-based reaction rate is calculated as

$\displaystyle k_{\mathrm{X}}=k_{\mathrm{ISC}}\hskip
1.9919pt\mathrm{exp}\left(-\Delta G^{\mathrm{X}}/k_{\mathrm{B}}T\right),$ (1)

where $G$ is given in SI Eqs. (S11)-(S15). $\Delta G^{\mathrm{X}}$ depends on
the energies, conformer ensembles, and vibrational frequencies of both the
reactant and the crossing geometry. We compute $k_{\mathrm{ISC}}$ using the
approach in Ref. 34, as described in SI Sec. S4. The expression depends on the
spin-orbit coupling, temperature, and forces on the singlet and triplet
surfaces at the crossing. The spin-orbit coupling is taken as constant,
$H_{\mathrm{SO}}\approx 20\ \mathrm{cm}^{-1}$, as described in SI Sec. S4. The
coupling is about 20 times larger than the typical value for planar aromatic
compounds, because of the $n-\pi^{*}$ character of the triplet state 7. This
corresponds to a 400-fold enhancement in the ISC rate.

We also compute the reaction rate from Eyring TS theory, given by

$\displaystyle
k_{\mathrm{TST}}=\frac{k_{\mathrm{B}}T}{h}\mathrm{exp}\left(-\Delta
G^{\dagger}/k_{\mathrm{B}}T\right).$ (2)

For ease of comparison we convert ISC-based reaction rates into the form of
Eq. (2), where $\Delta G^{\dagger}$ is replaced by the effective activation
free energy,

$\displaystyle\Delta G^{\mathrm{eff}}=-k_{\mathrm{B}}T\
\mathrm{log}\left(\frac{hk_{\mathrm{X}}}{k_{\mathrm{B}}T}\right)$
$\displaystyle=\Delta H^{\mathrm{X}}-T\Delta
S^{\mathrm{eff}}-\frac{1}{2}k_{\mathrm{B}}T.$ (3)

$\Delta S^{\mathrm{eff}}$ depends on both $k_{\mathrm{ISC}}$ and $\Delta
S^{\mathrm{X}}$, and is given in SI Eq. (S28).

Geometry type | Model type |  |  Singlet $\Delta{E}$
---
errora
|  Triplet $\Delta{E}$
---
errorb
|  Singlet $\vec{F}$
---
error
|  Triplet $\Delta\vec{F}$
---
errorc
| One model | 0.81 | 0.23 | 0.44 | 0.48
Optimized TSd |  Ensemblee | 0.73 | 0.19 | 0.34 | 0.44
|  Ensemble, lowest $95\%$ uncertaintyf | 0.66 | 0.16 | 0.31 | 0.42
| One model | 1.09 | 0.34 | 0.69 | 0.55
TS metadynamicsg | Ensemble | 0.99 | 0.31 | 0.54 | 0.50
|  Ensemble, lowest $95\%$ uncertainty | 0.86 | 0.27 | 0.49 | 0.48
Table 1: Model performance for 334 species outside the training set. Units are
kcal/mol for energies and kcal/mol/Å for forces. Forces are denoted by
$\vec{F}$.

### ML models

We train separate models to predict the $S_{0}$ energy, $T_{1}$ energy, and
$S_{0}/S_{1}$ gap. Our models use the PaiNN architecture 35, which predicts
molecular properties through equivariant message-passing. This approach
generates a feature vector for each atom that incorporates information from
its surrounding environment. The initial feature vector is generated from the
atomic number alone, and is then updated through a set of “messages”. These
messages incorporate the distance, orientation, and features of atoms within a
cutoff distance. The messages are then used to update the atomic feature
vectors. This is performed several times, which leads to information being
combined from increasingly distant atoms. Lastly, the atomic features are
mapped to per-atom energies using a neural network, which are summed to yield
the molecular energy. The forces are computed through automatic
differentiation of the energy.

We pre-train the models on 680,736 gas-phase SF/6-31G* calculations from non-
adiabiatc molecular dynamics (NAMD), which were previously generated in Ref.
6. We then refine the models using approximately 40,000 SF/6-31G* calculations
with a C-PCM model of water 36, 37, 38. Pre-training on existing gas-phase
data means that fewer new solvent calculations are required to reach a target
accuracy 39.

The geometries for SF-TDDFT/C-PCM calculations are generated through active
learning 39, 40, 6 based on our TS workflow (Fig. 3). In each round of active
learning, we train three $S_{0}$ models on previous SF-TDDFT/C-PCM data (the
gas-phase model is used in the first round). The difference in model
predictions is used to identify geometries that are poorly described by the
model. These high-uncertainty configurations, together with some geometries
that are sampled randomly or by energy, receive new quantum chemistry
calculations (see SI Sec. S1).

## Results

### Model performance

The model accuracy is shown in Table 1. Mean absolute errors (MAEs) are given
for the singlet and triplet models, for both optimized TSs and off-equilibrium
geometries sampled during TS conformer generation (see SI Sec. S3.2). All
geometries come from species outside the training set.

The model performance is excellent. The error in the barrier energy, $\Delta
E=E_{\mathrm{TS}}-E_{\mathrm{cis}}$, is 0.81 kcal/mol for one model and 0.73
kcal/mol for an ensemble of three models. The ensemble error falls to 0.66
kcal/mol after excluding the top 5% most uncertain geometries. These errors
are far below 1.0 kcal/mol, which is the typical definition of chemical
accuracy. The model error is significantly smaller than the SF-TDDFT error (SI
Sec. S12).

The force predictions are similarly accurate, with MAEs below 0.45 and 0.7
kcal/mol/Å for optimized and distorted TSs, respectively. The performance of
the triplet model is even better, with errors that are four times lower than
that of the singlet model. Lastly, the $S_{0}/S_{1}$ gap model has an MAE of
0.68 and 2.35 kcal/mol for optimized cis and trans geometries, respectively.
These are errors of 3.8 and 13.1 nm for a typical absorption wavelength of 400
nm.

### Comparison to experiment

The predicted and experimental activation free energies are compared in Fig.
4. The experimental data come from Refs. 41, 42, 43, 17, 44, 45, 46, 18, 20,
47, and can be found in the file containing the virtual screening results. 26
measurements were accessed in total, and 17 of these were used in Fig. 4. We
only used the measurements performed in solvent with dielectric constant
$\varepsilon\geq 30$. Our model was trained on implicit solvent calculations
using $\varepsilon=78.4$ for water. However, we found that the quantum
chemistry energies were quite similar when using $\varepsilon=30$, and so
included these measurements in the benchmark as well.

The results of TS and ISC theory are shown in panels (a) and (b),
respectively. The rotation mechanism was favored in TS theory for all
compounds shown, and so panel (a) equivalently shows $\Delta G^{\dagger}$ from
rotation. Panel (c) shows the results of TS theory when considering only the
two inversion mechanisms.

The correlation with experiment is quite good for both TS and ISC theory. Both
have Spearman rank coefficients $\rho$ near 0.6, and both have $R^{2}$ near
0.7 after linear regression. The MAEs after regression are 0.79 and 0.85
kcal/mol for TS and ISC theory, respectively. When including all 26
measurements performed in any solvent, we find that $\rho$ actually increases
to 0.67 for both methods. $R^{2}$ falls to 0.58 and 0.63 for TS and ISC
theory, respectively. The respective MAEs climb to 0.92 and 0.97 kcal/mol.

The species are properly separated into low-, medium-, and high-barrier groups
in Fig. 4. For example, species 1 and 2 are predicted to have low barriers,
the fluorinated derivatives 4-6 to have high barriers, and azobenzene to lie
in the middle. The models even have some success comparing fluorinated
derivatives to each other, with $\rho=0.33$ and $0.42$ among these species for
$\Delta G^{\dagger}$ and $\Delta G^{\mathrm{eff}}$, respectively. However, the
two approaches give $R^{2}=0.02$ and $0.04$, respectively. This means that the
numerical error is close to that of a random predictor, even though the
rankings are better than random.

Figure 4: Experimental vs. predicted activation free energies. Dotted orange
lines are linear regression results from predicted to experimental values.
$\rho$ denotes the Spearman rank correlation. $R^{2}$ is computed between the
regression results and the experimental data. Error bars are the standard
deviation of the energy predictions from three models. (a) Prediction accuracy
using TS theory. (b) Prediction accuracy using intersystem crossing. (c) As in
(a), but with only the inversion mechanism. (d) Selected compounds highlighted
in panels (a)-(c).

Both approaches overestimate the barriers, on average by 4.94 and 3.76
kcal/mol for TS theory and ISC, respectively. However, the overestimation is
largely systematic, as demonstrated by the high $R^{2}$ value after linear
regression. Further, it is a consequence of SF-TDDFT, not the models. Indeed,
as shown in SI Sec. S12, $\Delta H$ is well-reproduced by the accurate and
expensive method SF-EOM-CCSD(dT) 28. SF-TDDFT overestimates $\Delta H$ with
respect to both experiment and SF-EOM-CCSD(dT). Transfer learning to this
higher level of theory could therefore be of interest in the future.

$\Delta G^{\mathrm{eff}}$ is lower than $\Delta G^{\dagger}$ on average, but
otherwise their trends are quite similar. One reason is that rotation is the
predicted mechanism for all species. Since each singlet-triplet crossing is on
either side of a rotational TS, its energy is correlated with that of the TS.
Indeed, the correlation between the two methods is near unity, with
$\rho=0.97$ and $R^{2}=0.98$. This reflects the fact that
$E^{\dagger}_{\mathrm{rot}}-E^{\mathrm{X}}$ and $k_{\mathrm{ISC}}$ are nearly
constant among different species. However, we explain below that noticeable
differences arise when screening large virtual libraries.

The approaches’ strong performance should be contrasted with TS theory using
only inversion. These results are shown in panel (c). The performance is far
worse than using rotation or ISC, with $R^{2}$ reduced by over 50%, and $\rho$
reduced by over $35\%$. The MAEs after regression are 1.24 kcal/mol for polar
solvents and 1.43 kcal/mol for all data. Note that most works with DFT have
only produced inversion TSs. This is likely because of failed rotation
optimizations, which we attribute to the single-reference nature of DFT and
the associated TS cusps (SI Sec. S12). Similar difficulties were found in Ref.
9. Our results highlight the importance of rotation and multi-reference
effects.

Note that all species in Fig. 4 were in the training set. Hence the comparison
to experiment does not measure the model’s ability to generalize to new
compounds. Rather, it mainly measures the reliability of the workflow and the
underlying quantum chemistry. The models’ ability to generalize to new species
was shown in Table 1.

### Virtual screening

With reliable models and predictive workflows, we next screened a virtual
library of azobenzene derivatives for key properties in photopharmacology. We
ran the workflow of Fig. 3(b) for 25,000 compounds in all. The compounds were
generated using the common literature substitution patterns in Fig. 7(d),
folllowing the approach of Refs. 48, 6. The substituents are a combination of
literature groups and basic chemical moieties, and can be found with the
screening results online. After applying various filters, such as the right
number of imaginary frequencies, converged TSs for all mechanisms, and the
proper IRC endpoints, we were left with 19,000 species in total (see SI Sec.
S9).

#### Distributions

Figure 5: Distribution of various quantities among the 19,000 screened
derivatives. (a) Effective activation free energy. (b) Effective activation
entropy. (c) Model uncertainty in the singlet-triplet gap at MECPs. (d)-(e) As
in (a)-(b), but using TS theory. (f) Model uncertainty in the activation
energies.  Figure 6: Absorption wavelengths and thermal barriers. (a) Trans
vs. cis and absorption wavelengths. (b), (c) $\Delta G^{\mathrm{eff}}$ vs. cis
and trans absorption wavelengths, respectively. (d), (e) The two compounds of
interest. For each panel the graph is shown on the left, followed by the cis
geometry, the singlet-triplet crossing closer to cis, the TS, and the trans
geometry. Model and quantum chemistry (QC) predictions are shown below. The
compounds are circled in panels (a)-(c).

Figure 5(a) shows the distribution of $\Delta G^{\mathrm{eff}}$ from this
screen. The mean and median are 29.6 and 29.7 kcal/mol, respectively, while
unsubstituted azobenzene has a value of 28.9 kcal/mol. The average derivative
is thus more kinetically stable than azobenzene. The standard deviation of is
2.6 kcal/mol, which is a factor of 80 in the isomerization rate. 39% of
species (7,400) have a lifetime that is over 10$\times$ that of azobenzene;
19% (3,600) have a lifetime that is less than $1/10^{\mathrm{th}}$. We
conclude that the lifetime is highly tunable using the substitutions in this
work.

Figure 7: Chemical properties by substitution pattern. 3-$\sigma$ outliers
were removed for ease of visualization. (a) Effective activation free energy.
(b) Trans absorption wavelength. (c) Cis absorption wavelength. (d)
Substitution patterns. Figure 8: Relationships among substituent properties,
free energies, and absorption wavelengths for pattern A. The pattern is shown
in Fig. 7 (d). (a) $\lambda_{\mathrm{trans}}$ for donor/acceptor and ring/non-
ring $R_{1}$ substituents. (b) $\lambda_{\mathrm{trans}}$ as a function of
$R_{1}$ and $R_{2}$ donor scores. (c) $\Delta G^{\mathrm{eff}}$ as a function
of $R_{1}$ donor score. (d) $\Delta G^{\mathrm{eff}}$ as a function of both
$R_{1}$ and $R_{2}$ donor scores. (e) As in (d), but for $\Delta G^{\dagger}$.
(f) $\Delta G^{\mathrm{eff}}$ vs. $\Delta G_{\mathrm{rxn}}$.

Panel (b) shows the effective activation entropy. The mean is $-$24.5 J / mol
K and the median is $-$26.9 J / mol K. The calculated and experimental values
for azobenzene are $-$29.0 J / mol K and $-$50.2 J / mol K 21, respectively.
The associated error in the entropic free energy is 1.5 kcal/mol. Using the TS
approach gives $\Delta S^{\dagger}=4.7$ J / mol K, which has a much higher
entropic free energy error of 3.9 kcal/mol. Few other derivatives have
experimental activation entropies, and of those that do, most have values near
that of azobenzene. However, of the few with values far from azobenzene, none
were predicted accurately by the model 49, 20, 47. These cases should be
investigated in more detail in the future.

Panels (d) and (e) are analogous to (a) and (b), but with TS theory instead of
the ISC approach. Each distribution resembles its partner from ISC theory.
However, $\Delta G^{\dagger}$ is more asymmetric than $\Delta
G^{\mathrm{eff}}$, with a much steeper drop-off for higher barriers. The
reason is as follows. Within TS theory using SF-TDDFT in water, rotation is
more often the preferred mechanism (SI Sec. S6). However, inversion can become
preferred for species with high enough rotation barriers. This mechanism is
not available in ISC theory, since $T_{1}$ is always higher than $S_{0}$
during inversion. Hence inversion lowers the high barriers in TS theory, but
cannot do the same in ISC theory.

In principle one should calculate both $\Delta G^{\dagger}$ and $\Delta
G^{\mathrm{eff}}$, and use the lower value for the reaction rate. If $\Delta
G^{\dagger}_{\mathrm{inv}}$ were low enough, it would replace $\Delta
G^{\mathrm{eff}}$ in the high barrier regime, and so the steep drop-off of
panel (d) would be observed. However, as discussed in SI Sec. S6, SF-TDDFT is
not accurate enough to compare absolute $\Delta G^{\dagger}$ and $\Delta
G^{\mathrm{eff}}$ directly. Hence a more accurate treatment of this problem
would be of interest in the future.

Panels (c) and (f) show the model uncertainty in the MECP singlet-triplet gap
and the activation energy, respectively. Both are quite low, indicating high
model confidence for the derivatives studied here. The mean uncertainties are
0.21 kcal/mol and 0.41 kcal/mol for the singlet-triplet gap and activation
energy, respectively. This is consistent with the error trends in Table 1. The
uncertainty should be interpreted with caution, however, as neural network
ensembles tend to be overconfident 50. Large uncertainty necessarily means
high error, as demonstrated by the error reduction in Table 1 when excluding
the most uncertain geometries; however, low uncertainty does not guarantee low
error. A more detailed examination of uncertainty, including calibration to
the observed error and use of different architectures in the ensemble, is left
to future work.

#### Targeting desired properties

Absorption wavelength and thermal stability are two key properties in the
design of photoactive drugs. The preferred absorption range is generally 650
to 900 nm, since human tissue is transparent only in this narrow region of the
near-IR 4. For photoactive drugs one typically wants the isomerization barrier
to be as high as possible, so that the unstable isomer is active for as long
as possible. For ion channel blockers, by contrast, the target lifetime is
usually milliseconds 51. For reference, the half-life of azobenzene is 1.4
days in benzene solution at $35^{\circ}$C 52. Here we use the screening
results to identify redshifted derivatives with high or low barriers.

Figure 6 shows combinations of $\lambda_{\mathrm{cis}}$,
$\lambda_{\mathrm{trans}}$, and $\Delta G^{\mathrm{eff}}$, where $\lambda_{i}$
is the absorption wavelength of isomer $i$. We found that trans is the more
stable isomer for 99.6% of all species. The stable isomer is usually the one
activated by light, and so $\lambda_{\mathrm{trans}}$ is usually the quantity
of interest. Panels (a) and (c) show that redshifting to the near-IR is quite
difficult. Out of 19,000 compounds, only five have
$\lambda_{\mathrm{trans}}>600$ nm. Of these, only one has a barrier greater
than that of azobenzene. There are 1,641 species with
$\lambda_{\mathrm{trans}}>500$ nm (8.7%), including 475 with a barrier greater
than azobenzene (2.5%). However, the majority of predictions over 500 nm are
significant overestimates. As discussed below, most are actually closer to 470
nm.

Two redshifted species are shown in panels (d) and (e). Compound 7 has a high
barrier and compound 8 has a low barrier (top and bottom, respectively). We
confirmed these predictions using single-point SF-TDDFT calculations. The
quantum chemical activation free energies are shown below the model
predictions, and the two agree quite well. We used the model results for the
quasi-harmonic and conformational contributions to the enthalpy and entropy,
since a numerical Hessian with SF-TDDFT would be prohibitively expensive.

While the barriers agree quite well with SF-TDDFT, the trans absorption
wavelengths are overestimated. For example, the predicted and true wavelengths
are 567 nm and 467 nm for compound 7, and 574 and 533 nm for compound 8.
Moreover, this latter result is somewhat suspect, since the first excited
state with SF-TDDFT has square spin $\langle S^{2}\rangle=1.2$, indicating
high spin contamination. Indeed, a restricted TDDFT calculation with the
$\omega$B97X-D3 functional 53 and the def2-SVP basis 54 yielded
$\lambda_{\mathrm{trans}}=470$ nm. These results are common: after performing
calculations for the 40 species with the highest trans absorption wavelengths,
we found that all had either true values around 470 nm, or had significant
spin contamination leading to untrustworthy results.

We note that while the $S_{1}$ spin contamination was severe for some species
with ultra-high absorption wavelengths, it was otherwise low in general.
Indeed, the average square spin of the $S_{1}$ state was 0.37 in the training
set. This is higher than the mean value of 0.16 for the $S_{0}$ state, but
still reasonable. The maximum $\langle S^{2}\rangle$ allowed in the training
set was 1.5 for the $S_{1}$ state. To avoid any $S_{1}$ spin contamination,
one could always fine-tune the model with a small dataset of excitation
energies from spin-complete TDDFT and $\omega$B97X-D3. Multi-reference effects
would likely be small for equilibrium structures, and only a few new
calculations would be needed for fine-tuning 39, 6.

These results lead to several conclusions. First, redshifting is quite
difficult. This is shown by the fact that most model predictions above 500 nm
are actually error outliers. The associated wavelengths are still much higher
than the base compound, but quite far from the predicted values. Second, we
find that spin contamination is severe for many of the high-$\lambda$
predictions. Third, on a positive note, the model is able to identify
redshifted species with either high or low barriers, even though the redshift
is overestimated. This is encouraging for virtual screening in
photopharmacology.

The quantum chemistry absorption wavelengths also come with several sources of
uncertainty. They include errors in SF-TDDFT, implicit treatment of the
solvent, and use of static structures instead of thermally sampled geometries
6. Moreover, since the experimental absorption width is usually quite large,
compounds can often absorb at wavelengths 100-200 nm higher than their peak
55.

#### Graph-property relationships

Here we analyze the relationship between substitution and chemical properties.
Observing and explaining general trends will enable more focused candidate
screening in the future. Previous papers have also explored these
relationships computationally 17, 18, 20, 56, 12; we build on their
conclusions and extend them to other substitution patterns and groups.

Figure 7 shows barrier heights and absorption wavelengths by substitution
pattern. Panel (a) shows that motifs A, F, and G have bimodal $\Delta
G^{\mathrm{eff}}$ distributions, with high barriers around the second mode.
This is explored more below. Patterns B, C, and E have below-average barriers
and elongated distributions, while pattern H has a tight distribution and the
lowest mean barrier.

Panels (b) and (c) show $\lambda_{\mathrm{trans}}$ and
$\lambda_{\mathrm{cis}}$, respectively. Intriguingly, we see that the
distributions of classes E and H are very elongated for
$\lambda_{\mathrm{trans}}$. The same is true to a lesser extent for B and C.
While patterns like A, D, and F are tightly concentrated between 400 and 425
nm, class E samples from 400 to 525 nm with almost equal probability. The
$\lambda_{\mathrm{cis}}$ distributions are much tighter.

To better understand these results, we next analyze the relationship between
substituent properties and molecule properties for pattern A. Figure 8(a)
shows that $\lambda_{\mathrm{trans}}$ is maximized when $R_{1}$ is a non-ring
donor. In particular, panel (b) shows that this occurs when both $R_{1}$ and
$R_{2}$ are strong donors. Each box in this panel shows the root-mean-square
of $\lambda_{\mathrm{trans}}$ for the given ($R_{1}$, $R_{2}$) pair, computed
as
$\mathrm{min}\\{\lambda_{\mathrm{trans}}\\}+\mathrm{mean}\\{(\lambda_{\mathrm{trans}}-\mathrm{min}\\{\lambda_{\mathrm{trans}}\\})^{2}\\}^{0.5}$.
This gives a mean that is weighted towards higher values, reflecting our
interest in maximizing $\lambda_{\mathrm{trans}}$ better than a simple mean.
The results are somewhat unexpected: a simple picture of donors raising the
HOMO and acceptors lowering the LUMO would predict a redshift for strong
donors or acceptors. Yet the actual effect is only noticeable for donors.

Panel (c) shows that strong $R_{1}$ acceptors lead to very high barriers, with
narrow distributions centered around 33 kcal/mol. This is consistent with Ref.
18, which reported long thermal lifetimes for ortho-fluoro substituted
azobenzenes. Donor-substituted azobenzenes have noticeably lower barriers.
However, the associated barriers are tightly concentrated near 30 kcal/mol,
which is higher than that of unsubstituted azobenzene. This is encouraging,
since it means that substitution with two donors can redshift the absorption
wavelength without decreasing the barrier.

Panel (d) reinforces that $\Delta G^{\mathrm{eff}}$ has the opposite trend of
$\lambda_{\mathrm{trans}}$ (in this plot we use the mean for each box). Again,
however, we see that even the lowest barriers in the upper right corner are
similar to that of azobenzene. Intriguingly, panel (e) shows the opposite
trend for $\Delta G^{\dagger}$. This may be related to the absence of an
inversion mechanism for the ISC approach. This result reinforces the need for
high-accuracy quantum chemistry to accurately compare $\Delta G^{\dagger}$ and
$\Delta G^{\mathrm{eff}}$.

Panel (f) shows the relationship between the activation free energy and the
reaction free energy $\Delta G_{\mathrm{rxn}}$, where $\Delta
G_{\mathrm{rxn}}=G_{\mathrm{cis}}-G_{\mathrm{trans}}$. The Bell-Evans-Polanyi
principle 57, 58 states that the reaction enthalpy and activation enthalpy are
linearly related for reactions in the same family. This relationship was also
tested for azobenzene derivatives in Ref. 17. We see a moderate negative
correlation between the two quantities, with Spearman $\rho$ equal to $-0.27$.
Hence $\Delta G^{\mathrm{eff}}$ can be increased by making the cis isomer more
stable. However, the modest correlation means that this is not the full story,
and that the TS and MECP energies must be explicitly considered.

## Discussion

The main source of error in our approach is the underlying quantum chemistry
calculations. As discussed in SI Sec. S12, expensive wavefunction methods give
lower barriers than SF-TDDFT, but still give different answers from each
other. On balance the most accurate method seems to be SF-EOM-CCSD(dT), but
its prohibitive $N^{7}$ scaling 28 makes it a poor candidate for transfer
learning. Future work should focus on accurate quantum chemical approaches
that do not need manual setup, such as selection of active spaces, and that
are affordable enough for transfer learning.

Another limitation is that we have not considered azobenzene protonation and
azo-hydrazone tautomerism. These effects can be facilitated by substituents
such as $\mathrm{NH_{2}}$ and $\mathrm{OH}$, and by solvation in a protic
solvent, weakening the N=N double bond and lowering the isomerization barrier
59, 60, 61. Protonation from the solvent is not accounted for in a PCM
description. Incorporating automated protonation tools 62 into our workflow
would be of interest in the future. Similarly, the protein environment for a
given target in photopharmacology can also affect the isomerization rate 63,
23. Incorporating these effects for a specific target, as in Ref. 23, is of
interest for future work.

From the perspective of property optimization, the biggest remaining challenge
is redshifting. While it is straightforward to reach $\lambda=470$ nm for
trans isomers, it appears very difficult to reach $\lambda=$ 550-600 nm.
Averaging the gap over thermally sampled geometries may increase the
wavelength and improve prediction accuracy 6. Including bulky groups in all
four ortho positions may also increase the wavelength 6, but at the potential
cost of synthetic accessibility. A more targeted approach to wavelength
optimization could be of interest in the future. For example, methods such as
Monte Carlo tree search 64 could likely improve over virtual screening of
combinatorial libraries.

## Conclusions

We have presented a fast and automated method for predicting the isomerization
barriers of azobenzene derivatives. The approach can compute the activation
free energy through TS theory or ISC theory. We have demonstrated the accuracy
of the underlying ML model with respect to SF-TDDFT, reproduced trends in the
experimental isomerization rate, and argued for rotation-based ISC as the
reaction mechanism. Our software is fast, accurate, and easily accessible to
the community, making it a valuable tool for computational design of
photoactive molecules. Future work will focus on more accurate quantum
chemistry methods and more targeted molecular generation.

## Supporting Information

Extended methods, details of molecule generation, details of optimization,
non-adiabatic transition state theory, activation free energies by mechanism,
thermal isomerization rates, details of training, filtering protocol in
screening, collection of experimental data, note on mechanisms, and quantum
chemistry benchmark.

## Acknowledgements

We thank Dr. Johannes Dietschreit for bringing triplet-mediated isomerization
to our attention. Harvard Cannon cluster, MIT Engaging cluster, and MIT
Lincoln Lab Supercloud cluster at MGHPCC are gratefully acknowledged for
computational resources and support. Financial support from DARPA (Award
HR00111920025) is acknowledged.

## Supporting Information

## Thermal half-lives of azobenzene derivatives: virtual screening based on
intersystem crossing modeled with a machine learning potential

##  S1 Extended methods

All models used the PaiNN architecture 35 and were implemented in PyTorch 65.
As in our previous work, we used five convolutions instead of the three used
originally, as this substantially improved model performance 6. We also
allowed the $k$ values in the radial basis functions to be updated during
training. Remaining hyperparameters can be found in Ref. 6, and an in-depth
explanation of their meaning can be found in Ref. 35. Further hyperparameter
optimization is likely possible 66.

After training the model on data without empirical dispersion, we added D3BJ
dispersion 67, 68 to the total energy prediction:

$\displaystyle E=E_{\mathrm{model}}+E_{\mathrm{D3}}$ (S1)
$\displaystyle\mathbf{F}=-\nabla E_{\mathrm{model}}-\nabla E_{\mathrm{D3}}$
(S2)

The forces were computed with automatic differentation 69. Dispersion
parameters were taken from the PhysNet repository 70, 71. No cutoff was used.
D3 dispersion is a function only of the atom types and positions, and could
therefore be added analytically without approximation. The model performance
in the main text used SF-D3BJ TDDFT as the ground truth.

The geometries for SF-TDDFT/C-PCM calculations were generated through active
learning (Fig. 3(a)). In each round of active learning, the TS workflow in
Fig. 3(b) was performed for 1,000 (SMILES, mechanism) pairs using one of the
trained models. The pairs were randomly sampled from a set of 73,000 reactions
with four mechanisms each. We selected 3,000 generated geometries for DFT
calculations, added the new data to the training set, and repeated. 50% of the
geometries were selected by uncertainty, 30% by energy, and 20% randomly. For
uncertainty-based selection, we chose the geometries with the highest variance
in forces predicted by the three models. For energy selection, we sampled a
geometry $i$ with probability $p_{i}\propto e^{E_{i}/k_{\mathrm{B}}T}$. The
energy was computed relative to other geometries of the same species in the
same simulation. The energy selection sampled near-TS geometries from relaxed
scans, and high-energy MTD structures from conformer generation.

Active learning was repeated 15 times, yielding 43,938 calculations in total.
The final model was trained on 42,938 geometries, with 500 used for validation
and 500 for testing.

##  S2 Molecule generation

Molecules were generated using the patterns in Fig. 7(d). These are common
substitution patterns in the literature; see, for example, the synthesized
species in Refs. 72, 42, 43, 17, 44, 45, 46, 18, 20, 47. Some of the
functional groups were common azobenzene substituents aggregated in Ref. 6,
with new groups added from Refs. 73, 74, 75. Others were basic chemical
moieties, such as benzene, ethane, $\mathrm{\mathrm{R}-CF_{3}}$, and so on.

We aimed to avoid steric clashes in tetra-ortho substitution, since this can
make experimental synthesis more difficult. For this reason we separated
substituents by size. Those with four total atoms or fewer were labeled small,
and all others were labeled large. For all tetra-ortho substitutions, we only
allowed one side of the N=N bond to be substituted with large groups. The
other side could only be substituted with small groups. Substitution with
small groups at all positions was also allowed.

##  S3 Optimizations

###  S3.1 Relaxed scans

To perform the relaxed scans, we first identified the CNNC atoms using a
substructure search of azobenzene in RDKit 76. We then adjusted their angles
and/or dihedral angle by a constant amount in each step of the scan. The
angles and dihedrals were adjusted using the atomic simulation environment
(ASE) 77. For the inversion mechanism we set the final CNN or NNC angle to
179.5∘. We used 179.5∘ instead of 180∘ because the Cartesian derivative of the
latter is undefined, and hence the constraining forces described below would
also be undefined. For rotation we set the final dihedral angle to $\pm
90^{\circ}$. We also set the final CNN and NNC angles to 122∘. Without this
constraint, we found that some of the scans collapsed to inversion TSs instead
of rotation TSs. 122∘ is the angle of the rotational TS when optimizing
azobenzene with SF-TDDFT 24.

20 steps were used in each scan. After the angles and/or dihedrals were
adjusted, the geometry was optimized with forces given by
$\mathbf{F}+\mathbf{F}_{\mathrm{restrain}}$. The restraining forces are the
negative gradients of the following restraining energies:

$\displaystyle E_{\alpha,n}=k_{\alpha}(\alpha-\alpha_{n})^{2}$ (S3)
$\displaystyle E_{\omega,n}=-k_{\omega}\hskip 1.9919pt\mathrm{cos}\hskip
1.9919pt(\omega-\omega_{n}).$ (S4)

Here $\alpha$ denotes an angle, $\omega$ denotes a dihedral, $n$ denotes the
$n^{\mathrm{th}}$ step of the optimization, and $\alpha_{n},\ \omega_{n}$ are
target values at the $n^{\mathrm{th}}$ step. $E_{\alpha}$ was used for the
inversion mechanism, while both $E_{\alpha}$ and $E_{\omega}$ were used for
the rotation mechanism.

$k_{\alpha}$ and $k_{\omega}$ were each set to 1 Ha. The optimization was
performed with the BFGS algorithm 78, 79, 80, 81 in ASE. A rather loose
convergence threshold of $f_{\mathrm{max}}=0.05$ eV/Å was used. This was
because of the conformer generation performed after the scans. As described in
Sec. S3.2, each stage of conformer generation was restarted if a new conformer
had a lower energy than the seed conformer. We used a tolerance of
$f_{\mathrm{max}}=0.05$ eV/Å to minimize the optimization time spent on high-
energy conformers. Tight convergence criteria were only used once no lower
energy conformers were found. If the seed conformer were optimized with
tighter thresholds than the generated ones, it could have a lower energy
simply due to its stricter thresholds. Hence it needed to be optimized with
the same thresholds.

Angles and dihedrals were set with ASE at each step. ASE adjusts only these
internal coordinates, without rotating the rest of the molecule as a solid
body. RDKit uses a solid body rotation, and in principle should then require
fewer optimization steps at each stage of the scan. However, we found that
this led to unstable scans, and so we used ASE instead.

For all optimizations and dynamic simulations we updated the neighbor list
every ten steps. To account for neighbors coming into the cutoff radius
between updates, we computed the neighbors using a cutoff of
$r_{\mathrm{cut}}+r_{\mathrm{skin}}$, where $r_{\mathrm{cut}}=5.0$ Å is the
model cutoff radius and $r_{\mathrm{skin}}=2.0$ Å is the cutoff skin.
Distances between atoms and neighbors were computed at each step, and only
those within $r_{\mathrm{cut}}$ of each other were used in the model.

###  S3.2 Conformer generation

Figure S1: Outline of the conformer generation workflow used in this work. (a) Overview of the approach. (b) Details of the trimming procedure.  Parameter | Meaning | Value
---|---|---
$t_{\mathrm{new}}$ |  time between adding reference structures to $V_{\mathrm{bias}}$ | 1 ps
$\kappa$ |  turn-on parameter for adding new MTD Gaussians 82 | 0.03
$k^{\prime}$ |  $k=Nk^{\prime}$ is the biasing strength in Eq. (S5) 30 |  1.5 m$E_{\mathrm{h}}$
$\alpha$ |  Gaussian width in Eq. (S5) |  0.5 $\mathrm{Bohr}^{-2}$
$\Delta t$ | time step | 2 fs
$m_{\mathrm{H}}$ | hydrogen mass used in dynamics | 4 amu
thermostat | method for enforcing average temperature during dynamics | Nosé-Hoover
$Q$ | Nosé-Hoover effective mass | $(3N-6)\cdot\tau^{2}k_{\mathrm{B}}T$
$\tau$ | Nosé-Hoover relaxation time | 100 fs
$T$ | simulation temperature | $298.15\ \mathrm{K}$
$\Delta_{0}$ |  see Fig. S1(a) | 0.2 kcal/mol
$f_{\mathrm{max}}^{\mathrm{loose}}$ | maximum force component in loose opt | 0.2 eV/Å
$f_{\mathrm{max}}^{\mathrm{tight}}$ | maximum force component in tight opt | 0.05 eV/Å
$f_{\mathrm{max}}^{\mathrm{ultra}}$ | maximum force component in ultra-tight opt | 0.005 eV/Å
$E_{\mathrm{win}}^{\mathrm{loose}}$ | energy window for retaining conformers after loose opt |  15.0 $\mathrm{kcal}/\mathrm{mol}$
$E_{\mathrm{win}}^{\mathrm{tight}}$ | energy window for retaining conformers after tight opt |  8.0 $\mathrm{kcal}/\mathrm{mol}$
$E_{\mathrm{win}}^{\mathrm{ultra}}$ | energy window for retaining conformers after ultra-tight opt |  2.5 $\mathrm{kcal}/\mathrm{mol}$
$\Delta E_{\mathrm{thr}}$ |  two geometries cannot be duplicates or rotamers if their energies differ by more than $\Delta E_{\mathrm{thr}}$ | 0.15 kcal/mol
$\Delta R_{\mathrm{thr}}$ |  two geometries cannot be duplicates if their RMSD is more than $\Delta R_{\mathrm{thr}}$ | 0.175 Å
$\Delta B_{\mathrm{thr}}$ |  two geometries cannot be rotamers if their rotational constants differ by more than $\Delta B_{\mathrm{thr}}$ | 0.03 (3%)
Table S1: Parameters used in our NN conformer generation workflow. From top to
bottom, the sections contain MTD parameters, optimization parameters, and
cregen parameters. $N$ is the number of atoms in the system. The cregen
parameters are the defaults used in CENSO, a program that refines CREST
ensembles with DFT 83.

Conformer searches were performed for all reactants, products, and transition
states. We implemented a method based on CREST 30 using our NN. CREST combines
metadynamics, high-temperature molecular dynamics, and genetic structure
crossing 84 to sample conformational space. The sampled geometries are then
optimized with a multi-level filtering scheme. In this approach, optimizations
are performed with progressively tighter thresholds, and structures are
discarded if their energy is above a maximum value that is lowered in each
step. By default CREST uses the fast semi-empirical method GFN2-xTB 85 to
compute energies and forces.

The key component for phase space exploration is metadynamics (MTD). The
collective variables used for MTD are the root-mean-square displacements
(RMSDs) from previously visited structures. This drives the molecule away from
previously visited configurations and towards new regions of phase space. The
associated biasing potential is given by

$\displaystyle V_{\mathrm{bias}}=\sum_{i}^{n}k_{i}\hskip
1.9919pt\mathrm{exp}(-\alpha_{i}\Delta_{i}^{2}).$ (S5)

Here the sum is over $n$ previously visited structures, where a new structure
is added every 1 ps. $\Delta_{i}$ is the RMSD of the current structure with
respect to the $i^{\mathrm{th}}$ previous structure, $k_{i}$ is a pushing
strength parameter, and $\alpha_{i}$ is a width parameter. The forces at each
step are the sum of the actual forces and the negative gradient of
$V_{\mathrm{bias}}$. 14 MTD runs are performed in CREST, using different
combinations of $\alpha_{i}$ and $k_{i}$.

Initially we used CREST to perform conformer searches. However, given the
speed of relaxed scans and eigenvector following with our NN, we found that
CREST was by far the most time-consuming step in the workflow. For example,
with access to 24 Xeon-P8 CPU nodes and 1,152 cores in total, we ran 82 CREST
jobs at a time with 14 cores per job. The average job took 4.2 hours, which
meant a throughput of 470 geometries per day. Given four TS guesses per
species, plus the reactant and product geometries, this meant that only 78
species could be screened per day. Screening 25,000 species, as done in this
work, would have taken 11 months.

CREST has built-in options for conformer searches that are faster but less
extensive. These involve fewer MTD runs and tighter energy windows for each
level of optimization. However, given less extensive sampling and tighter
energy windows, it makes sense to use our NN instead xTB. For example, CREST
has a default energy window of 6.0 kcal/mol for the final ensemble, even
though conformers of energy greater than approximately $2.5$ kcal/mol are not
populated at room temperature. This is to compensate for errors in xTB, so
that low-lying conformers mistakenly assigned a high energy by xTB are not
discarded. These conformers can be subsequently re-ranked or re-optimized with
a higher level of theory. Therefore, it is more reasonable to use a tight
energy window for the NN, which is specifically trained to reproduce SF-TDDFT
for azobenzene derivatives. Similarly, the sampling should be more extensive
for xTB than for the NN. This is because errors in xTB mean that
conformational space must be thoroughly sampled up to an energy of 6.0
kcal/mol. The NN, by contrast, only needs to thoroughly sample the space up to
an energy of 2.5 kcal/mol. Alternatively, one could retain the extensive
sampling but use the ultra-fast GFN force-field 86. However, we found that the
force-field produced poor results for TSs.

We therefore implemented our own conformer search with the NN, and used
tighter energy windows and fewer dynamics simulations. The workflow is shown
schematically in Fig. S1. The approach closely follows that of CREST, but uses
only one MTD run per stage, and does not use MD or genetic structure crossing.
As in CREST, we determined the MTD time using a flexibility measure derived
from the chemical graph (the current formula in CREST differs from that of the
original publication; we used the up-to-date version in the source code 87,
commit 6bb6355). To partially compensate for the decreased total simulation
time, we used a minimum MTD time of 15 ps, instead of the 5 ps used in CREST.
The parameters used in our workflow are given in Table S1.

Several points require further discussion. First, we used the CREST cregen
tool to remove both duplicates and rotamers in the ensemble. Duplicates are
pairs of geometries with $\Delta E<\Delta E_{\mathrm{thr}}$,
$\mathrm{RMSD}<R_{\mathrm{thr}}$, and $B<B_{\mathrm{thr}}$, where $\Delta$ is
the difference between the two quantities, $E$ is the energy, $B$ is the
rotational constant, and $\mathrm{thr}$ denotes a threshold. Rotamers are
structures with $\Delta E<\Delta E_{\mathrm{thr}}$,
$\mathrm{RMSD}>R_{\mathrm{thr}}$, and $B<B_{\mathrm{thr}}$. We removed
rotamers because the difference in rotamer count per conformer should be
small, and because accounting for all rotamers is quite difficult. This is
especially true when MD and genetic structure crossing are not included, since
they are included in CREST primarily to find rotamers 30.

Second, we used the L-BFGS algorithm 88 in ASE to perform optimizations. The
convergence criterion is $\max_{i,\alpha}|f_{i,\alpha}|<f_{\mathrm{max}}$,
where $f_{i,\alpha}$ is a force component, $i\in\\{1,N\\}$ is the atom index,
$\alpha\in\\{x,y,z\\}$ is the Cartesian index, and $f_{\mathrm{max}}$ is a
threshold. For equilibrium geometries we used
$f_{\mathrm{max}}^{\mathrm{ultra}}=0.005$ eV/Å. For TS geometries we used
$f_{\mathrm{max}}^{\mathrm{ultra}}=0.01$ eV/Å, since the average change in
energy between the two thresholds is only 0.1 kcal/mol, while the number of
extra steps can be significant. We used $f_{\mathrm{max}}=0.005$ eV/Å for
eigenvector following performed on the five lowest-energy conformers. Note
also that the algorithm is implemented in Cartesian coordinates. This makes it
less efficient than the internal optimizer in xTB, which uses internal
coordinates. Use of an internal coordinate optimizer may be of interest in the
future.

Third, unlike in CREST 82, we did not constrain bond lengths with SHAKE 89.
Typically SHAKE is used to allow longer time steps and therefore accelerate
the dynamics. Indeed, CREST uses a default time step of 5 fs, while the
maximum value for unconstrained dynamics is typically 0.5 fs. However, the
SHAKE implementation in ASE is rather slow, and therefore became a bottleneck
instead of reducing run times. One can also increase the time step by
artificially increasing the mass of hydrogen, since, as the lightest element,
its motion is usually the fastest in the system. We followed the default in
CREST and set the hydrogen mass to 4 amu. This allowed us to use a time step
of 2 fs. Interestingly, we found that CREST automatically reduced the time
step from 5 fs to 2 fs for TS conformer searches, since trial MTD runs with
larger time steps all failed. This happened even though SHAKE is used in
CREST. The same thing did not happen for conformer searches of equilibrium
geometries. Hence our time step was the same as the CREST time step for TSs,
even though we did not use SHAKE. Yet an efficient implementation of SHAKE is
still of interest, since fixing the bond lengths during MTD might decrease the
number of subsequent optimization steps needed.

Fourth, we fixed the CNNC atoms in each molecule when performing TS conformer
searches. We experimented with Hookean constraints on the CNN angles and CNNC
dihedrals, but this led to unstable dynamics. This was likely because the
Cartesian gradient of an angle is undefined when the angle is $180^{\circ}$.
We excluded the fixed atoms from the RMSD computation in MTD. For equilibrium
geometries we excluded the CNNC atoms from the RMSD, but did not fix the
atoms. Excluding these atoms ensured that cis did not isomerize to trans.

Figure S2: Energy differences between TSs generated with NN conformer searches
and with CREST.

To evaluate our approach, we performed conformer searches with CREST and with
our method on 100 relaxed scan TS geometries. We then optimized each CREST
geometry with the NN, removed duplicates, and retained the geometries with
energies under 2.5 kcal/mol. The top 5 conformers from each method were then
optimized with eigenvector following. We checked that the mechanism of each TS
matched the target mechanism of the relaxed scan (e.g. rotational TS for a
rotation relaxed scan). We removed one species that had relaxed to an
inversion TS with CREST after being seeded with a rotational TS. The
distribution of TS energy differences is shown in Fig. S2(a). The mean value
of $E_{\mathrm{TS}}^{\mathrm{NN,\ MTD}}-E_{\mathrm{TS}}^{\mathrm{CREST}}$ is
0.28 kcal/mol. The mean increases to 0.48 kcal/mol after removing outliers
with absolute energy differences over 6 kcal/mol. The agreement between the
two methods is quite good, despite using only one MTD simulation in the NN
approach instead of 14. Further, the error in the $\Delta E^{\dagger}$ is
expected to be smaller than the error in $E_{\mathrm{TS}}$. This is because
the conformer search can only overestimate the lowest energy, not
underestimate it, and so some error cancellation can be expected.

We also compared the completeness of the ensembles in each of the two methods.
To do so we computed the conformational free energy, defined in SI Sec. S3.6.
Figure S2(b) shows the difference in conformational free energy between the
two methods. Interestingly, we see that the NN generates a more complete
ensemble on average. The mean value of $G_{\mathrm{conf}}^{\mathrm{NN,\
MTD}}-G_{\mathrm{conf}}^{\mathrm{CREST}}$ is $-0.36$ kcal/mol. This trend
persists even when the lowest-energy conformers are quite close in energy
between the two methods. For example, restricting ourselves to species in
which the $|E_{\mathrm{TS}}^{\mathrm{NN,\
MTD}}-E_{\mathrm{TS}}^{\mathrm{CREST}}|\leq 0.5$ kcal/mol, we find that the
difference in $G_{\mathrm{conf}}$ is $-0.60$ kcal/mol. We conclude that the NN
ensembles are quite complete, and that the associated $G_{\mathrm{conf}}$ is
reliable.

The NN approach is faster than CREST, mainly because fewer MTD runs are
performed. NN energy and force calculations also tend to be faster than xTB
for large molecules. A precise comparison is difficult, because xTB is run on
CPUs and NNs are run on GPUs. This means that one cannot simply compare wall-
clock times using equivalent hardware. Indeed, to compare the number of
species that can be screened per day, one must also consider the number of
GPUs available to the user vs. the number of CPUs. Nevertheless, some rough
comparisons can still be made. For example, running 189 ps of MTD for a
molecule with 102 atoms took 24 hours and 48 minutes with xTB on one Xeon-P8
core. If the xTB calculations were parallelized over all 48 cores on the node,
then the run would take 31 minutes given perfect parallelization. Hence 1.9
MTD runs could be performed per node per hour. By contrast, the corresponding
NN MTD run took 4 hours with one NVIDIA A100 GPU, using a batch of 50 species.
This corresponds to 12.5 MTD runs per GPU per hour, or 25 per node per hour,
assuming 2 GPUs per node. In this case, NN MTD is effectively $25/1.9\approx
13$ times faster than xTB MTD.

Lastly, since the conformer search is the rate-limiting step in our workflow,
any additional speed-ups would increase our throughput. While MTD enables
exploration of conformational space, it is fundamentally limited by the fact
that all steps must be taken in series. Uncorrelated geometries can only be
generated after taking several hundred steps. Stochastic conformer generators
are attractive because they can rapidly produce many uncorrelated structures.
Examples include the ETKDG method in RDKit 90 and the commercial software
Omega 91, 92. To adapt these generators for TSs, one would have to generate
conformers of equilibrium geometries, and perform a relaxed scan from the
equilibrium structures to the TS. Even if TS geometries would be generated
directly, the structures would still have to be optimized. In either case,
optimization would incur significant cost. Indeed, structure optimizations
take the most time of any step in CREST and NN conformer generation.
Therefore, avoiding MTD would not reduce the computation time to that of a
stochastic conformer generator. Another approach is to train a generative
model to produce TS geometries 93. ML for equilibrium conformer generation has
rapidly progressed in recent years 94, 95, 96, 97, 98, 99, 100, 101, 102, and
many of the methods could also be applied to TS generation. These avenues
would certainly be of interest in future work.

###  S3.3 Eigenvector following

EVF was implemented with Baker’s rational function optimization method 103. In
the first step, the numerical Hessian was computed with finite differences,
using the ASE vibrations package. We used a step size of 0.005 Å, and checked
that the results were indistinguishable from the analytical Hessian computed
with PyTorch. We used finite differences rather than automatic differentiation
to minimize the consumed GPU memory. In subsequent steps the Hessian was
updated using Powell’s method 104. We used a convergence threshold of
$f_{\mathrm{max}}=0.005$ eV/Å.

###  S3.4 Intrinsic reaction coordinate

Parameter | Meaning | Value
---|---|---
init_displ_de | Target energy change for the first step | 0.25 kcal/mol
scale_displ_sd | Factor for scaling the first SD step | 0.15
adapt_scale_disp |  Modify scale_displ_sd when the step size becomes smaller or larger | True
sd_parabolic_fit | Do a parabolic fit for finding the optimal SD step length | True
interpolate_only | Only allow interpolation for parabolic fit, not extrapolation | True
do_sd_corr | Apply a correction to the first SD step | True
scale_displ_sd_corr | Factor for scaling the correction to the SD step | 0.33
sd_corr_parabolic_fit | Do a parabolic fit for finding the optimal SD step length | True
tol_max_g | Maximum gradient for convergence |  $2\times 10^{-3}$ Ha/Bohr
tol_rms_g | RMS gradient for convergence |  $5\times 10^{-4}$ Ha/Bohr
Table S2: Parameters used in IRC searches. SD stands for steepest descent.

We implemented an IRC algorithm based on the approach in Orca 105. This
approach is itself based on the method of Morokuma and coworkers 33. Details
can be found in the Orca 5.0.2 manual. The parameters that we used are given
in Table S2. Note that the force tolerances are rather large, since only a
loose optimization is needed to see if the IRC has the product and reactant on
either end.

###  S3.5 Singlet-triplet crossing searches

Singlet-triplet crossings were located with a two-step procedure. In the first
step, we performed IRC to locate the crossings one either side of the TS. In
the second step, we optimized each of the two crossings by minimizing the
energy while keeping the singlet-triplet gap close to zero. This yielded MECPs
on each side of the TS.

In the first step we performed the normal IRC algorithm on the $S_{0}$ state.
We also kept track of the $T_{1}$ energy, and stopped the optimization once
$E_{T_{1}}-E_{S_{0}}$ changed sign. We then used the final geometry as the
starting point for the MECP optimization.

The MECPs were optimized with the approach of Ref 106. In a regular
optimization, the objective function is the energy $E$, and so its negative
gradient is $\mathbf{F}$. In an MECP optimization, the objective function is
the energy subject to the constraint that the singlet-triplet gap is zero. We
therefore used the following effective forces in the optimization:

$\displaystyle\mathbf{F}^{\prime}=\mathbf{P}\bar{\mathbf{F}}+\frac{1}{\alpha}\hskip
1.9919pt\Delta E\hskip 1.9919pt\Delta\mathbf{F}.$ (S6)

The first term is the mean of the singlet and triplet forces,
$\bar{\mathbf{F}}$, with the gradient of the energy gap $\Delta E$ projected
out by $\mathbf{P}$. Its effect is to minimize the average singlet and triplet
energies without increasing their gap. The second term contains the force
difference $\Delta\mathbf{F}$. Since it is scaled by the energy difference
$\Delta E$, its effect is to minimize the singlet-triplet gap. $\alpha$ is a
constant with units of energy. The terms are given by

$\displaystyle\Delta E=E_{T_{1}}-E_{S_{0}}$ (S7)
$\displaystyle\Delta\mathbf{F}=\mathbf{F}_{T_{1}}-\mathbf{F}_{S_{0}}$ (S8)
$\displaystyle\bar{\mathbf{F}}=\frac{1}{2}\left(\mathbf{F}_{T_{1}}+\mathbf{F}_{S_{0}}\right)$
(S9)
$\displaystyle\mathbf{P}\bar{\mathbf{F}}=\bar{\mathbf{F}}-\frac{\Delta\mathbf{F}^{T}\bar{\mathbf{F}}}{\Delta\mathbf{F}^{T}\Delta\mathbf{F}}\hskip
1.9919pt\Delta\mathbf{F},$ (S10)

where a superscript $T$ denotes transposition.

We implemented this approach in ASE, using a calculator that produced the
forces of Eq. (S6). We then performed the optimization using the BFGS
algorithm in ASE. The maximum force tolerance for convergence was set to 0.005
eV/Å. The constant $\alpha$ was set to 12.55 kcal/mol following a related MECP
algorithm for conical intersections 107. We found that the optimized MECPs
were insensitive to the choice of $\alpha$. For example, nearly identical MECP
energies were obtained when using $\alpha=1.0$ kcal/mol.

###  S3.6 Free energy calculations

The free energy was computed for reactants, products, and TSs, using

$\displaystyle G=H-TS,$ (S11)

where $G$ is the Gibbs free energy, $H$ is the enthalpy, and $S$ is the
entropy. The enthalpy was computed as

$\displaystyle H=E_{\mathrm{el}}+E_{\mathrm{ZPE}}+E_{\mathrm{vib}}(T)+4\hskip
1.9919ptk_{\mathrm{B}}T+E_{\mathrm{conf}},$ (S12)

where the last term is the average conformer energy,

$\displaystyle E_{\mathrm{conf}}=\sum_{i}p_{i}E_{i},$ (S13)

and the sum is over all conformers. $p_{i}$ is the Boltzmann probability for
conformer $i$, given by $p_{i}\propto\mathrm{exp}(-E_{i}/k_{\mathrm{B}}T)$.
$E_{\mathrm{el}}$ is the electronic energy, $E_{\mathrm{ZPE}}$ is the zero-
point energy, and $E_{\mathrm{vib}}(T)$ is the thermal vibrational energy.
Standard expressions for $E_{\mathrm{ZPE}}$ and $E_{\mathrm{vib}}(T)$ can be
found in Ref. 17. Apart from $E_{\mathrm{conf}}$, all terms in Eq. (S12) were
computed using the ASE thermochemistry package.

The entropy was computed as

$\displaystyle S=S_{\mathrm{mRRHO}}+S_{\mathrm{conf}},$ (S14)

where $S_{\mathrm{mRRHO}}$ is the modified rigid rotor-harmonic oscillator
entropy (see below), and $S_{\mathrm{conf}}$ is the conformational entropy:

$\displaystyle S_{\mathrm{conf}}=-k_{\mathrm{B}}\sum_{i}p_{i}\hskip
1.9919pt\mathrm{ln}\hskip 1.9919ptp_{i}.$ (S15)

The sum is over all conformers of all four mechanisms when using TS theory,
but only over rotational conformers when using the ISC approach. Note that if
two mechanisms are the same by symmetry, then the conformational entropy is
increased by $kT\hskip 1.9919pt\mathrm{log}\hskip 1.9919pt2$, while $\Delta G$
is lowered by this amount. This is the same factor that would arise from using
the symmetry number $\sigma=2$ in the calculations 19.

The modified rigid rotor harmonic oscillator approximation 108 interpolates
between a free rotor at low frequencies and a harmonic oscillator at high
frequencies. This is more physically sound than a standard harmonic
approximation, and avoids the divergent entropy of low-frequency harmonic
modes. The rotor cutoff frequency was set to 50 $\mathrm{cm}^{-1}$, following
the default in xTB 85. All other parameters were unchanged from Ref. 108.

Harmonic modes and frequencies were computed in the normal way 109 by
diagonalizing the mass-weighted Hessian after projecting out rotation and
translation. For reactants and products we checked that there were no
imaginary frequencies. For TSs we checked that there was only one imaginary
frequency, with magnitude $\geq 200$ $\mathrm{cm}^{-1}$. Species with
reactants or products not satisfying these conditions were removed. Species
were kept if they had at least one TS geometry from each mechanism satisfying
these conditions.

For $G^{\mathrm{X}}$, $H^{\mathrm{X}}$, and $S^{\mathrm{X}}$ we used
thermostatistical quantities from the rotational TSs. In principle the
vibrational terms should be computed with the average Hessian of the singlet
and triplet states at the crossing 106. However, we found that this approach
led to a sensitive dependence on $\alpha$ in Eq. (S6), and so we used the
rotational TS results instead.

##  S4 Non-adiabatic transition state theory

The reaction rate from non-adiabatic transition state theory (NA-TST) is given
by 106

$\displaystyle k_{\mathrm{NA-
TST}}=\frac{Z_{\mathrm{X}}}{hZ_{\mathrm{R}}}\int_{0}^{\infty}d\varepsilon\hskip
1.9919ptP(\varepsilon)\ \mathrm{exp}(-\beta\varepsilon).$ (S16)

Here $\beta=1/(k_{\mathrm{B}}T)$, $Z_{\mathrm{X}}$ and $Z_{\mathrm{R}}$ are
the rovibrational partition functions at the crossing point and reactant
geometry, respectively, and $P(\varepsilon)$ is the probability of
transitioning between the two electronic states at energy $\varepsilon$.
$P(\varepsilon)$ can be computed within the Wentzel-Kramers-Brillouin (WKB)
approximation 110, 111. The result is 34

$\displaystyle k_{\mathrm{NA-TST}}=k_{\mathrm{ISC}}\hskip
1.9919pt\mathrm{exp}(-\beta\Delta G^{X}),$ (S17)

where the intersystem crossing rate is

$\displaystyle
k_{\mathrm{ISC}}=\frac{\pi^{3/2}\alpha}{2h\sqrt{\lambda/(k_{\mathrm{B}}T)}}\left[1+\frac{1}{2}\hskip
1.9919pt\mathrm{exp}\left(\frac{1}{12\hskip 1.9919pt\alpha^{2}\hskip
1.9919pt(k_{\mathrm{B}}T\lambda)^{3}}\right)\right],$ (S18)

and

$\displaystyle\alpha=\frac{4H_{\mathrm{SO}}^{3/2}}{\hbar}\left(\frac{\mu}{F_{g}|\Delta\mathbf{F}|}\right)^{1/2}$
(S19) $\displaystyle\lambda=\frac{|\Delta\mathbf{F}|}{2F_{g}H_{\mathrm{SO}}}$
(S20) $\displaystyle
F_{g}=\bigg{|}\sum_{j=1}^{N}\sum_{n=1}^{3}\left(\mathbf{F}_{S_{0}}\right)_{jn}\left(\mathbf{F}_{T_{1}}\right)_{jn}\bigg{|}^{1/2}.$
(S21)

Here $H_{\mathrm{SO}}$ is the spin-orbit coupling, $\mu$ is the reduced mass
of the reaction coordinate, $F_{g}$ is the geometric mean of the singlet and
triplet forces at the crossing, $N$ is the number of atoms, and
$\Delta\mathbf{F}$ is the force difference (Eq. (S8)). In NA-TST, the reaction
coordinate is the direction of $\Delta\mathbf{F}$. The reduced mass is then
given by

$\displaystyle\mu=\left(\frac{1}{|\Delta\mathbf{F}|^{2}}\sum_{j=1}^{N}\sum_{n=1}^{3}\Delta\mathbf{F}_{jn}^{2}m_{j}^{-1}\right)^{-1},$
(S22)

where $m_{j}$ is the mass of the $j^{\mathrm{th}}$ atom. Notice that Eq. (S18)
scales quadratically with $H_{\mathrm{SO}}$, as expected from Fermi’s golden
rule 112. In re-expressing the equation from Ref. 34 we have used the relation
$G_{\mathrm{rv}}=-k_{\mathrm{B}}T\hskip 1.9919pt\mathrm{log}Z_{\mathrm{rv}}$
for the rovibrational free energy, and written
$G^{\mathrm{X}}=E^{\mathrm{X}}+G_{\mathrm{rv}}^{\mathrm{X}}$. We have also
corrected two typos (Eq. (7) in Ref. 34 should not contain $k_{\mathrm{B}}T$,
and $F_{g}$ involves a square root rather than a square).

In this work we have approximated $H_{\mathrm{SO}}$ as constant among all
azobenzene derivatives. We tested this approximation by computing
$H_{\mathrm{SO}}$ at the rotational TS of several derivatives. We included a
derivative with chlorine, since, as a relatively heavy third row element, it
could increase the coupling. Couplings are available in Q-Chem for TDFFT, but
not for SF-TDDFT or CASPT2. They are also not available for CASPT2 in Orca. We
therefore used TDDFT with both the B3LYP and BHHLYP functionals, and found
that $H_{\mathrm{SO}}\approx 40\ \mathrm{cm}^{-1}$, with a relative range of
25% between the smallest and biggest values. The coupling is within a factor
of two of the (14, 12) CASPT2 result 7, and the small variation indicates that
$H_{\mathrm{SO}}\approx\mathrm{const}.$ is a good approximation. It also
suggests that the couplings from CASPT2 would be close to constant. In our
calculations we therefore set $H_{\mathrm{SO}}=20\ \mathrm{cm}^{-1}$,
following the (presumably) more accurate CASPT2 result.

It is informative to estimate the impact of this constant coupling
approximation. To do so, we note that the ISC rate scales quadratically with
the spin-orbit coupling. Thus a relative change of 25% leads to a 56% change
in the rate. Converting this to an effective activation entropy, and using
$t_{\mathrm{ISC}}$ = 4.2 ps as discussed below, gives a change of 0.26
kcal/mol to $\Delta G^{\mathrm{eff}}$. Thus the maximum error from this
approximation is only 0.26 kcal/mol, and is thus negligible compared to other
sources of error.

The ISC rate must be multiplied by two for the singlet $\to$ triplet
transition, and by three for the triplet $\to$ singlet transition 7. Using Eq.
(S18) together with $H_{\mathrm{SO}}=20\ \mathrm{cm}^{-1}$, we found that
$t_{\mathrm{ISC}}=5.7$ ps for $S\to T$ and 4.2 ps for $T\to S$. These are
within a factor of two of the results quoted in Ref. 7 using Fermi’s golden
rule, thus validating both their approach and ours. Note that we use the same
spin-orbit coupling constant as Ref. 7, and so the differences come solely
from the different rate formulas. The benefit of Eq. (S18) is that it does not
contain the Franck-Condon factors used in Fermi’s golden rule expressions 7,
113. Such terms are computationally expensive, and cannot be easily computed
for anharmonic modes such as the reaction coordinate. We found that
$t_{\mathrm{ISC}}$ was essentially constant among derivatives: the minimum
time was 4.8 ps, and the maximum time was 5.8 ps.

##  S5 Connection between non-adiabatic and Eyring transition state theory

The intersystem crossing rate in NA-TST can be connected to the experimental
$\Delta S^{\dagger}$. Experimentally one starts by assuming an Arrhenius
transition rate,

$\displaystyle k_{\mathrm{A}}=Ae^{-E_{a}/(k_{\mathrm{B}}T)}.$ (S23)

Here $E_{a}$ is the activation energy and $A$ is the Arrhenius prefactor. The
two parameters are obtained experimentally from a plot of $\mathrm{log}\hskip
1.9919ptk$ vs. $1/T$ 19. Differentiating $\mathrm{log}\hskip 1.9919ptk$ with
respect to $1/T$ in Eqs. (2) and (S23), equating the two results, and invoking
Eq. (2) yields 19

$\displaystyle\Delta
S^{\dagger}=k_{\mathrm{B}}\left(\mathrm{log}\left(A\frac{h}{k_{\mathrm{B}}T}\right)-1\right),$
(S24)

where the slight temperature dependence of $S$ and $H$ has been neglected. If
the process is actually mediated by ISC, then we can apply the same logic to
Eq. (S17), which yields

$\displaystyle E_{a}=\Delta
H^{\mathrm{X}}+\frac{k_{\mathrm{B}}T}{2}-\frac{3\eta
k_{\mathrm{B}}}{T^{2}}\frac{\mathrm{e}^{\eta/T^{3}}}{2+\mathrm{e}^{\eta/T^{3}}}$
(S25) $\displaystyle\eta=\frac{1}{12\alpha^{2}k_{\mathrm{B}}^{3}\lambda^{3}}.$
(S26)

Inserting Eqs. (S25) and (S23) into (S17) then gives

$\displaystyle A_{\mathrm{NA-TST}}=k_{\mathrm{ISC}}\hskip
1.9919pt\mathrm{exp}\left(\frac{\Delta
S^{\mathrm{X}}}{k_{\mathrm{B}}}+\frac{1}{2}\right)\hskip
1.9919pt\mathrm{exp}\left(\frac{-3\eta\hskip
1.9919pte^{\eta/T^{3}}}{(2+e^{\eta/T^{3}})T^{3}}\right).$ (S27)

Typically $\eta/T^{3}\ll 1$, and so the last exponential in Eq. (S27) is of
order one. Neglecting this term and inserting the result into Eq. (S24) gives
the effective activation entropy,

$\displaystyle\Delta S^{\mathrm{eff}}=k_{\mathrm{B}}\left(\hskip
1.9919pt\mathrm{log}\hskip
1.9919pt\frac{hk_{\mathrm{ISC}}}{k_{\mathrm{B}}T}-\frac{1}{2}\right)+\Delta
S^{\mathrm{X}}.$ (S28)

This is the entropy that would be inferred from experiment, given a triplet-
mediated process with intersystem crossing rate $k_{\mathrm{ISC}}$, and an
entropy difference $\Delta S^{\mathrm{X}}$ between the intersection geometry
and the reactant geometry.

##  S6 Activation free energies by mechanism

SI Fig. S3 compares the activation free energies of the different mechanisms.
The plots are for the 19,000 compounds screened in the main text. Panel (a)
compares triplet-mediated isomerization with standard $S_{0}$ isomerization.
$\Delta G^{\mathrm{eff}}$ is lower than $\Delta G^{\dagger}$ in 65% of cases.
This indicates that, at the SF-TDDFT level of theory, isomerization proceeds
through $T_{1}$ for roughly 2/3 of the molecules. However, this result should
be interpreted with caution. As shown in Table S3, SF-TDDFT overestimates the
energy at both the MECP and the TS relative to the highly accurate SF-EOM-
CCSD(dT). The overestimate is 5.4 kcal/mol at the MECP and 3.9 kcal/mol at the
TS. This means that $\Delta G^{\mathrm{eff}}-\Delta G^{\dagger}$ is
overestimated by 1.5 kcal/mol. We therefore expect that, for the 19,000
compounds screened, $\Delta G^{\mathrm{eff}}$ should be lower than $\Delta
G^{\dagger}$ in more than 65% of cases. Indeed, if we assume that these
overestimates are constant among species, and thus subtract 1.5 kcal/mol from
all $\Delta G^{\mathrm{eff}}$, we find that the triplet mechanism is preferred
in 88% of cases. This result should also be interpreted with caution, however,
since it neglects errors in $\Delta S^{\mathrm{eff}}$ and $\Delta
S^{\dagger}$, and assumes that the energy error is constant among all species.

Panel (b) compares the rotation and inversion mechanisms in TST. $\Delta
G^{\dagger}_{\mathrm{rot}}$ is lower than $\Delta G^{\dagger}_{\mathrm{inv}}$
in 55% of all cases. Interestingly, this means that there is no strong
preference for either mechanism at the SF-TDDFT level of theory. This should
be contrasted with the species in Fig. 4, for which rotation was preferred in
every case.

Figure S3: Reaction barriers for different mechanisms. (a) $\Delta
G^{\mathrm{eff}}$ from ISC vs. $\Delta G^{\dagger}$ from TST on the $S_{0}$
surface. (b) $\Delta G^{\dagger}$ from TST for rotation vs. inversion.

##  S7 Thermal isomerization rates

SI Fig. S4 compares predicted and experimental thermal lifetimes, rather than
$\Delta G$ values. It contains the same information as Fig. 4 in the main
text, but with $\Delta G$ converted to a rate using Eq. (2), with $T=298.15$
K. The lifetime is the reciprocal of the rate. Notice that we first converted
experimental rates to $\Delta G$, given the temperature of the experiment, and
then back to a rate using $T=298.15$ K. This ensures that the lifetimes are
all compared at the same temperature. The line of best fit was computed for
$\Delta G$, and then converted to a lifetime.

We see that the experimental lifetimes span five orders of magnitude. Species
1 has the shortest lifetime at 10 minutes, while species 6 has the longest
lifetime at 2.5 years. The systematic overestimation of the energies at the
MECP and the TS leads to an overestimate of the lifetimes. The predicted
lifetimes are respectively 4, 5, and 7 orders of magnitude too high for TST,
ISC, and TST with only inversion. However, since the errors are systematic,
the model still has a high Spearman rank correlation with experiment.
Converting from $\Delta G$ to lifetime does not affect the ranking, and so the
Spearman rank coefficients are the same as in Fig. 4. $R^{2}$ is negative
because of the exponential dependence of lifetime errors on $\Delta G$ errors.

Figure S4: Experimental vs. predicted thermal lifetimes. Dotted orange lines
are linear regression results from predicted to experimental $\Delta G$,
converted to a lifetime. $\rho$ denotes the Spearman rank correlation. $R^{2}$
is computed between the regression results and the experimental data. (a)
Prediction accuracy using TS theory. (b) Prediction accuracy using intersystem
crossing. (c) As in (a), but with only the inversion mechanism. (d) Selected
compounds highlighted in panels (a)-(c).

##  S8 Training

###  S8.1 Singlet models

Three different models were pre-trained on 680,736 gas-phase SF-TDDFT/6-31G*
calculations from Ref. 6. 95% of the data was used for training, 4% was used
for validation, and 1% for testing. Each model was initialized with a
different random seed and trained on different train/validation/test splits.
Each model was then fine-tuned with SF-TDDFT/6-31G* using a C-PCM model of
water 36, 37, 38. The final models were trained on 42,938 geometries, with 500
used for validation and 500 for testing. Different splits were used to fine-
tune each of the final models.

The singlet ground-state corresponds to an excitation from a triplet reference
state in SF-TDDFT. We identified two singlets as the two states with the
lowest $\langle S^{2}\rangle$ of the three-lowest energy excitations 6. The
lower energy singlet was then taken as the ground singlet. Geometries with
$\langle S^{2}\rangle>1$ in the ground singlet state were discarded. The spin
contamination in the training set was fairly low, with $\mathrm{mean}\langle
S^{2}\rangle=0.16$.

Approximately 6% of species in the training set had non-zero charge ($+1$,
$+2$, or $+4$). The model does not use charge in the input; however, all
charge states in this work can be inferred from the number of bonds for each
atom. For example, if nitrogen has a single bond to four different atoms, then
it must have a partial charge of $+1$ to have a full octet. Since bonds can be
inferred from atom type and distances 67, the model can learn the effect of
different charge states from positions and atomic numbers alone.

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 computed with SF-TDDFT in vacuum. The set of
vacuum reference energies was used for both vacuum and solvent training.

Configurations with 10-$\sigma$ energy outliers were removed prior to
training. Those with forces $\geq$450 kcal/mol/Å from the mean or energies
$\geq$900 kcal/mol from the mean were also removed. Geometries were not
removed on the basis of $X$-$\sigma$ force deviations, since many geometries
were close to TSs and hence had near-zero forces.

Models were trained with the Adam algorithm, using a batch size of 60 for the
gas-phase models and 20 for the fine-tuned models. We used an MAE loss for the
forces and energies:

$\displaystyle\mathcal{L}\ =\ $ $\displaystyle\frac{1}{M}\hskip
1.9919pt\sum_{i=1}^{M}\rho_{E}\hskip 1.9919pt|E_{i}-\hat{E}_{i}|\ +\
\frac{1}{3\sum_{i=1}^{M}N_{i}}\hskip
1.9919pt\sum_{i=1}^{M}\sum_{j=1}^{N_{i}}\sum_{n=1}^{3}\rho_{F}\hskip
1.9919pt|F_{j,n}-\hat{F}_{j,n}|.$ (S29)

Here $\rho_{E}=0.2$ is the energy weight, $\rho_{F}=1.0$ is the force weight,
$M$ is the number of geometries in a batch, $N_{i}$ is the number of atoms in
geometry $i$, $\hat{X}$ is a predicted quantity and $X$ is its true value. We
used an MAE loss instead of a mean-squared-error loss, since this led to
significantly better model performance.

The learning rate was initialized to $10^{-4}$ and reduced by a factor of two
if the validation loss had not improved in $X$ epochs. $X$ was set to 10 for
the gas-phase models and 50 for the fine-tuned models. Training was stopped
when the learning rate fell below $10^{-5}$. 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 approximately 2 days for each of the
models.

###  S8.2 Triplet models

Triplet models were trained on the $S_{0}$/$T_{1}$ gap. They were pre-trained
using energies from Ref. 6. Gradients were not used since they had not been
computed. For fine-tuning we computed gradients of both the singlet and
triplet state for all geometries, and thus used both the gap and its gradient
for training. The same loss tradeoff was used as for the singlet model. The
learning rate was set $5\times 10^{-5}$ for pre-training and $10^{-4}$ for
fine-tuning; the former was chosen because of instabilities in pre-training.
To identify the triplet state, we first identified the two singlet states from
the three lowest-energy excitations using the method above. The triplet state
was then the remaining one. Data points with $\langle S^{2}\rangle<2.0$ or
$\langle S^{2}\rangle>2.4$ were discarded.

SF-TDDFT uses an $M_{\mathrm{S}}=1$ triplet reference state. The excited
$M_{\mathrm{S}}=0$ triplet and the reference $M_{\mathrm{S}}=1$ triplet should
have the same energy. We have found this to be true for spin-flip coupled-
cluster methods, but not for SF-TDDFT. This non-zero energy difference was
also observed in Ref. 114. The authors explained that the difference is due to
orbital relaxation, since the orbitals are relaxed for the reference state but
not for the TDDFT states.

It is not immediately clear which triplet energy should be used for ISC. We
tested results using the reference triplet, the excited triplet, and an
average of the two. We compared $E^{\mathrm{X}}$ to that of spin-flip coupled
cluster for azobenzene (SI Sec. S12), and found that the best agreement was
for the triplet average. The excited triplet overestimated $E^{\mathrm{X}}$ by
5 kcal/mol, and the reference triplet underestimated it by nearly the same
amount. However, we found that the correlation with experiment was slightly
better when we used the SF-TDDFT triplet. Further, this is more consistent
with the principle of treating singlet and triplet states on “equal footing”.
We therefore trained our models on the excited-state triplet energies. The
same network architecture and training parameters were used as for the singlet
models.

###  S8.3 Singlet excitation model

The singlet excitation model was trained directly on the $S_{0}$/$S_{1}$ gap.
It was first pre-trained using energies and gradients from Ref. 6. It was then
fine-tuned using energies only, since we did not compute $S_{1}$ gradients on
the new geometries.

##  S9 Filtering the screening results

We only kept the species that satisfied the following conditions:

  1. 1.

Endpoints are done. Both cis and trans must have conformer ensembles and a
minimum-energy geometry with no imaginary frequencies.

  2. 2.

Endpoints maintain cis/trans isomerism. The optimized cis and trans isomers
must actually be cis and trans, respectively. This is checked with RDKit.

  3. 3.

TS conformer generation is finished for each mechanism. A TS conformer
ensemble must be generated for each of the four mechanisms.

  4. 4.

At least one TS is converged for each mechanism. To be converged, the TS must
have only one imaginary frequency, and its magnitude must exceed 200
$\mathrm{cm}^{-1}$.

  5. 5.

Graph is unchanged. All of the TS and endpoint geometries must have the same
molecular graph as the input SMILES. This is checked using D3 coordination
numbers. Atom pairs with coordination numbers $\geq 0.95$ are assigned a bond.
The resulting bond list is checked against that of the original SMILES string.

  6. 6.

Vibrational frequencies are reliable. A frequency calculation is deemed
reliable if $|\Delta G^{\dagger}-\Delta E^{\dagger}|\leq 10$ kcal/mol.

  7. 7.

For the two singlet-triplet MECPs, one is closer to cis and the other to
trans. This is checked using the root-mean-square displacement of the CNNC
atoms after alignment.

  8. 8.

Each rotational TS connects to cis and trans through the IRC. This is true if
condition 7 is satisfied.

25,000 species were screened in total, and 18,877 remained after applying
these filters.

##  S10 Experimental data

Most papers reported either the thermal half-life $\tau_{1/2}$, or the thermal
lifetime $\tau$. These are related through $\tau=(e/2)\hskip
1.9919pt\tau_{1/2}$. In some cases the authors did not provide $\tau$ or
$\tau_{1/2}$, but did plot the absorption spectrum at different time
intervals. In these cases we inferred the thermal lifetime from the absorption
plots.

In all cases we computed the activation free energy from Eyring TST using
$\tau$:

$\displaystyle\Delta G^{\dagger}=k_{\mathrm{B}}T\hskip
1.9919pt\mathrm{log}\frac{k_{\mathrm{B}}T\tau}{h}.$ (S30)

Note that different works used different temperatures, and so $T$ depends on
the source.

##  S11 Note on mechanisms

###  S11.1 Concerted inversion

In our own calculations we have found that concerted inversion has two
imaginary frequencies, and therefore is not a true TS. This is because it is a
combination of two different inversions, each of which has one imaginary
frequency. It is therefore higher in energy than a single inversion, not a
true TS, and can be excluded from possible thermal mechanisms.

###  S11.2 Inversion-assisted rotation

A relaxed scan from $\alpha\approx 120^{\circ}$ to $\alpha\approx 180^{\circ}$
will either end in a pure inversion TS or an inversion-assisted rotation TS,
depending on the substituents. One can also be converted into the other
through a conformational search with fixed CNNC atoms. Therefore, we grouped
inversion and inversion-assisted rotation together.

##  S12 Quantum chemistry benchmark

Method | Basis | Restricted? |  | Treatment of
---
Static Correlation
| Treatment of
---
Dynamic Correlation
$\Delta E^{\dagger}$ | $\Delta E^{\mathrm{X}}$
DFT | cc-pVDZ | Yes | Poor | Strong |  35.6 | n/a
DLPNO-UCCSD(T) | F12-cc-pV(D, T)Z extrap. | No | Poor | Very strong | 42.1 | n/a
SF-TDDFT | 6-31G* | No | Moderate | Strong | 31.5 | 27.8
SF-TDDFT (D3BJ) | 6-31G* | No | Moderate | Strong | 32.6 | 28.9
CASPT2 (14, 12) | cc-pVDZ | Yes | Strong | Moderate | 25.4 | 18.7
SF-EOM-CCSD(dT) | 6-31G* | Yes | Moderate | Very strong | 28.7 |  23.5
Table S3: Summary of the methods and results from our benchmark. $\Delta
E^{\dagger}$ is the energy difference between the TS and the reactant, and
$\Delta E^{\mathrm{X}}$ is the energy difference between the singlet-triplet
crossing and the reactant. They are shown schematically in Fig. 2 in the main
text. DFT uses the B3LYP functional with D3BJ dispersion, while SF-TDDFT uses
the BHHLYP functional, both with and without D3BJ dispersion. Energy
differences are given in kcal/mol.

Here we compare the torsional energy profiles of unsubstituted azobenzene
predicted by different electronic structure methods. We analyze results from
DFT, DLPNO-UCCSD(T) (domain based local pair-natural orbitals 115, 116, 117,
118, 119, 120, 121, 122, 123, 124 for unrestricted coupled cluster with
single, double, and perturbative triple excitations), SF-TDDFT, CASPT2
(complete active space with second-order perturbation theory), and SF-EOM-
CCSD(dT) (spin-flip equation-of-motion coupled-cluster with single, double,
and perturbative triple excitations 28). The methods are summarized in Table
S3. For DFT we use the B3LYP functional 125 with D3 67 Becke-Johnson (BJ)
dispersion 68. For SF-TDDFT we use the BHHLYP functional both with and without
D3BJ dispersion.

For B3LYP-D3BJ we use the double-zeta correlation-consistent (cc) basis set
cc-pVDZ 126. This model chemistry predicts activation free energies of
azoarene isomerization that are in good agreement with experiment 11. For
DLPNO-UCCSD(T) we extrapolate to the complete basis-set limit using explicit
correlation (F12) 127 with double- and triple-zeta basis sets (cc-pV(D, T)Z).
Such high-accuracy basis sets can be used because of the computational
efficiency of DLPNO. For BHHLYP SF-TDDFT we use the 6-31G* basis 128, since
this is a medium-cost, medium-accuracy basis set that can be used to generate
large amounts of data for ML. Similar quality basis-sets are used for CASPT2
and SF-EOM-CCSD(dT) due to computational constraints. B3LYP-D3BJ, CASPT2, and
DLPNO-CCSD(T) calculations are performed with Orca 5.2 105. For DFT and
CCSD(T) calculations in Orca we use the resolution of identity 129, 130, 131,
132, 133, 134, 135 and chain-of-spheres 136 approximations (RIJCOSX). The
TightPNO setting is used for DLPNO calculations. All other calculations are
performed with Q-Chem 5.3.

Figure S5: Torsional energy profile of azobenzene computed with different
methods. Singlet and triplet energies are shown with full and dashed lines,
respectively. Singlet and triplet energies computed with the same method are
shown with the same color. (a) Results from all methods. For ease of
visualization, different subsets of methods are compared in panels (b)-(d).
(b) The multi-reference methods SF-EOM-CCSD(dT) and CASPT2 are compared with
the single-reference method CCSD(T). (c) The multi-reference methods SF-EOM-
CCSD(dT) and SF-TDDFT are compared with the single-reference method CCSD(T).
(d) The multi-reference method SF-EOM-CCSD(dT) is compared with the single-
reference methods B3LYP and CCSD(T).

We performed a relaxed scan from $\omega=0$ to $\omega=180^{\circ}$ in
$1^{\circ}$ steps using our trained NN potential. We constrained the dihedral
angle at each step using a force constant of 1 Ha. We then performed single-
point gas-phase calculations with the methods given above. Calculations were
performed on all geometries, except for SF-EOM-CCSD(dT), which was performed
on only nine geometries due to high computational cost. Those points are shown
as stars, and are connected by a smooth line using a quadratic interpolation.
Triplet calculations with DLPNO-CCSD(T) failed to converge and are hence not
reported.

Results are shown in Fig. S5. The plots contain several salient features.
First, the single-reference methods DFT and UCCSD(T) predict a cusp near
$\omega=90^{\circ}$, while the multi-reference methods produce a smooth
maximum of lower energy. This can be seen in Fig. S5(d), where DFT (red)
closely tracks UCCSD(T) (gray) until $\omega=90^{\circ}$. Both methods have a
cusp at $\omega=90^{\circ}$, but the energy from DFT is lower. Similar results
have been found for ethylene torsion 10.

As shown in Table S3, DFT gives an activation energy of 35.6 kcal/mol, while
UCCSD(T) gives 42.1 kcal/mol. The cusp in the DFT energy explains why previous
work could not optimize the rotational TS with DFT 19. Optimization methods
such as EVF 103 use a quadratic expansion of the PES to locate critical
points. This expansion fails in the vicinity of a cusp, since the Hessian is
undefined at this point.

Previous calculations with DFT found that inversion is the preferred mechanism
in gas phase 17. However, given that the DFT activation energy is 10 kcal/mol
higher than CASPT2, such conclusions are not to be trusted. Indeed, the
inversion activation energy was found to be 31 kcal/mol with CASPT2 (10, 8)
22, which is 6 kcal/mol higher than the CASPT2 rotation barrier here. Thus
according to CASPT2, isomerization would proceed by rotation even if it were
not mediated by ISC.

It is also interesting that the maximum value of the CCSD(T)
$T_{1}$-diagnostic 137 is 0.013, which occurs at $\omega=90^{\circ}$. The
minimum value is 0.01. A value above 0.02 indicates that the problem has
multi-reference character, and that CCSD(T) should not be trusted. The
diagnostic stays well below 0.02 despite the clear multi-reference nature of
the problem. We also note that the unrestricted CCSD(T) results collapsed to
the restricted results, with the square spin $\langle S^{2}\rangle=0$ for all
geometries.

A second important result is that all multi-reference methods predict smooth
TS maxima with singlet-triplet crossings on either side. Examples are CASPT2
and SF-EOM-CCSD(dT) in Fig. S5(b), and SF-EOM-CCSD(dT) and SF-TDDFT in Fig.
S5(c). However, the $\Delta E^{\dagger}$ and $\Delta E^{\mathrm{X}}$ predicted
by the different methods are in quantitative disagreement. The activation
energies predicted by CASPT2, SF-EOM-CCSD(dT), and SF-TDDFT (D3BJ) are 25.4,
28.7, and 32.6 kcal/mol, respectively. The singlet-triplet crossing energies
are 18.7, 23.5, and 28.9 kcal/mol respectively. Given the disagreement between
CASPT2 and SF-EOM, it is not immediately clear which method should be taken as
the “ground truth”.

On balance, however, it seems that SF-EOM-CCSD(dT) is likely more accurate.
One reason is that its predictions are in quantitative agreement with MR-
CISD+Q (multi-reference configuration interaction with singles and doubles,
plus a Davidson correction) for a standard rotational conical intersection
benchmark 138. Its treatment of static correlation should therefore be
sufficient for azobenzene, while its treatment of dynamic correlation is
better than that of CASPT2. Moreover, its predicted activation enthalpy is in
better agreement with experiment. The experimental activation enthalpy is 21.1
kcal/mol 21, which corresponds to $\Delta H^{\mathrm{X}}$ in the ISC picture.
Our NN calculations indicate that $\Delta H^{\mathrm{X}}\approx\Delta
E^{\mathrm{X}}-2.5$ kcal/mol for azobenzene. Therefore, $\Delta
H^{\mathrm{X}}\approx 21.0$ kcal/mol with SF-EOM-CCSD(dT), while $\Delta
H^{\mathrm{X}}\approx 16.2$ kcal/mol with CASPT2. The former agrees better
with experiment than the latter.

A third important result is that SF-TDDFT’s predictions are in qualitative
agreement with other multi-reference methods, but not quantitative agreement.
As in the other methods, SF-TDDFT predicts a smooth maximium with singlet-
triplet crossings on either side of the TS. However, its predicted barrier is
3.9 kcal/mol higher than that of SF-EOM, and 7.2 kcal/mol higher than that of
CASPT2. Its singlet-triplet crossing energy is 5.4 kcal/mol higher than SF-
EOM, and 10.2 kcal/mol higher than CASPT2. These errors make it difficult to
compare the ISC reaction rate with the Eyring reaction rate using SF-TDDFT.
Hence in the main text we simply assume that all azobenzene derivatives
isomerize through ISC.

Lastly, the plots show the importance of adding dispersion corrections to DFT.
Figure S5(a) includes results of SF-TDDFT with and without dispersion. These
are shown in green and cyan, respectively. When dispersion is included, all
single- and multi-reference methods are in good agreement away from
$\omega=90^{\circ}$. When dispersion is not included, SF-TDDFT overestimates
the cis energy, because the attractive force between the two benzene rings is
underestimated. Similar conclusions have been made in previous work 19, 11.
There is still some residual error in the cis energy, even with dispersion,
which does not occur with B3LYP-D3BJ. This is likely because the D3 parameters
for BHHLYP were chosen for ground-state DFT, not SF-TDDFT.

## References

  * Orgiu and Samorì 2014 Orgiu, E.; Samorì, P. $25^{\mathrm{th}}$ anniversary article: organic electronics marries photochromism: generation of multifunctional interfaces, materials, and devices. _Advanced Materials_ 2014, _26_ , 1827–1845
  * Kolpak and Grossman 2011 Kolpak, A. M.; Grossman, J. C. Azobenzene-functionalized carbon nanotubes as high-energy density solar thermal fuels. _Nano letters_ 2011, _11_ , 3156–3162
  * Broichhagen et al. 2015 Broichhagen, J.; Frank, J. A.; Trauner, D. A roadmap to success in photopharmacology. _Accounts of chemical research_ 2015, _48_ , 1947–1960
  * Lerch et al. 2016 Lerch, M. M.; Hansen, M. J.; van Dam, G. M.; Szymanski, W.; Feringa, B. L. Emerging targets in photopharmacology. _Angewandte Chemie International Edition_ 2016, _55_ , 10978–10999
  * Pyzer-Knapp et al. 2015 Pyzer-Knapp, E. O.; Suh, C.; Gómez-Bombarelli, R.; Aguilera-Iparraguirre, J.; Aspuru-Guzik, A. What is high-throughput virtual screening? A perspective from organic materials discovery. _Annual Review of Materials Research_ 2015, _45_ , 195–216
  * Axelrod et al. 2022 Axelrod, S.; Shakhnovich, E.; Gómez-Bombarelli, R. Excited state non-adiabatic dynamics of large photoswitchable molecules using a chemically transferable machine learning potential. _Nature Communications_ 2022, _13_ , 3440
  * Cembran et al. 2004 Cembran, A.; Bernardi, F.; Garavelli, M.; Gagliardi, L.; Orlandi, G. On the mechanism of the cis-trans isomerization in the lowest electronic states of azobenzene: $S_{0}$, $S_{1}$, and $T_{1}$. _Journal of the American Chemical Society_ 2004, _126_ , 3234–3243
  * Heindl and Wegner 2020 Heindl, A. H.; Wegner, H. A. Rational design of azothiophenes—substitution effects on the switching properties. _Chemistry–A European Journal_ 2020, _26_ , 13730–13737
  * Kuntze et al. 2022 Kuntze, K.; Viljakka, J.; Titov, E.; Ahmed, Z.; Kalenius, E.; Saalfrank, P.; Priimagi, A. Towards low-energy-light-driven bistable photoswitches: ortho-fluoroaminoazobenzenes. _Photochemical & Photobiological Sciences_ 2022, _21_ , 159–173
  * Shao et al. 2003 Shao, Y.; Head-Gordon, M.; Krylov, A. I. The spin–flip approach within time-dependent density functional theory: Theory and applications to diradicals. _The Journal of chemical physics_ 2003, _118_ , 4807–4818
  * Adrion et al. 2021 Adrion, D. M.; Kaliakin, D. S.; Neal, P.; Lopez, S. A. Benchmarking of density functionals for Z-azoarene half-lives via automated transition state search. _The Journal of Physical Chemistry A_ 2021, _125_ , 6474–6485
  * Mukadum et al. 2021 Mukadum, F.; Nguyen, Q.; Adrion, D. M.; Appleby, G.; Chen, R.; Dang, H.; Chang, R.; Garnett, R.; Lopez, S. A. Efficient discovery of visible light-activated azoarene photoswitches with long half-lives using active search. _Journal of Chemical Information and Modeling_ 2021, _61_ , 5524–5534
  * Griffiths et al. 2022 Griffiths, R.-R.; Greenfield, J. L.; Thawani, A. R.; Jamasb, A. R.; Moss, H. B.; Bourached, A.; Jones, P.; McCorkindale, W.; Aldrick, A. A.; Fuchter, M. J., et al. Data-driven discovery of molecular photoswitches with multioutput Gaussian processes. _Chemical Science_ 2022, _13_ , 13541–13551
  * Blaiszik et al. 2016 Blaiszik, B.; Chard, K.; Pruyne, J.; Ananthakrishnan, R.; Tuecke, S.; Foster, I. The materials data facility: data services to advance materials science research. _Jom_ 2016, _68_ , 2045–2052
  * Blaiszik et al. 2019 Blaiszik, B.; Ward, L.; Schwarting, M.; Gaff, J.; Chard, R.; Pike, D.; Chard, K.; Foster, I. A data ecosystem to support machine learning in materials science. _MRS Communications_ 2019, _9_ , 1125–1133
  * Bandara and Burdette 2012 Bandara, H. D.; Burdette, S. C. Photoisomerization in different classes of azobenzene. _Chemical Society Reviews_ 2012, _41_ , 1809–1825
  * Dokic et al. 2009 Dokic, J.; Gothe, M.; Wirth, J.; Peters, M. V.; Schwarz, J.; Hecht, S.; Saalfrank, P. Quantum chemical investigation of thermal cis-to-trans isomerization of azobenzene derivatives: substituent effects, solvent effects, and comparison to experimental data. _The Journal of Physical Chemistry A_ 2009, _113_ , 6763–6773
  * Knie et al. 2014 Knie, C.; Utecht, M.; Zhao, F.; Kulla, H.; Kovalenko, S.; Brouwer, A. M.; Saalfrank, P.; Hecht, S.; Bléger, D. ortho‐Fluoroazobenzenes: Visible light switches with very long‐lived Z isomers. _Chemistry–A European Journal_ 2014, _20_ , 16492–16501
  * Rietze et al. 2017 Rietze, C.; Titov, E.; Lindner, S.; Saalfrank, P. Thermal isomerization of azobenzenes: on the performance of Eyring transition state theory. _Journal of Physics: Condensed Matter_ 2017, _29_ , 314002
  * Schweighauser et al. 2015 Schweighauser, L.; Strauss, M. A.; Bellotto, S.; Wegner, H. A. Attraction or repulsion? London dispersion forces control azobenzene switches. _Angewandte Chemie International Edition_ 2015, _54_ , 13436–13439
  * Asano et al. 1981 Asano, T.; Okada, T.; Shinkai, S.; Shigematsu, K.; Kusano, Y.; Manabe, O. Temperature and pressure dependences of thermal cis-to-trans isomerization of azobenzenes which evidence an inversion mechanism. _Journal of the American Chemical Society_ 1981, _103_ , 5161–5165
  * Casellas et al. 2016 Casellas, J.; Bearpark, M. J.; Reguero, M. Excited-state decay in the photoisomerisation of azobenzene: a new balance between mechanisms. _ChemPhysChem_ 2016, _17_ , 3068–3079
  * Liang and Bakhtiiari 2022 Liang, R.; Bakhtiiari, A. Multiscale simulation unravels the light-regulated reversible inhibition of dihydrofolate reductase by phototrexate. _The Journal of Chemical Physics_ 2022, _156_ , 245102
  * Yue et al. 2018 Yue, L.; Liu, Y.; Zhu, C. Performance of TDDFT with and without spin-flip in trajectory surface hopping dynamics: cis$\rightleftharpoons$trans azobenzene photoisomerization. _Physical Chemistry Chemical Physics_ 2018, _20_ , 24123–24139
  * Lee et al. 2019 Lee, S.; Shostak, S.; Filatov, M.; Choi, C. H. Conical intersections in organic molecules: benchmarking mixed-reference spin–flip time-dependent DFT (MRSF-TD-DFT) vs spin–flip TD-DFT. _The Journal of Physical Chemistry A_ 2019, _123_ , 6455–6462
  * Becke 1993 Becke, A. D. A new mixing of Hartree–Fock and local density-functional theories. _The Journal of chemical physics_ 1993, _98_ , 1372–1377
  * Francl et al. 1982 Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. _The Journal of Chemical Physics_ 1982, _77_ , 3654–3665
  * Manohar and Krylov 2008 Manohar, P. U.; Krylov, A. I. A noniterative perturbative triples correction for the spin-flipping and spin-conserving equation-of-motion coupled-cluster methods with single and double substitutions. _The Journal of chemical physics_ 2008, _129_ , 194105
  * Axelrod et al. 2022 Axelrod, S.; Schwalbe-Koda, D.; Mohapatra, S.; Damewood, J.; Greenman, K. P.; Gómez-Bombarelli, R. Learning matter: Materials design with machine learning and atomistic simulations. _Accounts of Materials Research_ 2022, _3_ , 343–357
  * Pracht et al. 2020 Pracht, P.; Bohle, F.; Grimme, S. Automated exploration of the low-energy chemical space with fast quantum chemical methods. _Physical Chemistry Chemical Physics_ 2020, _22_ , 7169–7192
  * Axelrod and Gomez-Bombarelli 2022 Axelrod, S.; Gomez-Bombarelli, R. GEOM, energy-annotated molecular conformations for property prediction and molecular generation. _Scientific Data_ 2022, _9_ , 1–14
  * Axelrod and Gomez-Bombarelli 2020 Axelrod, S.; Gomez-Bombarelli, R. Molecular machine learning with conformer ensembles. _arXiv preprint arXiv:2012.08452_ 2020,
  * Ishida et al. 1977 Ishida, K.; Morokuma, K.; Komornicki, A. The intrinsic reaction coordinate. An ab initio calculation for HNC $\to$ HCN and $\mathrm{H}^{-}+\mathrm{CH_{4}}$ $\to$ $\mathrm{CH_{4}}^{+}\mathrm{H}^{-}$. _The Journal of Chemical Physics_ 1977, _66_ , 2153–2156
  * Liu et al. 2014 Liu, S.; Srinivasan, S.; Tao, J.; Grady, M. C.; Soroush, M.; Rappe, A. M. Modeling spin-forbidden monomer self-initiation reactions in spontaneous free-radical polymerization of acrylates and methacrylates. _The Journal of Physical Chemistry A_ 2014, _118_ , 9310–9318
  * Schütt et al. 2021 Schütt, K.; Unke, O.; Gastegger, M. Equivariant message passing for the prediction of tensorial properties and molecular spectra. International Conference on Machine Learning. 2021; pp 9377–9388
  * Truong and Stefanovich 1995 Truong, T. N.; Stefanovich, E. V. A new method for incorporating solvent effect into the classical, ab initio molecular orbital and density functional theory frameworks for arbitrary shape cavity. _Chemical Physics Letters_ 1995, _240_ , 253–260
  * Barone and Cossi 1998 Barone, V.; Cossi, M. Quantum calculation of molecular energies and energy gradients in solution by a conductor solvent model. _The Journal of Physical Chemistry A_ 1998, _102_ , 1995–2001
  * Cossi et al. 2003 Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Energies, structures, and electronic properties of molecules in solution with the C-PCM solvation model. _Journal of computational chemistry_ 2003, _24_ , 669–681
  * Ang et al. 2021 Ang, S. J.; Wang, W.; Schwalbe-Koda, D.; Axelrod, S.; Gómez-Bombarelli, R. Active learning accelerates ab initio molecular dynamics on reactive energy surfaces. _Chem_ 2021, _7_ , 738
  * Wang et al. 2020 Wang, W.; Yang, T.; Harris, W. H.; Gómez-Bombarelli, R. Active learning and neural network potentials accelerate molecular screening of ether-based solvate ionic liquids. _Chemical Communications_ 2020, _56_ , 8920
  * Gegiou et al. 1968 Gegiou, D.; Muszkat, K.; Fischer, E. Temperature dependence of photoisomerization. V. Effect of substituents on the photoisomerization of stilbenes and azobenzenes. _Journal of the American Chemical Society_ 1968, _90_ , 3907–3918
  * Nagamani et al. 2005 Nagamani, S. A.; Norikane, Y.; Tamaoki, N. Photoinduced hinge-like molecular motion: studies on xanthene-based cyclic azobenzene dimers. _The Journal of organic chemistry_ 2005, _70_ , 9304–9313
  * Sierocki et al. 2006 Sierocki, P.; Maas, H.; Dragut, P.; Richardt, G.; Vögtle, F.; De Cola, L.; Brouwer, F.; Zink, J. I. Photoisomerization of azobenzene derivatives in nanostructured silica. _The Journal of Physical Chemistry B_ 2006, _110_ , 24390–24398
  * Bandara et al. 2010 Bandara, H. D.; Friss, T. R.; Enriquez, M. M.; Isley, W.; Incarvito, C.; Frank, H. A.; Gascon, J.; Burdette, S. C. 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_ 2010, _75_ , 4817–4827
  * Bandara et al. 2011 Bandara, H. D.; Cawley, S.; Gascón, J. A.; Burdette, S. C. Short-circuiting azobenzene photoisomerization with electron-donating substituents and reactivating the photochemistry with chemical modification. _Chemistry – A European Journal_ 2011, _2011_ , 2916–2919
  * Beharry et al. 2011 Beharry, A. A.; Sadovski, O.; Woolley, G. A. Azobenzene photoswitching without ultraviolet light. _Journal of the American Chemical Society_ 2011, _133_ , 19684–19687
  * Moreno et al. 2016 Moreno, J.; Gerecke, M.; Grubert, L.; Kovalenko, S. A.; Hecht, S. Sensitized two-NIR-photon Z$\to$E isomerization of a visible-light-addressable bistable azobenzene derivative. _Angewandte Chemie International Edition_ 2016, _55_ , 1544–1547
  * Gómez-Bombarelli et al. 2016 Gómez-Bombarelli, R.; Aguilera-Iparraguirre, J.; Hirzel, T. D.; Duvenaud, D.; Maclaurin, D.; Blood-Forsythe, M. A.; Chae, H. S.; Einzinger, M.; Ha, D.-G.; Wu, T., et al. Design of efficient molecular organic light-emitting diodes by a high-throughput virtual screening and experimental approach. _Nature materials_ 2016, _15_ , 1120–1127
  * Wazzan et al. 2010 Wazzan, N. A.; Richardson, P. R.; Jones, A. C. Cis-trans isomerisation of azobenzenes studied by laser-coupled NMR spectroscopy and DFT calculations. _Photochemical & Photobiological Sciences_ 2010, _9_ , 968–974
  * Kahle and Zipoli 2022 Kahle, L.; Zipoli, F. Quality of uncertainty estimates from neural network potential ensembles. _Physical Review E_ 2022, _105_ , 015311
  * Mourot et al. 2011 Mourot, A.; Kienzler, M. A.; Banghart, M. R.; Fehrentz, T.; Huber, F. M.; Stein, M.; Kramer, R. H.; Trauner, D. Tuning photochromic ion channel blockers. _ACS chemical neuroscience_ 2011, _2_ , 536–543
  * Talaty and Fargo 1967 Talaty, E. R.; Fargo, J. C. Thermal cis–trans-isomerization of substituted azobenzenes: a correction of the literature. _Chemical Communications (London)_ 1967, 65–66
  * Chai and Head-Gordon 2008 Chai, J.-D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections. _Physical Chemistry Chemical Physics_ 2008, _10_ , 6615–6620
  * Weigend and Ahlrichs 2005 Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. _Physical Chemistry Chemical Physics_ 2005, _7_ , 3297–3305
  * Dong et al. 2017 Dong, M.; Babalhavaeji, A.; Collins, C. V.; Jarrah, K.; Sadovski, O.; Dai, Q.; Woolley, G. A. Near-infrared photoswitching of azobenzenes under physiological conditions. _Journal of the American Chemical Society_ 2017, _139_ , 13483–13486
  * Liu et al. 2018 Liu, X.-M.; Jin, X.-Y.; Zhang, Z.-X.; Wang, J.; Bai, F.-Q. Theoretical study on the reaction mechanism of the thermal cis–trans isomerization of fluorine-substituted azobenzene derivatives. _RSC advances_ 2018, _8_ , 11580–11588
  * Bell 1936 Bell, R. P. The theory of reactions involving proton transfers. _Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences_ 1936, _154_ , 414–429
  * Evans and Polanyi 1937 Evans, M.; Polanyi, M. On the introduction of thermodynamic variables into reaction kinetics. _Transactions of the Faraday Society_ 1937, _33_ , 448–452
  * Kojima et al. 2005 Kojima, M.; Nebashi, S.; Ogawa, K.; Kurita, N. Effect of solvent on cis-to-trans isomerization of 4-hydroxyazobenzene aggregated through intermolecular hydrogen bonds. _Journal of physical organic chemistry_ 2005, _18_ , 994–1000
  * Matazo et al. 2008 Matazo, D. R.; Ando, R. A.; Borin, A. C.; Santos, P. S. Azo- Hydrazone Tautomerism in Protonated Aminoazobenzenes: Resonance Raman Spectroscopy and Quantum-Chemical Calculations. _The Journal of Physical Chemistry A_ 2008, _112_ , 4437–4443
  * Rickhoff et al. 2022 Rickhoff, J.; Arndt, N. B.; Böckmann, M.; Doltsinis, N. L.; Ravoo, B. J.; Kortekaas, L. Reversible, Red-Shifted Photoisomerization in Protonated Azobenzenes. _The Journal of organic chemistry_ 2022, _87_ , 10605–10612
  * Pracht et al. 2017 Pracht, P.; Bauer, C. A.; Grimme, S. Automated and efficient quantum chemical determination and energetic ranking of molecular protonation sites. _Journal of Computational Chemistry_ 2017, _38_ , 2618–2631
  * Gaspari et al. 2017 Gaspari, R.; Prota, A. E.; Bargsten, K.; Cavalli, A.; Steinmetz, M. O. Structural basis of cis-and trans-combretastatin binding to tubulin. _Chem_ 2017, _2_ , 102–113
  * Dieb et al. 2019 Dieb, T. M.; Ju, S.; Shiomi, J.; Tsuda, K. Monte Carlo tree search for materials design and discovery. _MRS Communications_ 2019, _9_ , 532–536
  * Paszke et al. 2019 Paszke, A. et al. _Advances in Neural Information Processing Systems 32_ ; 2019; pp 8024–8035
  * Frey et al. 2022 Frey, N.; Soklaski, R.; Axelrod, S.; Samsi, S.; Gomez-Bombarelli, R.; Coley, C.; Gadepally, V. Neural Scaling of Deep Chemical Models. 2022,
  * Grimme et al. 2010 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _The Journal of chemical physics_ 2010, _132_ , 154104
  * Grimme et al. 2011 Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. _Journal of computational chemistry_ 2011, _32_ , 1456–1465
  * Paszke et al. 2017 Paszke, A.; Gross, S.; Chintala, S.; Chanan, G.; Yang, E.; DeVito, Z.; Lin, Z.; Desmaison, A.; Antiga, L.; Lerer, A. Automatic differentiation in pytorch. 2017,
  * Unke and Meuwly 2019 Unke, O. T.; Meuwly, M. PhysNet: A neural network for predicting energies, forces, dipole moments, and partial charges. _Journal of chemical theory and computation_ 2019, _15_ , 3678–3693
  * 71 PhysNet. https://github.com/MMunibas/PhysNet
  * Malkin and Fischer 1962 Malkin, S.; Fischer, E. Temperature dependence of photoisomerization. Part II. Quantum yields of cis-trans isomerizations in azo-compounds. _The Journal of Physical Chemistry_ 1962, _66_ , 2482–2486
  * Gutzeit et al. 2021 Gutzeit, V. A.; Acosta-Ruiz, A.; Munguba, H.; Häfner, S.; Landra-Willm, A.; Mathes, B.; Mony, J.; Yarotski, D.; Börjesson, K.; Liston, C., et al. A fine-tuned azobenzene for enhanced photopharmacology in vivo. _Cell Chemical Biology_ 2021, _28_ , 1648–1663
  * Lv et al. 2021 Lv, S.; Li, X.; Yang, L.; Ren, H.; Jiang, J. Computational design of photoswitchable anion receptors: Red-shifted and bistable di-ortho-fluoro di-ortho-chloro azobenzene derivatives. _Chemical Physics_ 2021, _548_ , 111246
  * Konrad et al. 2020 Konrad, D. B.; Savasci, G.; Allmendinger, L.; Trauner, D.; Ochsenfeld, C.; Ali, A. M. Computational design and synthesis of a deeply red-shifted and bistable azobenzene. _Journal of the American Chemical Society_ 2020, _142_ , 6538–6547
  * 76 RDKit: Open-source cheminformatics. http://www.rdkit.org
  * Larsen et al. 2017 Larsen, A. H. et al. The atomic simulation environment—a Python library for working with atoms. _Journal of Physics: Condensed Matter_ 2017, _29_ , 273002
  * Broyden 1970 Broyden, C. G. The convergence of a class of double-rank minimization algorithms 1. general considerations. _IMA Journal of Applied Mathematics_ 1970, _6_ , 76–90
  * Fletcher 1970 Fletcher, R. A new approach to variable metric algorithms. _The computer journal_ 1970, _13_ , 317–322
  * Goldfarb 1970 Goldfarb, D. A family of variable-metric methods derived by variational means. _Mathematics of computation_ 1970, _24_ , 23–26
  * Shanno 1970 Shanno, D. F. Conditioning of quasi-Newton methods for function minimization. _Mathematics of computation_ 1970, _24_ , 647–656
  * Grimme 2019 Grimme, S. Exploration of chemical compound, conformer, and reaction space with meta-dynamics simulations based on tight-binding quantum chemical calculations. _Journal of chemical theory and computation_ 2019, _15_ , 2847–2862
  * Grimme et al. 2021 Grimme, S.; Bohle, F.; Hansen, A.; Pracht, P.; Spicher, S.; Stahn, M. Efficient quantum chemical calculation of structure ensembles and free energies for nonrigid molecules. _The Journal of Physical Chemistry A_ 2021, _125_ , 4039–4054
  * Grimme et al. 2017 Grimme, S.; Bannwarth, C.; Dohm, S.; Hansen, A.; Pisarek, J.; Pracht, P.; Seibert, J.; Neese, F. Fully automated quantum-chemistry-based computation of spin–spin-coupled nuclear magnetic resonance spectra. _Angewandte Chemie International Edition_ 2017, _56_ , 14763–14769
  * Bannwarth et al. 2019 Bannwarth, C.; Ehlert, S.; Grimme, S. GFN2-xTB–An accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions. _Journal of chemical theory and computation_ 2019, _15_ , 1652–1671
  * Spicher and Grimme 2020 Spicher, S.; Grimme, S. Robust atomistic modeling of materials, organometallic, and biochemical systems. _Angewandte Chemie International Edition_ 2020, _59_ , 15665–15673
  * 87 Conformer-Rotamer Ensemble Sampling Tool. https://github.com/grimme-lab/crest
  * Liu and Nocedal 1989 Liu, D. C.; Nocedal, J. On the limited memory BFGS method for large scale optimization. _Mathematical programming_ 1989, _45_ , 503–528
  * Ryckaert et al. 1977 Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. _Journal of computational physics_ 1977, _23_ , 327–341
  * Riniker and Landrum 2015 Riniker, S.; Landrum, G. A. Better informed distance geometry: using what we know to improve conformation generation. _Journal of chemical information and modeling_ 2015, _55_ , 2562–2574
  * Hawkins et al. 2010 Hawkins, P. C.; Skillman, A. G.; Warren, G. L.; Ellingson, B. A.; Stahl, M. T. Conformer generation with OMEGA: algorithm and validation using high quality structures from the Protein Databank and Cambridge Structural Database. _Journal of chemical information and modeling_ 2010, _50_ , 572–584
  * Hawkins and Nicholls 2012 Hawkins, P. C.; Nicholls, A. Conformer generation with OMEGA: learning from the data set and the analysis of failures. _Journal of chemical information and modeling_ 2012, _52_ , 2919–2936
  * Pattanaik et al. 2020 Pattanaik, L.; Ingraham, J. B.; Grambow, C. A.; Green, W. H. Generating transition states of isomerization reactions with deep learning. _Physical Chemistry Chemical Physics_ 2020, _22_ , 23618–23626
  * Xu et al. 2020 Xu, M.; Luo, S.; Bengio, Y.; Peng, J.; Tang, J. Learning neural generative dynamics for molecular conformation generation. International Conference on Learning Representations. 2020
  * Shi et al. 2021 Shi, C.; Luo, S.; Xu, M.; Tang, J. Learning gradient fields for molecular conformation generation. International Conference on Machine Learning. 2021; pp 9558–9568
  * Ganea et al. 2021 Ganea, O.; Pattanaik, L.; Coley, C.; Barzilay, R.; Jensen, K.; Green, W.; Jaakkola, T. Geomol: Torsional geometric generation of molecular 3D conformer ensembles. _Advances in Neural Information Processing Systems_ 2021, _34_
  * Xu et al. 2021 Xu, M.; Wang, W.; Luo, S.; Shi, C.; Bengio, Y.; Gomez-Bombarelli, R.; Tang, J. An end-to-end framework for molecular conformation generation via bilevel programming. International Conference on Machine Learning. 2021; pp 11537–11547
  * Roney et al. 2021 Roney, J. P.; Maragakis, P.; Skopp, P.; Shaw, D. E. Generating realistic 3D molecules with an equivariant conditional likelihood model. 2021,
  * Luo et al. 2021 Luo, S.; Shi, C.; Xu, M.; Tang, J. Predicting molecular conformation via dynamic graph score matching. _Advances in Neural Information Processing Systems_ 2021, _34_
  * Jing et al. 2022 Jing, B.; Corso, G.; Barzilay, R.; Jaakkola, T. S. Torsional diffusion for molecular conformer generation. ICLR2022 Machine Learning for Drug Discovery. 2022
  * Xu et al. 2021 Xu, M.; Yu, L.; Song, Y.; Shi, C.; Ermon, S.; Tang, J. GeoDiff: A Geometric diffusion model for molecular conformation generation. International Conference on Learning Representations. 2021
  * Hoogeboom et al. 2022 Hoogeboom, E.; Satorras, V. G.; Vignac, C.; Welling, M. Equivariant diffusion for molecule generation in 3D. _arXiv preprint arXiv:2203.17003_ 2022,
  * Baker 1986 Baker, J. An algorithm for the location of transition states. _Journal of Computational Chemistry_ 1986, _7_ , 385–395
  * Powell 1971 Powell, M. J. Recent advances in unconstrained optimization. _Mathematical Programming_ 1971, _1_ , 26–57
  * Neese et al. 2020 Neese, F.; Wennmohs, F.; Becker, U.; Riplinger, C. The ORCA quantum chemistry program package. _The Journal of Chemical Physics_ 2020, _152_ , 224108
  * Lykhin et al. 2016 Lykhin, A. O.; Kaliakin, D. S.; dePolo, G. E.; Kuzubov, A. A.; Varganov, S. A. Nonadiabatic transition state theory: Application to intersystem crossings in the active sites of metal-sulfur proteins. _International Journal of Quantum Chemistry_ 2016, _116_ , 750–761
  * Levine et al. 2008 Levine, B. G.; Coe, J. D.; Martínez, T. J. Optimizing conical intersections without derivative coupling vectors: application to multistate multireference second-order perturbation theory (MS-CASPT2). _The Journal of Physical Chemistry B_ 2008, _112_ , 405–413
  * Grimme 2012 Grimme, S. Supramolecular binding thermodynamics by dispersion-corrected density functional theory. _Chemistry–A European Journal_ 2012, _18_ , 9955–9964
  * Ochterski 1999 Ochterski, J. W. Vibrational analysis in Gaussian. 1999; https://gaussian.com/vib/
  * Nikitin and Kearsley 1974 Nikitin, E. E.; Kearsley, M. J. _Theory of Elementary Atomic and Molecular Processes in Gases_ ; Clarendon Press: Oxford, 1974
  * Delos 1973 Delos, J. B. On the reactions of $\mathrm{N}_{2}$ with O. _The Journal of Chemical Physics_ 1973, _59_ , 2365–2369
  * Fermi 1995 Fermi, E. _Notes on quantum mechanics_ ; University of Chicago Press, 1995
  * Valiev et al. 2018 Valiev, R.; Cherepanov, V.; Baryshnikov, G. V.; Sundholm, D. First-principles method for calculating the rate constants of internal-conversion and intersystem-crossing transitions. _Physical Chemistry Chemical Physics_ 2018, _20_ , 6121–6133