# Inverse molecular design and parameter optimization with Hückel theory using
automatic differentiation

Rodrigo A. Vargas–Hernández  Chemical Physics Theory Group, Department of
Chemistry, University of Toronto, Toronto, Ontario, M5S 3H6, Canada.  Vector
Institute for Artificial Intelligence, 661 University Ave. Suite 710, Toronto,
Ontario M5G 1M1, Canada.  current affiliation: Department of Chemistry and
Chemical Biology, McMaster University, 1280 Main Street West, Hamilton,
Ontario, L8S 4M1 Canada<EMAIL_ADDRESS>Kjell Jorner  Chemical Physics
Theory Group, Department of Chemistry, University of Toronto, Toronto,
Ontario, M5S 3H6, Canada.  Department of Computer Science, University of
Toronto, 40 St. George St, Ontario M5S 2E4, Canada.  Department of Chemistry
and Chemical Engineering, Chalmers University of Technology, Kemigården 4,
SE-41258, Gothenburg, Sweden.  Robert Pollice  Chemical Physics Theory Group,
Department of Chemistry, University of Toronto, Toronto, Ontario, M5S 3H6,
Canada.  Department of Computer Science, University of Toronto, 40 St. George
St, Ontario M5S 2E4, Canada.  current affiliation: Stratingh Institute for
Chemistry, University of Groningen, Nijenborgh 4, Groningen, 9747 AG, The
Netherlands.  Alán Aspuru–Guzik  Chemical Physics Theory Group, Department of
Chemistry, University of Toronto, Toronto, Ontario, M5S 3H6, Canada.
Department of Computer Science, University of Toronto, 40 St. George St,
Ontario M5S 2E4, Canada.  Department of Chemical Engineering & Applied
Chemistry, 200 College St., University of Toronto, Ontario M5S 3E5, Canada.
Department of Materials Science & Engineering, 184 College St., University of
Toronto, Ontario M5S 3E4, Canada.  Vector Institute for Artificial
Intelligence, 661 University Ave. Suite 710, Toronto, Ontario M5G 1M1, Canada.
Lebovic Fellow, Canadian Institute for Advanced Research (CIFAR), 661
University Ave., Toronto, Ontario M5G 1M1, Canada.

###### Abstract

Semi-empirical quantum chemistry has recently seen a renaissance with
applications in high-throughput virtual screening and machine learning. The
simplest semi-empirical model still in widespread use in chemistry is Hückel’s
$\pi$-electron molecular orbital theory. In this work, we implemented a Hückel
program using differentiable programming with the JAX framework, based on
limited modifications of a pre-existing NumPy version. The auto-differentiable
Hückel code enabled efficient gradient-based optimization of model parameters
tuned for excitation energies and molecular polarizabilities, respectively,
based on as few as 100 data points from density functional theory simulations.
In particular, the facile computation of the polarizability, a second-order
derivative, via auto-differentiation shows the potential of differentiable
programming to bypass the need for numeric differentiation or derivation of
analytical expressions. Finally, we employ gradient-based optimization of atom
identity for inverse design of organic electronic materials with targeted
orbital energy gaps and polarizabilities. Optimized structures are obtained
after as little as 15 iterations, using standard gradient-based optimization
algorithms.

##  1 Introduction

Mathematical models that are both predictive and provide insight are a
cornerstone of the physical sciences. However, accurate models for complicated
processes often have no analytical solution and require large computational
resources to solve numerically. At the same time, they also tend to be hard to
interpret, as highlighted by Mulliken’s famous quote ”the more accurate the
calculations became, the more the concepts tended to vanish into thin air”
[1]. Approximate models with problem-specific parameters are therefore used in
practice, but finding optimal values for these parameters can be non-trivial.
Parameter optimization normally requires considerable amounts of reference
data and is done either manually or with algorithms that do not take advantage
of first or higher order derivatives as the corresponding analytical
expressions are often unavailable.

In chemistry, the Schrödinger equation is an archetype of such a mathematical
model that describes the interactions between nuclei and electrons in both
atoms and molecules. However, (near) exact solutions are too computationally
expensive for most molecules of interest. Quantum chemistry is an entire
research field dedicated to finding computationally efficient solutions to the
Schrödinger equation by introducing prudent approximations or reformulations
[2]. One approach that was extremely successful in the early days of quantum
chemistry is the use of so-called semiempirical (SE) approximations [3]. The
central idea is the use of problem-specific parameters to simplify the
mathematical form of the Schrödinger equation. One of the earliest SE models
was Hückel’s method to treat the $\pi$-electrons in organic molecules [4, 5,
6, 7]. Traditionally, the parameters in the Hückel method were derived
manually by human scientists with the aim to reproduce properties for well-
known reference molecules [8, 9], or they were derived from more accurate
calculations [10]. Over the years, the Hückel method has been used for
pedagogical purposes and for obtaining physical insight into problems in
organic chemistry [11] and photochemistry and photophysics [12, 13]. However,
it can also be used as a fast method for the prediction of molecular
properties [14], and for inverse design of molecules with desired target
properties [15, 16].

The recent upsurge in machine learning (ML), and specifically deep neural
networks, created a need for robust and efficient algorithms to co-optimize a
very large number of model parameters for various architectures. This problem
is now solved by automatic differentiation (AD), a technique to evaluate the
derivatives of mathematical expressions via the chain rule [17]. Importantly,
AD removes the need to determine analytic expressions for derivatives and
makes complicated mathematical models amenable to gradient-based optimization,
allowing them to be applied in the same way as general supervised machine
learning models. Regular machine learning approaches like deep neural networks
are meant to be very general mathematical models with a large number of
parameters. Through learning, they can adapt to essentially any problem
provided sufficient training data is available. In contrast, physics-based
mathematical models have expressions that are specific to a certain type of
problem to be solved and feature a much smaller number of parameters.
Implementing physical models such as quantum chemistry within AD frameworks
enables the use of default learning algorithms for parameter optimization with
a potentially much smaller training data requirement. Along these lines,
autodifferentiable versions of Hartree-Fock [18, 19], density functional
theory (DFT) [20, 21, 22, 19], excited state mean-field theory [23]. For semi-
empirical methods [24, 25, 26, 27], and other applications [28, 20, 29, 30,
30, 31, 32, 33, 34, 35, 36, 37], AD has been used to accelerate the
calculation of gradients physical methods and to blend with ML algorithms.

In this work, we developed an auto-differentiable implementation of the Hückel
method, by minimal adaptation of an initially developed NumPy [38] version
into the JAX [39] AD framework. We use this model to demonstrate the ease and
efficiency of parameter fitting based on computational reference data sets for
both excitation energies and molecular polarizabilities, a property calculated
via a second order derivative. Additionally, we demonstrate that our AD model
allows for gradient-based inverse design by regarding the atomic composition
of a molecular system as an adjustable parameter to find molecules with
targeted properties.[16] The corresponding code is made publicly available,
allowing it to be applied to a large variety of chemical problems. As the
Hückel calculations are extremely fast, our workflow allows for facile
development of property-specific models that can be readily used in molecular
generative models that require on the order of 105–106 property evaluations.

The paper is structured as follows: we first present a short introduction on
automatic differentiation and the Hückel model (Sections 2.1 and 2.2).
Following that, we execute inverse design of molecules as a fully
differentiable procedure (Section 3.1) and perform optimization of the Hückel
model parameters using modern gradient-based methods (Section 3.2).

##  2 Methods

###  2.1 Automatic differentiation

Gradients and high-order derivatives are at the core of physical simulations.
For physical models, common approaches to evaluate derivatives of any order
are closed-form solutions, symbolic differentiation, and numerical
differentiation, i.e., finite differences [40, 17]. For any function
represented as a computer program, AD [17] is an alternative way to compute
gradients and higher order derivatives. AD makes use of the chain rule for
differentiation to create a program that computes the gradients during
evaluation. There are two main modes in AD, forward and reverse mode. For
scalar functions, reverse mode is more efficient as differentiation requires a
single evaluation of the function to fully compute the Jacobian. An example of
reverse mode differentiation is the backpropagation algorithm that is used for
training neural networks. For more details about AD, we refer the reader to
Ref. [17].

The optimization of ML models is mostly done with methods that require the
gradient of the loss or error function (${\cal L}$) with respect to the model
parameters ($\boldsymbol{\theta}$), $\nabla_{\boldsymbol{\theta}}{\cal
L}(\boldsymbol{\theta})$. All contemporary ML libraries, e.g., Tensorflow
[41], PyTorch [42] and JAX [39], are built on top of an AD engine which
computes $\nabla_{\boldsymbol{\theta}}{\cal L}(\boldsymbol{\theta})$ for any
ML model. Given the robustness of AD libraries, differentiating physical
models [40] could be done similarly to modern ML algorithms.

###  2.2 Hückel model

The Hückel model, a well-known semi-empirical quantum chemistry model [4, 5,
6, 7], was first proposed to describe the interactions of $\pi$-electrons in
conjugated unsaturated hydrocarbons. In the Hückel model, this interaction is
restricted to electrons centered at nearest neighbour atoms. Generally, the
Hückel model is considered a tight-binding type Hamiltonian (Eq. 1) where the
on-site and hopping parameters are commonly denoted in the literature as
$\alpha_{\ell}$ and $\beta_{\ell,k}$, respectively. The matrix elements of the
Hückel Hamiltonian are given by

$H_{\ell,k}=\begin{cases}\alpha_{\ell},&\ell=k\\\ \beta_{\ell,k},&\ell\text{
and }k\text{ are adjacent }\\\ 0,&\text{ otherwise, }\end{cases}$ (1)

where the $\alpha_{\ell}$ parameters roughly represent the energy of an
electron in a 2p orbital, and the $\beta_{\ell,k}$ parameters describe the
energy of an electron in the bond $\ell-k$. Extensions of the Hückel
Hamiltonian are possible and can, for instance, incorporate distance-
dependence via $\beta_{\ell,k}=\beta^{0}_{\ell,k}\;g(\mathbf{R}_{\ell,k})$
(cf. Section 3.2). For more details, we refer the reader to standard quantum
chemistry textbooks [43, 44].

Notably, any molecular property computed with the Hückel method depends
directly on the $\alpha_{\ell}$ and $\beta_{\ell,k}$ parameters. Therefore, by
tuning their values, one can either construct a more accurate Hückel model for
a given molecule and property (cf. Section 3.2), or search for atomic
compositions that optimize target properties given a preset connectivity (cf.
Section 3.1). In the following sections, we demonstrate how AD can be used to
facilitate both these types of problems.

##  3 Results and Discussion

###  3.1 Inverse molecular design

Inverse molecular design can be carried out via gradient-based optimization
methods, as shown in Ref. [16]. The Hückel model can be extended to search for
the molecular structure with a desired property. Both the diagonal and off-
diagonal elements of the Hückel Hamiltonian matrix can be described by a
weighted average of atom types at each site,

$H_{\ell,k}=\begin{cases}\sum_{i}^{M}b_{\ell}^{i}\alpha^{i}_{\ell},&\ell=k\\\
\sum_{i}^{M}\sum_{j}^{M}b_{\ell}^{i}b_{k}^{j}\beta^{ij}_{\ell,k},&\ell\text{
and }k\text{ are adjacent }\\\ 0,&\text{ otherwise },\end{cases}$ (2)

where $b_{\ell}^{i}$ is the weight of the atom of type $i$ for site $\ell$.
For a meaningful description, the weights of each site must be normalized,
i.e., $\sum_{i}^{M}b_{\ell}^{i}=1$. $M$ is the total number of atom types
considered in the search. As a proof of concept, we consider eight different
molecular frameworks [16], which are displayed in Fig. 1. The $x$-symbol
indicates the sites with variable atom types to be optimized. We only
considered carbon ($C$), nitrogen ($N$) and phosphorus ($P$), i.e., $M=3$, as
these atom types each contribute one electron, assuming that the remaining
valences of carbon will be satisfied with a bond to an implicit hydrogen atom,
and can be incorporated interchangeably at all sites with two neighbors in the
$\pi$-framework (Fig. 1). Therefore we defined the following vector of atom
type weight parameters:
$\mathbf{b}_{\ell}=[b_{\ell}^{C},b_{\ell}^{N},b_{\ell}^{P}]$. For clarity,
$\mathbf{b}$ jointly describes the $\mathbf{b}_{\ell}$ parameters for all
search sites in a molecule, i.e., $\mathbf{b}=\\{\mathbf{b}_{\ell}\\}^{N}$.
For all results presented, the values of the $\alpha_{\ell}$ and
$\beta^{0}_{\ell,k}$ parameters were previously optimized with respect to the
desired property (cf. Section 3.2).

For the set of eight molecules considered (cf. Fig. 1), we search for the type
of atom $\mathbf{b}$ at each site that gives the lowest HOMO-LUMO gap (Eq. 3),
denoted as $\epsilon_{HL}$, and the maximum polarizability denoted as
$\langle\tilde{\alpha}\rangle$, (Eq, 4). $\epsilon_{HL}$ is defined as,

$\epsilon_{HL}=\epsilon_{LUMO}-\epsilon_{HOMO},$ (3)

where $\epsilon_{HOMO}$ and $\epsilon_{LUMO}$ are the eigenvalues of the
highest occupied molecular orbital (HOMO), and lowest unoccupied molecular
orbital (LUMO), respectively. The polarizability function is defined as

$\langle\tilde{\alpha}\rangle=\frac{1}{3}\left(\tilde{\alpha}_{xx}+\tilde{\alpha}_{yy}+\tilde{\alpha}_{zz}\right),$
(4)

where the $\tilde{\alpha}_{ij}$ elements are the polarizability components
defined as

$\tilde{\alpha}_{ij}=-\frac{\partial^{2}E}{\partial F_{i}\partial F_{j}}.$ (5)

The $F_{i}$ terms are the components of the electric field,
$\vec{\boldsymbol{F}}=[F_{x},F_{y},F_{z}]$, and $E$ is the electronic energy
of the system. The elements of the polarizability tensor are usually computed
using a finite-difference (FD) approach [16],

$\tilde{\alpha}_{ii}=\frac{2E(0)-\left[E(-F_{i})+E(+F_{i})\right]}{F_{i}^{2}},$
(6)

where the electronic energy is evaluated several times, typically three times
for each diagonal element, and four times for each cross term. Notably, if the
parameters $\mathbf{b}$ are to be optimized using a gradient-based method, the
Jacobians $\nabla_{\mathbf{b}}\epsilon_{HL}$ and
$\nabla_{\mathbf{b}}\langle\tilde{\alpha}\rangle$ will also be constructed
using an FD approach. However, this will increase the number of energy
evaluations needed, especially for $\langle\tilde{\alpha}\rangle$ as the
elements of $\nabla_{\mathbf{b}}\langle\tilde{\alpha}\rangle$ are third-order
derivatives:
$\frac{\partial}{\partial\mathbf{b}^{i}}\frac{\partial^{2}\tilde{\alpha}_{kk}}{\partial
F^{2}_{k}}$. Thus, for a single element of
$\nabla_{\mathbf{b}}\langle\tilde{\alpha}\rangle$, using FD will require 18
energy calculations, ${\cal O}(18\times\|\mathbf{b}\|)$, where
$\|\mathbf{b}\|$ is the total number of parameters in $\mathbf{b}$. For
$\epsilon_{HL}$, using FD, we only require ${\cal O}(2\times\|\mathbf{b}\|)$,
as $\nabla_{\mathbf{b}}\epsilon_{HL}$ is a first-order derivative.

In contrast, using modern AD frameworks, we can efficiently compute, with a
single energy calculation (i.e., one forward pass), the Jacobian of
$\epsilon_{HL}$ with respect to $\mathbf{b}$. For
$\langle\tilde{\alpha}\rangle$, the number of total energy evaluations depends
on the dimension of the external field to construct the diagonal elements of
the Hessian (Eq. 5), which we also computed using AD. The Jacobian of
$\langle\tilde{\alpha}\rangle$ with respect to $\mathbf{b}$, a third order
derivative, can be constructed from only three energy evaluations using AD
[17], a drastic reduction from the 18 required for FD.

After implementation of the Hückel model using the JAX ecosystem [39], we
could fully differentiate both observables, $\epsilon_{HL}$ and
$\langle\tilde{\alpha}\rangle$. Importantly, using JAX allowed us to convert
our existing Python-based Hückel code very easily by replacing calls to NumPy
with almost equivalent calls to the JAX.Numpy package. For optimization, we
used the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm via the JaxOpt
library [45]. Instead of using a constrained optimization scheme to satisfy
the site-normalization restriction for $\mathbf{b}$, we used the softmax
function,

$b^{i}=\frac{\exp(\underline{b}^{i})}{\sum_{i}^{M}\exp(\underline{b}^{i})},$
(7)

where $\underline{b}$ represents the unnormalized $b$ parameters.

Given the flexibility of the JAX ecosystem, we were able to test other
gradient optimization algorithms such as Adam [46] and canonical gradient-
descent, but found the BFGS to be most efficient as it required, on average,
fifteen or less total iterations to reach convergence (cf. Fig. 2). The Adam
and gradient-descent algorithms, with an exponential learning rate decay, each
needed more than thirty iterations to minimize $\epsilon_{HL}$ or
$\langle\tilde{\alpha}\rangle$. Notably, we initialized the values for all
$\underline{b}^{i}$ parameters by sampling a uniform distribution,
$\underline{b}^{i}\sim{\cal U}(-1,1)$. Instead of using literature Hückel
parameters, we used our optimized parameters for each target observable where
5,000 training molecules were used. More details are described in Section 3.2.
As seen in Figs. 3 and 4, we found that our random initialization of
$\mathbf{\underline{b}}$ allows us to sample a wide range of molecules with a
broad range of values for both objectives, $\epsilon_{HL}^{\mathrm{initial}}$
and $\langle\tilde{\alpha}\rangle_{\mathrm{initial}}$.

Because of the statistical description of the molecules by the $\mathbf{b}$
parameters, the optimal parameters ($\mathbf{b}^{*}$) found by optimization
are not one-hot vectors that correspond to only one atom type per site in the
molecule but rather a linear combination of multiple atom types. We define the
observable value for this unphysical molecule as $y_{\mathrm{virtual}}$ and
that for the real molecule as $y_{\mathrm{feasible}}$ (i.e., the most probable
atom is picked for each site to define the real molecule). An example of this
is displayed for framework 3 in Fig. 5 where we show the change in
$\mathbf{b}$ throughout the optimization for both objectives ($\epsilon_{HL}$
and $\langle\tilde{\alpha}\rangle$). As we can observe, the change in
$\epsilon_{HL}$ from the initial random molecule to the optimal one is close
to 1 eV. The optimizations in Fig. 5 were done using Adam only to properly
illustrate the change in $\mathbf{b}$ as the change between iterations is
smoother and more readily discernable. Importantly, this shows that the
generative model can shift the distribution of properties towards the target
property with little dependence on the random initialization of $\mathbf{b}$
(Figs. 3 and 4). We also observe that, for the majority of the optimized
molecules, $y_{virtual}$ and $y_{feasible}$ are linearly correlated,
indicating that the optimization converged essentially to feasible molecules.
The property distributions of $y_{\mathrm{virtual}}$ and
$y_{\mathrm{feasible}}$ are reasonably close, even in the few cases when the
correlations are poor.

Fig. 6 displays the molecules with the lowest $\epsilon_{HL}$ and maximum
$\langle\tilde{\alpha}\rangle$ from the ensemble of different optimizations.
First, we notice that there is a higher amount of phosphorus atoms in the
molecules when $\langle\tilde{\alpha}\rangle$ was the target property. This is
not unexpected as molecular polarizabilities, while not simply a sum of the
atomic polarizabilities, are strongly influenced by the atomic
polarizabilities of the constituent atoms [47, 48]. As phosphorus is a third-
row element in the same group with nitrogen and atomic polarizability
increases significantly when increasing the row number, its atomic
polarizability, both in free atoms [49] and in molecules [50], is
significantly larger than both nitrogen and carbon. Therefore, incorporating a
large number of phosphorus atoms is expected to be a viable strategy to
maximize the molecular polarizability in all of the molecular frameworks
considered. For the molecules with the lowest $\epsilon_{HL}$, we see
extensive incorporation of both nitrogen and phosphorous atoms. This can be
understood in terms of the effect of heteroatom substitution on
$\epsilon_{HL}$ within the HMO framework [51]. For alternant compounds such as
1–8, $\epsilon_{HL}$ is unaffected by changing $\alpha_{\ell}$. The main
effect comes from changing $\beta_{\ell,k}$, and is expected to be largest for
bonds that feature a bonding interaction in the HOMO and an anti-bonding
interaction in the LUMO. A lowered $\beta_{\ell,k}$ leads to decreased bonding
interactions in the HOMO and consequently, a higher HOMO energy. For the LUMO,
decreasing the antibonding interactions by a lowered $\beta_{\ell,k}$ leads to
a lowering of the energy. The net effect by raising the HOMO and lowering the
LUMO is a decrease in the $\epsilon_{HL}$ gap. We can therefore expect
optimization to favor atom pairs with a low $\beta_{\ell,k}$ for bonds that
feature a bonding interaction in the HOMO and an antibonding interaction in
the LUMO. Inspection of the molecular orbitals of the optimized frameworks
(cf. Figure S1) indeed reveals that these bonds are dominantly between two N
atoms, which feature the by far lowest $\beta_{\ell,k}$ at 0.159 (the next
lowest is for P–P at 0.539). Control optimizations with the original parameter
set by Van-Catledge [10] instead gives molecules with P–P at those bonds (cf.
Figure S2), consistent with the fact that $\beta_{P,P}$ = 0.63 is the lowest
for this parameter set.

For this proof of principle work, we picked two distinct molecular target
properties. However, based on the framework employed, a significant number of
alternative properties could also be predicted and, thus, used for inverse
design via gradient-based optimization. Additionally, any combined objective
that is derived from multiple target properties can equally be optimized for
via the same types of algorithms out of the box. This is particularly
interesting for properties where Hückel models are known to provide reasonable
prediction accuracies such as HOMO-LUMO gaps. The use of gradient-based
optimization algorithms enables fast convergence towards the closest local
optimum solution reducing the number of evaluations and leading to
significantly increased computation time. This is particularly important as
one of the main bottlenecks in current approaches to inverse molecular design
is the number of property evaluations needed to find an optimal structure
[52]. Going beyond single-objective optimization, one possible extension of
our presented approach is targeting multiple objectives via genuine gradient-
based multi-objective optimization, for example, both $\epsilon_{HL}$ and
$\langle\tilde{\alpha}\rangle$. The standard approach to perform multi-
objective optimization is via via property concatenation into a single
function to use standard single-objective algorithms, where algorithms like
Bayesian optimization are used [53, 54, 55, 56, 57]. However, gradient-based
multi-objective optimization algorithms [58] have been developed and they,
together with automatic differentiation, could be employed for both parameter
optimization and inverse molecular design in order to explore the
corresponding Pareto fronts in a systematic manner.

From a conceptual point of view, representing chemical structure subspaces in
a parameterized form can greatly facilitate inverse design [16] as it allows
the use of well-established approaches for parameter optimization to be used
for the design of molecules. This is particularly effective when used in
combination with AD due to its numerical stability and computational
efficiency compared to alternative means to compute gradients. Consequently,
this also makes the molecular size that can still be feasibly treated in such
an approach larger and thus, essentially, expands the chemical subspace the
generative model can explore. However, one of the main downsides of the
approach implemented in this work is the reliance on fixed molecular
frameworks, which is common for alchemical formulations [59] strongly limiting
the structural space considered in the optimization. Simple extensions would
be i) the combination of methods to change the molecular framework without
relying on gradients with the method presented here to modify the atom
identities within the respective framework, or ii) differentiable supermatrix
structure where atom vacancies are allowed [60]. Ideally, future extensions
should aim to find prudent ways allowing for framework modifications based on
gradients as this potentially can lead to a dramatic reduction in the number
of structure optimization steps and, thus, the number of property evaluations
necessary. The extended Hückel model is also compatible with the proposed
methodology, even with ML learned parameters [25], by considering a
description of the overlap integrals between different atoms types, similar to
Eq. 2.

Figure 1: The eight molecular frameworks considered for inverse molecular
design in this work [16]. The $x$-symbol represents atomic sites whose
identity is optimized.

Figure 2: Average learning curve for $\epsilon_{HL}$ (left panel) and
$\langle\tilde{\alpha}\rangle$ (right panel) for 250 random initial molecules
based on the eight different molecular frameworks. We use the BFGS algorithm
to optimize both observables. More details about the random initialization are
provided in the text.

Figure 3: Histograms of the optimized HOMO-LUMO gap, Eq. 3, for 250 random
initial molecules ($\epsilon_{HL}^{\mathrm{initial}}$). The inset panels
compare the similarity between $\epsilon_{HL}$ computed with the value of
$\mathbf{b}$ at the end of the optimization protocol
($\epsilon_{HL}^{\mathrm{virtual}}$), and the values of $\epsilon_{HL}$
selecting the most probable atoms given $\mathbf{b}^{*}$
($\epsilon_{HL}^{\mathrm{feasible}}$). Curves represent the derived histograms
using kernel density estimation, (solid) $\epsilon_{HL}^{\mathrm{feasible}}$,
(dashed) $\epsilon_{HL}^{\mathrm{virtual}}$, and (dotdash)
$\epsilon_{HL}^{\mathrm{initial}}$. All molecules were optimized using the
BFGS algorithm. Molecular frameworks are displayed in Fig. 1.

Figure 4: Histograms of the optimized $\langle\tilde{\alpha}\rangle$, Eq. 4,
for 250 random initial molecules
($\langle\tilde{\alpha}\rangle^{\mathrm{initial}}$). The inset panels compare
the similarity between $\langle\tilde{\alpha}\rangle$ computed with the value
of $\mathbf{b}$ at the end of the optimization protocol
($\langle\tilde{\alpha}\rangle^{\mathrm{virtual}}$), and the values of
$\langle\tilde{\alpha}\rangle$ selecting the most probable atoms given
$\mathbf{b}^{*}$ ($\langle\tilde{\alpha}\rangle^{\mathrm{feasible}}$). Curves
represent the derived histograms using kernel density estimation, (solid)
$\langle\tilde{\alpha}\rangle^{\mathrm{feasible}}$, (dashed)
$\langle\tilde{\alpha}\rangle^{\mathrm{virtual}}$, and (dotdash)
$\langle\tilde{\alpha}\rangle^{\mathrm{initial}}$. All molecules were
optimized using the BFGS algorithm. Molecular frameworks are displayed in Fig.
1.

Figure 5: Change in the parameters $\mathbf{b}$ for $\epsilon_{HL}$ (left
panel) and $\langle\tilde{\alpha}\rangle$ (right panel) during the
optimization of a single random initial virtual molecule based on framework 3.
For the initial, one intermediate and the final $\mathbf{b}$, we plot the
values of $\mathbf{b}$ for molecular framework 3. For both panels, the filled
markers represent the values of $\epsilon_{HL}$ and
$\langle\tilde{\alpha}\rangle$ computed with $\mathbf{b}$ at each iteration
(virtual), and the empty markers represent the observable values computed only
with the most probable atoms given $\mathbf{b}$ at each iteration (feasible).
For each search site we only considered three different atom types, namely
$C$, $N$ and $P$;
$\mathbf{b}_{\ell}=[b_{\ell}^{C},b_{\ell}^{N},b_{\ell}^{P}]$. The optimization
of $\epsilon_{HL}$ and $\langle\tilde{\alpha}\rangle$ w.r.t. $\mathbf{b}$ was
carried out with Adam using a learning rate of 0.2. More details are provided
in the main text.

Figure 6: The molecular structures with the lowest HOMO-LUMO gap and maximum
polarizability for each of the eighth molecules considered (Fig. 1).

###  3.2 Parameter optimization

Another important task that might sometimes be underappreciated in
computational chemistry is model parameter optimization. Here, we leverage the
flexibility of AD and optimize all free parameters of the Hückel model in the
same way as it is done for modern ML algorithms. Originally, the Hückel model
is solely based on electronic interactions between nearest-neighbour atoms,
which is typically also referred to as the tight-binding approximation (cf.
$\beta_{\ell,k}$ parameters in Eq. 1). Beyond the standard Hückel model, one
can introduce atomic distance-dependence of the corresponding interactions via
$\beta_{\ell,k}=\beta^{0}_{\ell,k}g(\mathbf{R}_{\ell,k})$. For example, based
on previous work by Longuett-Higgins and Salem [61], $\beta_{\ell,k}$ has an
exponential dependence on $\mathbf{R}_{\ell,k}$,

$\beta_{\ell,k}^{exp}=-\beta^{0}_{\ell,k}\exp^{-\frac{\Delta
R_{\ell,k}}{y_{\ell,k}}}.$ (8)

A second functional form, which is based on the work of Su, Schrieffer and
Heege [62, 63] on conducting polymers, uses a linear distance-dependence of
the interactions,

$\beta_{\ell,k}^{lr}=-\beta^{0}_{\ell,k}\left(1-y^{-1}_{\ell,k}\Delta
R_{\ell,k}\right).$ (9)

For both expressions (see Eqs. 8–9), $\Delta R_{\ell,k}$ is the difference
with respect to the reference bond length distance $R^{0}_{\ell,k}$, and
$y_{\ell,k}$ is a length scale parameter. By including $R^{0}_{\ell,k}$ and
$y_{\ell,k}$ in the set of parameters for the Hückel model, the complete set
of parameters becomes
$\boldsymbol{\theta}=[\alpha_{\ell},\beta^{0}_{\ell,k},y_{\ell,k},R^{0}_{\ell,k}]$.

For this work, all initial $\alpha_{\ell}$ and $\beta^{0}_{\ell,k}$ parameters
were taken from Van-Catledge [10], and the initial $R^{0}_{\ell,k}$ parameters
were approximated from tables of standard bond lengths [64]. The length scale
parameters ($y_{\ell,k}$) were initially set to 0.3 Å, which corresponds to
the value that has been used for C–C in the literature [65, 66].

We used a subset of the GDB-13 data [67] set that only consists of molecules
with $\pi$-systems for fitting our model parameters (see Supplementary
Material for details on how the dataset was generated). Note that some
molecules in the dataset could have n-$\pi$* transitions as their lowest
excited state. We used a pool of 60,000 molecules and randomly sampled 100,
1,000, and 5,000 molecules from this set, and used 1,000 additional molecules
as validation set to monitor the optimization procedure. To optimize
$\boldsymbol{\theta}$, we used the mean squared error as a loss function,

${\cal
L}(\boldsymbol{\theta})=\frac{1}{2}\sum_{i}^{N}\left(\hat{\epsilon}_{HL}({\cal
M}_{i})-\epsilon_{HL}(\boldsymbol{\theta};{\cal M}_{i})\right)^{2},$ (10)

where ${\cal M}_{i}$ is a single molecule of the training set, and
$\hat{\epsilon}_{HL}$ is the vertical excitation energy between the ground
state and the first excited singlet state computed at the TDA-
SCS-$\omega$PBEPP86/def2-SVP level of theory [68, 69]. At the Hückel level of
theory, this excitation energy simply corresponds to the HOMO-LUMO gap
($\epsilon_{HL}$) due to the disregard for electron correlation. To compare
the prediction of $\epsilon_{HL}$ with the DFT reference values properly, we
linearly transformed the results of the Hückel model using two additional
parameters, $w_{0}$ and $w_{1}$ (Eq. 11),

$\epsilon_{HL}(\boldsymbol{\theta};{\cal
M}_{i})=w_{1}\times\epsilon_{HL}(\alpha_{\ell},\beta^{0}_{\ell,k},y_{\ell,k},R^{0}_{\ell,k};{\cal
M}_{i})+w_{0},$ (11)

where $\boldsymbol{\theta}$ jointly represents all parameters of the model,
i.e.,
$\boldsymbol{\theta}=[\alpha_{\ell},\beta^{0}_{\ell,k},y_{\ell,k},R^{0}_{\ell,k},w_{0},w_{1}]$.

For the optimization of all free parameters, we used the AdamW optimization
algorithm [46], as implemented in the Optax library [70], with a learning rate
of $0.02$, and a weight decay of $10^{-4}$. Notably, we considered various
training scenarios that included different values for the weight decay, and
the regularization of different sets of the Hückel parameters. However, we
found no impact on the accuracy of the model. The initial model parameters
were gathered from Refs. [10, 64].

We optimized the parameters of three different Hückel models, i) the original
one where $\beta_{\ell,k}$ is distance-independent
($\beta_{\ell,k}=\beta^{0}_{\ell,k}$), and both ii) the exponential (Eq. 8),
and iii) the linear (Eq. 9) distance-dependence functional forms. We want to
emphasize that any other analytic form for $\beta_{\ell,k}$ could be
considered as well as AD makes any of these expressions fully differentiable.
Following the convention in the literature, we scaled the parameters
$\beta^{0}_{\ell,k}$ and $\alpha_{\ell}$ with respect to the carbon atom
parameters according to $\alpha_{\ell}=\alpha_{\ell}-\alpha_{C}$, and
$\beta^{0}_{\ell,k}=\beta^{0}_{\ell,k}/\beta^{0}_{C,C}$. Notably, at least for
our results in this work, we found that including a regularization term in the
loss function did not impact the accuracy of the model. Finally, we found 20
epochs to be enough to minimize the loss function when the parameters are
initialized with values from Refs. [10, 64].

In Figures 7–10 we display the optimized values of the parameters for the
three different Hückel models considered. From the optimized parameters, we
observe that $\alpha_{O}$ (i.e., the 2p orbital energy parameter for oxygen),
for all three models, is the one that differs the most from the literature
[10, 64]. While there is no good reference data for $y_{\ell,k}$ to compare
to, we observe, nevertheless, that the C–C parameter value changes
considerably from the initial value of 0.3. For the $\beta^{0}_{\ell,k}$
parameters, only the values for N–C resemble the literature values.
Furthermore, the optimal values of $R^{0}_{\ell,k}$ change the least from the
values in Ref. [64].

Using these optimized parameters, we predicted $\epsilon_{HL}$ for 40,000
additional test set molecules and compare the results with DFT reference data.
The direct comparison is depicted in Figures 11–13. By optimizing all
parameters there is a significant improvement in the prediction of
$\epsilon_{HL}$ using our semi-empirical model. Notably, we also found that
considering a larger training data set does not impact the accuracy of the
Hückel model which suggests that either the corresponding molecules do not
provide any additional information with respect to the relevant interaction
parameters or that the model already is close to its best expected performance
and cannot be improved further. Another important observation in that regard
is that the analytical form of $g(\mathbf{R}_{\ell,k})$ in $\beta_{\ell,k}$
does not impact the accuracy of the model when optimized parameters are used.

Next, we also optimized the parameters with respect to the polarizability for
the distance-independent Hückel model. Here, the gradients needed for training
are of third-order, e.g.,
$\frac{\partial}{\partial\alpha_{\ell}}\frac{\partial^{2}\tilde{\alpha}_{kk}}{\partial
F^{2}_{k}}$ or
$\frac{\partial}{\partial\beta_{k,\ell}}\frac{\partial^{2}\tilde{\alpha}_{kk}}{\partial
F^{2}_{k}}$, and can be computed more efficiently via AD, illustrating the
potential of this approach. The molecular polarizabilities that were used as
reference data were computed using dftd4 (version 3.4.0) [71, 72, 73] via the
default methodology summing atomic polarizabilities. Even though we observe a
higher prediction accuracy when the optimized parameters are used compared to
the model before parameter refinement, as depicted in Fig. 14, the accuracy of
the model still remains relatively low and does not improve anymore when more
training data is used. We suspect that this prediction task is particularly
challenging for the Huc̈kel model as the simulated polarizability only
corresponds to the contribution from $\pi$-electrons, while that of the
reference data accounts for all the electrons in the molecules. Even though we
expect a significant portion of the molecular polarizabilities to stem from
the $\pi$-electrons, the contributions of the $\sigma$-electrons cannot be
neglected and can dominate this property. Nevertheless, this proof-of-concept
application example demonstrates the operational ease of conducting parameter
refinement of a given physics-based prediction model based on reference data,
even when derivative properties are targeted.

Figure 7: Optimized $\alpha_{\ell}$ parameters for different Hückel models, a)
$\beta^{0}_{\ell,k}$, b) $\beta^{lr}_{\ell,k}$ (Eq. 9), and c)
$\beta^{exp}_{\ell,k}$ (Eq. 8). All parameters reported were averaged over ten
different training data sets, and different number of training molecules.
Colored symbols and bars represent the mean and standard deviation of the
optimized parameters averaged over ten different data sets. The reference
parameters ($\blacksquare$-symbol) were taken from Refs. [10], and used as the
initial parameters for all optimizations. We refer the reader to the main text
for the optimization details.

Figure 8: Optimized $\beta^{0}_{\ell,k}$ parameters for different Hückel
models, a) $\beta^{0}_{\ell,k}$, b) $\beta^{lr}_{\ell,k}$ (Eq. 9), and c)
$\beta^{exp}_{\ell,k}$ (Eq. 8). All parameters reported were averaged over ten
different training data sets, and different number of training molecules.
Colored symbols and bars represent the mean and standard deviation of the
optimized parameters averaged over ten different data sets. The reference
parameters ($\blacksquare$-symbol) were taken from Refs. [10], and used as the
initial parameters for all optimizations. We refer the reader to the main text
for the optimization details.

Figure 9: Optimized $y_{\ell,k}$ parameters for different Hückel models, a)
$\beta^{lr}_{\ell,k}$ (Eq. 9), and b) $\beta^{exp}_{\ell,k}$ (Eq. 8). All
parameters reported were averaged over ten different training data sets, and
different number of training molecules. Colored symbols and bars represent the
mean and standard deviation of the optimized parameters averaged over ten
different data sets.The reference parameters ($\blacksquare$-symbol) were set
to 0.3 Å, and used as the initial parameters for all optimizations. We refer
the reader to the main text for the optimization details.

Figure 10: Optimized $R^{0}_{\ell,k}$ parameters for different Hückel models,
a) $\beta^{lr}_{\ell,k}$ (Eq. 9), and b) $\beta^{exp}_{\ell,k}$ (Eq. 8). All
parameters reported were averaged over ten different training data sets, and
different number of training molecules. Colored symbols and bars represent the
mean and standard deviation of the optimized parameters averaged over ten
different data sets. The reference parameters ($\blacksquare$-symbol) were
taken from Refs. [64], and used as the initial parameters for all
optimizations. We refer the reader to the main text for the optimization
details.

Figure 11: $\epsilon_{HL}$ predicted with the Hückel models and with DFT level
for 40K molecules not considered during training. The atom-atom interaction of
the Hückel model is described by a distance-independent parameter,
$\beta^{0}_{\ell,k}$. Results computed with a Hückel with parameters taken
from the literature (a), and parameters optimized with N=100 (b) and N=5,000
(c) data points. We refer the reader to the main text for the optimization
details.

Figure 12: $\epsilon_{HL}$ predicted with the Hückel models and with DFT level
for 40K molecules not considered during training. The atom-atom interaction of
the Hückel model is described by a distance-dependent parameter,
$\beta^{lr}_{\ell,k}$ (Eq. 9). Results computed with a Hückel with parameters
taken from the literature (a), and parameters optimized with N=100 (b) and
N=5,000 (c) data points. We refer the reader to the main text for the
optimization details.

Figure 13: $\epsilon_{HL}$ predicted with the Hückel models and with DFT level
for 40K molecules not considered during training. The atom-atom interaction of
the Hückel model is described by a distance-dependent parameter,
$\beta^{exp}_{\ell,k}$ (Eq. 8). Results computed with a Hückel model with
parameters taken from the literature (a), and parameters optimized with N=100
(b) and N=5,000 (c) data points. We refer the reader to the main text for the
optimization details.

Figure 14: $\langle\tilde{\alpha}\rangle$ predicted with the Hückel models and
with DFT level for 40K molecules not considered during training. The atom-atom
interaction of the Hückel model is described by a distance-independent
parameter, $\beta^{0}_{\ell,k}$. Results computed with a Hückel model with
parameters taken from the literature (a), and (b) parameters optimized with
N=100 data points. We refer the reader to the main text for the optimization
details.

##  4 Conclusions

In this work, we demonstrate the power of automatic differentiation to enable
the efficient use of physics-inspired models for gradient-based optimization
problems in the realm of molecular chemistry via semi-empirical Hückel models.
In particular, we showcase inverse molecular design via an alchemical problem
formulation using fixed molecular frameworks. This allows us to perform
structure optimization requiring only a very small number of intermediate
structures to find local minima with respect to the properties of interest
utilizing gradients with respect to atom identities at specific sites. While
our approach is currently limited to a fixed molecular framework, performing
optimizations over the molecular composition space alone is far from trivial.
Compared to various alternative approaches, our implementation shows a
remarkably high molecular sampling efficiency due to efficient utilization of
gradient information in combination with powerful gradient-based optimization
algorithms. Additionally, we showcase the ease of generating calibrated
physics-based property prediction models using high quality reference training
data of relatively modest size, again allowing for quick convergence of model
parameters. This is particularly important as most physical models that rely
on empirical parameters such as semi-empirical quantum chemistry models and
density functional approximations are still largely optimized by hand, making
the corresponding procedures tedious. Thus, we believe that our work will
serve as an inspiration for the field of computational chemistry in order to
adopt the readily available AD capabilities of mature ML programming
frameworks allowing to accelerate the construction of ever more accurate
physics-based property simulation models.

## Acknowledgments

R.P. acknowledges funding through a Postdoc.Mobility fellowship by the Swiss
National Science Foundation (SNSF, Project No. 191127). A.A.-G. thanks Anders
G. Frøseth for his generous support. A.A.-G. acknowledges the generous support
of Natural Resources Canada and the Canada 150 Research Chairs program. We
also thank the SciNet HPC Consortium for support regarding the use of the
Niagara supercomputer. SciNet is funded by the Canada Foundation for
Innovation, the Government of Ontario, Ontario Research Fund - Research
Excellence, and the University of Toronto.

## Data Availability

The data that support the findings of this study are available within the
article and its supplementary material and in GitHub (inverse molecular
design) https://github.com/RodrigoAVargasHdz/huxel_molecule_desing and
(parameter optimization) https://github.com/RodrigoAVargasHdz/huxel.

## References

## References

  * [1] R. S. Mulliken, “Molecular scientists and molecular science: Some reminiscences,” Journal of Chemical Physics, vol. 43, pp. S2–S11, Nov. 1965.
  * [2] P.-O. Löwdin, “Present situation of quantum chemistry,” Journal of Physical Chemistry, vol. 61, pp. 55–68, Jan. 1957.
  * [3] W. Thiel, “Semiempirical quantum–chemical methods,” WIREs Computational Molecular Science, vol. 4, no. 2, pp. 145–157, 2014.
  * [4] E. Hückel, “Quantentheoretische beiträge zum benzolproblem,” Zeitschrift für Physik, vol. 70, pp. 204–286, Mar 1931.
  * [5] E. Hückel, “Quanstentheoretische beiträge zum benzolproblem,” Zeitschrift für Physik, vol. 72, pp. 310–337, May 1931.
  * [6] E. Hückel, “Quantentheoretische beiträge zum problem der aromatischen und ungesättigten verbindungen. iii,” Zeitschrift für Physik, vol. 76, pp. 628–648, Sep 1932.
  * [7] E. Hückel, “Die freien radikale der organischen chemie,” Zeitschrift für Physik, vol. 83, pp. 632–668, Sep 1933.
  * [8] A. Streitwieser, Molecular Orbital Theory for Organic Chemists.  New York,: Wiley, 1961.
  * [9] B. Andes Hess, L. Schaad, and C. Holyoke, “On the aromaticity of annulenones,” Tetrahedron, vol. 28, pp. 5299–5305, Jan. 1972.
  * [10] F. A. Van-Catledge, “A Pariser-Parr-Pople-based set of Hückel molecular orbital parameters,” The Journal of Organic Chemistry, vol. 45, pp. 4801–4802, Nov. 1980.
  * [11] I. Fleming, Molecular Orbitals and Organic Chemical Reactions: Reference Edition.  Chichester, United Kingdom: Wiley, 2010.
  * [12] P. Klán and J. Wirz, Photochemistry of Organic Compounds: From Concepts to Practice.  Wiley, 2009.
  * [13] C. S. Anstöter and J. R. R. Verlet, “A Hückel model for the excited-sate dynamics of a protein chromophore developed using photoelectron imaging,” Accounts of Chemical Research, vol. 55, no. 9, pp. 1205–1213, 2022.
  * [14] J. A. N. F. Gomes and R. B. Mallion, “Aromaticity and ring currents,” Chemical Reviews, vol. 101, no. 5, pp. 1349–1384, 2001.
  * [15] B. Sanchez-Lengeling and A. Aspuru-Guzik, “Inverse molecular design using machine learning: Generative models for matter engineering,” Science, vol. 361, pp. 360–365, July 2018.
  * [16] D. Xiao, W. Yang, and D. N. Beratan, “Inverse molecular design in a tight-binding framework,” The Journal of Chemical Physics, vol. 129, no. 4, p. 044106, 2008.
  * [17] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, “Automatic differentiation in machine learning: a survey,” J. Mach. Learn. Res., vol. 18, 2018.
  * [18] T. Tamayo-Mendoza, C. Kreisbeck, R. Lindh, and A. Aspuru-Guzik, “Automatic differentiation in quantum chemistry with applications to fully variational hartree–fock,” ACS Central Science, vol. 4, no. 5, pp. 559–566, 2018\.
  * [19] X. Zhang and G. K.-L. Chan, “Differentiable quantum chemistry with pyscf for molecules and materials at the mean-field level and beyond,” 2022.
  * [20] M. F. Kasim, S. Lehtola, and S. M. Vinko, “Dqc: a python program package for differentiable quantum chemistry,” 2021.
  * [21] L. Li, S. Hoyer, R. Pederson, R. Sun, E. D. Cubuk, P. Riley, and K. Burke, “Kohn-sham equations as regularizer: Building prior knowledge into machine-learned physics,” Phys. Rev. Lett., vol. 126, p. 036401, Jan 2021\.
  * [22] M. F. Kasim and S. M. Vinko, “Learning the exchange-correlation functional from nature with fully differentiable density functional theory,” Phys. Rev. Lett., vol. 127, p. 126403, Sep 2021.
  * [23] L. Zhao and E. Neuscamman, “Excited state mean-field theory without automatic differentiation,” The Journal of Chemical Physics, vol. 152, no. 20, p. 204112, 2020.
  * [24] G. Zhou, B. Nebgen, N. Lubbers, W. Malone, A. M. N. Niklasson, and S. Tretiak, “Graphics processing unit-accelerated semiempirical Born Oppenheimer molecular dynamics using PyTorch,” Journal of Chemical Theory and Computation, vol. 16, pp. 4951–4962, Aug. 2020.
  * [25] T. Zubatiuk, B. Nebgen, N. Lubbers, J. S. Smith, R. Zubatyuk, G. Zhou, C. Koh, K. Barros, O. Isayev, and S. Tretiak, “Machine learned hückel theory: Interfacing physics and deep neural networks,” The Journal of Chemical Physics, vol. 154, no. 24, p. 244108, 2021.
  * [26] C. H. Pham, R. K. Lindsey, L. E. Fried, and N. Goldman, “High-accuracy semiempirical quantum models based on a minimal training set,” The Journal of Physical Chemistry Letters, vol. 13, no. 13, pp. 2934–2942, 2022\.  PMID: 35343698.
  * [27] G. Zhou, N. Lubbers, K. Barros, S. Tretiak, and B. Nebgen, “Deep learning of dynamically responsive chemical hamiltonians with semiempirical quantum mechanics,” Proceedings of the National Academy of Sciences, vol. 119, no. 27, p. e2120333119, 2022.
  * [28] A. Yachmenev and S. N. Yurchenko, “Automatic differentiation method for numerical construction of the rotational-vibrational hamiltonian as a power series in the curvilinear internal coordinates using the eckart frame,” The Journal of Chemical Physics, vol. 143, no. 1, p. 014105, 2015.
  * [29] A. S. Abbott, B. Z. Abbott, J. M. Turney, and H. F. Schaefer, “Arbitrary-order derivatives of quantum chemical methods via automatic differentiation,” The Journal of Physical Chemistry Letters, vol. 12, no. 12, pp. 3232–3239, 2021\.  PMID: 33764068.
  * [30] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, M. S. Alam, S. Ahmed, J. M. Arrazola, C. Blank, A. Delgado, S. Jahangiri, K. McKiernan, J. J. Meyer, Z. Niu, A. Száva, and N. Killoran, “Pennylane: Automatic differentiation of hybrid quantum-classical computations,” 2020.
  * [31] F. Pavošević and S. Hammes-Schiffer, “Automatic differentiation for coupled cluster methods,” 2020.
  * [32] G. F. von Rudorff, “Arbitrarily accurate quantum alchemy,” The Journal of Chemical Physics, vol. 155, no. 22, p. 224103, 2021.
  * [33] N. Fedik, R. Zubatyuk, M. Kulichenko, N. Lubbers, J. S. Smith, B. Nebgen, R. Messerly, Y. W. Li, A. I. Boldyrev, K. Barros, O. Isayev, and S. Tretiak, “Extending machine learning beyond interatomic potentials for predicting molecular properties,” Nature Reviews Chemistry, Aug 2022.
  * [34] R. A. Vargas-Hernández, R. T. Q. Chen, K. A. Jung, and P. Brumer, “Fully differentiable optimization protocols for non-equilibrium steady states,” New Journal of Physics, vol. 23, p. 123006, dec 2021.
  * [35] R. A. Vargas-Hernández, R. T. Q. Chen, K. A. Jung, and P. Brumer, “Inverse design of dissipative quantum steady-states with implicit differentiation,” 2020\.
  * [36] N. Yoshikawa and M. Sumita, “Automatic differentiation for the direct minimization approach to the hartree–fock method,” The Journal of Physical Chemistry A, vol. 126, no. 45, pp. 8487–8493, 2022.  PMID: 36346835.
  * [37] S. Bac, A. Patra, K. J. Kron, and S. Mallikarjun Sharada, “Recent advances toward efficient calculation of higher nuclear derivatives in quantum chemistry,” The Journal of Physical Chemistry A, vol. 126, no. 43, pp. 7795–7805, 2022.  PMID: 36282088.
  * [38] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” Nature, vol. 585, pp. 357–362, Sept. 2020.
  * [39] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, and S. Wanderman-Milne, “JAX: composable transformations of Python+NumPy programs,” 2018.
  * [40] B. Ramsundar, D. Krishnamurthy, and V. Viswanathan, “Differentiable physics: A position piece,” 2021.
  * [41] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015.  Software available from tensorflow.org.
  * [42] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” 2019.
  * [43] K. YATES, “Ii - hÜckel molecular orbital theory,” in Hückel Molecular Orbital Theory (K. YATES, ed.), pp. 27–87, Academic Press, 1978.
  * [44] J. Pople and D. Beveridge, Approximate Molecular Orbital Theory.  McGraw-Hill series in advanced chemistry, McGraw-Hill, 1970.
  * [45] M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-López, F. Pedregosa, and J.-P. Vert, “Efficient and modular implicit differentiation,” arXiv preprint arXiv:2105.15183, 2021.
  * [46] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” in International Conference on Learning Representations, 2019.
  * [47] F. Eisenlohr, “Eine neuberechnung der atomrefraktionen. i,” Zeitschrift für Physikalische Chemie, vol. 75U, no. 1, pp. 585–607, 1911.
  * [48] A. L. von Steiger, “Ein beitrag zur summationsmethodik der molekularrefraktionen, besonders bei aromatischen kohlenwasserstoffen,” Berichte der deutschen chemischen Gesellschaft (A and B Series), vol. 54, no. 6, pp. 1381–1393, 1921.
  * [49] P. Schwerdtfeger and J. K. Nagle, “2018 table of static dipole polarizabilities of the neutral elements in the periodic table,” Molecular Physics, vol. 117, no. 9-12, pp. 1200–1225, 2019.
  * [50] A. Krawczuk, D. Pérez, and P. Macchi, “PolaBer: a program to calculate and visualize distributed atomic polarizabilities based on electron density partitioning,” Journal of Applied Crystallography, vol. 47, pp. 1452–1458, Aug 2014.
  * [51] E. Heilbronner and H. Bock, The HMO Model and Its Application: Basis and Manipulation.  Wiley, 1976.
  * [52] W. Gao, T. Fu, J. Sun, and C. W. Coley, “Sample efficiency matters: A benchmark for practical molecular optimization,” arXiv preprint arXiv:2206.12411, 2022.
  * [53] F. Häse, L. M. Roch, and A. Aspuru-Guzik, “Chimera: enabling hierarchy based multi-objective optimization for self-driving laboratories,” Chem. Sci., vol. 9, pp. 7642–7655, 2018.
  * [54] R. J. Hickman, M. Aldeghi, F. Häse, and A. Aspuru-Guzik, “Bayesian optimization with known experimental and design constraints for chemistry applications,” 2022.
  * [55] M. Seifrid, R. J. Hickman, A. Aguilar-Granda, C. Lavigne, J. Vestfrid, T. C. Wu, T. Gaudin, E. J. Hopkins, and A. Aspuru-Guzik, “Routescore: Punching the ticket to more efficient materials development,” ACS Central Science, vol. 8, no. 1, pp. 122–131, 2022.
  * [56] H. Tao, T. Wu, S. Kheiri, M. Aldeghi, A. Aspuru-Guzik, and E. Kumacheva, “Self-driving platform for metal nanoparticle synthesis: Combining microfluidics and machine learning,” Advanced Functional Materials, vol. 31, no. 51, p. 2106725, 2021.
  * [57] R. A. Vargas-Hernández, C. Chuang, and P. Brumer, “Multi-objective optimization for retinal photoisomerization models with respect to experimental observables,” The Journal of Chemical Physics, vol. 155, no. 23, p. 234109, 2021.
  * [58] J.-A. Désidéri, “Multiple-gradient descent algorithm (mgda) for multiobjective optimization,” Comptes Rendus Mathematique, vol. 350, no. 5, pp. 313–318, 2012.
  * [59] G. F. von Rudorff and O. A. von Lilienfeld, “Alchemical perturbation density functional theory,” Phys. Rev. Research, vol. 2, p. 023220, May 2020.
  * [60] C. Shen, M. Krenn, S. Eppel, and A. Aspuru-Guzik, “Deep molecular dreaming: inverse machine learning for de-novo molecular design and interpretability with surjective representations,” Machine Learning: Science and Technology, vol. 2, p. 03LT02, jul 2021.
  * [61] H. Longuet-Higgins and L. Salem, “The alternation of bond lengths in long conjugated chain molecules,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, vol. 251, pp. 172–185, May 1959.
  * [62] W. P. Su, J. R. Schrieffer, and A. J. Heeger, “Solitons in polyacetylene,” Phys. Rev. Lett., vol. 42, pp. 1698–1701, Jun 1979.
  * [63] A. J. Heeger, S. Kivelson, J. R. Schrieffer, and W. P. Su, “Solitons in conducting polymers,” Rev. Mod. Phys., vol. 60, pp. 781–850, Jul 1988\.
  * [64] F. H. Allen, O. Kennard, D. G. Watson, L. Brammer, A. G. Orpen, and R. Taylor, “Tables of bond lengths determined by x-ray and neutron diffraction. part 1. bond lengths in organic compounds,” J. Chem. Soc., Perkin Trans. 2, pp. S1–S19, 1987.
  * [65] L. Z. Stolarczyk and T. M. Krygowski, “Augmented hückel molecular orbital model of $\pi$-electron systems: from topology to metric. i. general theory,” Journal of Physical Organic Chemistry, vol. 34, no. 3, p. e4154, 2021.  e4154 poc.4154.
  * [66] J. H. Kwapisz and L. Z. Stolarczyk, “Applications of hückel-su-schrieffer-heeger method,” Structural Chemistry, vol. 32, pp. 1393–1406, Aug 2021.
  * [67] L. C. Blum and J.-L. Reymond, “970 million druglike small molecules for virtual screening in the chemical universe database gdb-13,” Journal of the American Chemical Society, vol. 131, no. 25, pp. 8732–8733, 2009.
  * [68] M. Casanova-Páez and L. Goerigk, “Time-dependent long-range-corrected double-hybrid density functionals with spin-component and spin-opposite scaling: A comprehensive analysis of singlet–singlet and singlet–triplet excitation energies,” Journal of Chemical Theory and Computation, vol. 17, no. 8, pp. 5165–5186, 2021.  PMID: 34291643.
  * [69] F. Weigend and R. Ahlrichs, “Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for h to rn: Design and assessment of accuracy,” Phys. Chem. Chem. Phys., vol. 7, pp. 3297–3305, 2005.
  * [70] M. Hessel, D. Budden, F. Viola, M. Rosca, E. Sezener, and T. Hennigan, “Optax: composable gradient transformation and optimisation, in jax!,” 2020.
  * [71] E. Caldeweyher, C. Bannwarth, and S. Grimme, “Extension of the d3 dispersion coefficient model,” The Journal of Chemical Physics, vol. 147, no. 3, p. 034112, 2017.
  * [72] E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth, and S. Grimme, “A generally applicable atomic-charge dependent london dispersion correction,” The Journal of Chemical Physics, vol. 150, no. 15, p. 154122, 2019.
  * [73] E. Caldeweyher, J.-M. Mewes, S. Ehlert, and S. Grimme, “Extension and evaluation of the d4 london-dispersion model for periodic systems,” Phys. Chem. Chem. Phys., vol. 22, pp. 8499–8512, 2020.