Inverse molecular design and parameter optimization with H\"uckel theory using automatic differentiation

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\"uckel's $\pi$-electron molecular orbital theory. In this work, we implemented a H\"uckel program using differentiable programming with the JAX framework, based on limited modifications of a pre-existing NumPy version. The auto-differentiable H\"uckel 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.


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 nontrivial. 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 π-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 meanfield 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 10 5 -10 6 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).

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 (L) with respect to the model parameters (θ), ∇ θ L(θ). All contemporary ML libraries, e.g., Tensorflow [41], PyTorch [42] and JAX [39], are built on top of an AD engine which computes ∇ θ L(θ) for any ML model. Given the robustness of AD libraries, differentiating physical models [40] could be done similarly to modern ML algorithms.

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 π-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 α and β ,k , respectively. The matrix elements of the Hückel Hamiltonian are given by where the α parameters roughly represent the energy of an electron in a 2p orbital, and the β ,k parameters describe the energy of an electron in the bond −k. Extensions of the Hückel Hamiltonian are possible and can, for instance, incorporate distance-dependence via β ,k = β 0 ,k g(R ,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 α and β ,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.

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, where b i is the weight of the atom of type i for site . For a meaningful description, the weights of each site must be normalized, i.e., M i b 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 π-framework (Fig. 1). Therefore we defined the following vector of atom type weight parameters: For clarity, b jointly describes the b parameters for all search sites in a molecule, i.e., b = {b } N . For all results presented, the values of the α and β 0 ,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 b at each site that gives the lowest HOMO-LUMO gap (Eq. 3), denoted as HL , and the maximum polarizability denoted as α , (Eq, 4). HL is defined as, where HOM O and LU M O are the eigenvalues of the highest occupied molecular orbital (HOMO), and lowest unoccupied molecular orbital (LUMO), respectively. The polarizability function is defined as where theα ij elements are the polarizability components defined as The F i terms are the components of the electric field, 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], 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 b are to be optimized using a gradient-based method, the Jacobians ∇ b HL and ∇ b α will also be constructed using an FD approach. However, this will increase the number of energy evaluations needed, especially for α as the elements of ∇ b α are third-order derivatives: Thus, for a single element of ∇ b α , using FD will require 18 energy calculations, In contrast, using modern AD frameworks, we can efficiently compute, with a single energy calculation (i.e., one forward pass), the Jacobian of HL with respect to b. For α , 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 α with respect to 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, HL and α . 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 b, we used the softmax function, where 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 HL or α . Notably, we initialized the values for all b i parameters by sampling a uniform distribution, b i ∼ 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 b allows us to sample a wide range of molecules with a broad range of values for both objectives, initial HL and α initial . Because of the statistical description of the molecules by the b parameters, the optimal parameters (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 virtual and that for the real molecule as y 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 b throughout the optimization for both objectives ( HL and α ). As we can observe, the change in 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 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 b (Figs. 3 and 4). We also observe that, for the majority of the optimized molecules, y virtual and y f easible are linearly correlated, indicating that the optimization converged essentially to feasible molecules. The property distributions of y virtual and y feasible are reasonably close, even in the few cases when the correlations are poor. Fig. 6 displays the molecules with the lowest HL and maximum α from the ensemble of different optimizations. First, we notice that there is a higher amount of phosphorus atoms in the molecules when α 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 HL , we see extensive incorporation of both nitrogen and phosphorous atoms. This can be understood in terms of the effect of heteroatom substitution on HL within the HMO framework [51]. For alternant compounds such as 1-8, HL is unaffected by changing α . The main effect comes from changing β ,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 β ,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 β ,k leads to a lowering of the energy. The net effect by raising the HOMO and lowering the LUMO is a decrease in the HL gap. We can therefore expect optimization to favor atom pairs with a low β ,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 β ,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 β 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 HL and α . 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 X X X X X X X X X 1 2 4 3 5 8 7 6 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. 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.  . Histograms of the optimized α , Eq. 4, for 250 random initial molecules ( α initial ). The inset panels compare the similarity between α computed with the value of b at the end of the optimization protocol ( α virtual ), and the values of α selecting the most probable atoms given b * ( α feasible ). Curves represent the derived histograms using kernel density estimation, (solid) α feasible , (dashed) α virtual , and (dotdash) α initial . All molecules were optimized using the BFGS algorithm. Molecular frameworks are displayed in Fig. 1.

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. β ,k parameters in Eq. 1). Beyond the standard Hückel model, one can introduce atomic distance-dependence of the corresponding interactions via β ,k = β 0 ,k g(R ,k ). For example, based on previous work by Longuett-Higgins and Salem [61], β ,k has an exponential dependence on R ,k , 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, For both expressions (see Eqs. [8][9], ∆R ,k is the difference with respect to the reference bond length distance R 0 ,k , and y ,k is a length scale parameter. By including R 0 ,k and y ,k in the set of parameters for the Hückel model, the complete set of parameters becomes θ = [α , β 0 ,k , y ,k , R 0 ,k ]. For this work, all initial α and β 0 ,k parameters were taken from Van-Catledge [10], and the initial R 0 ,k parameters were approximated from tables of standard bond lengths [64]. The length scale parameters (y ,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 π-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-π* 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 θ, we used the mean squared error as a loss function, where M i is a single molecule of the training set, andˆ HL is the vertical excitation energy between the ground state and the first excited singlet state computed at the TDA-SCS-ω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 ( HL ) due to the disregard for electron correlation. To compare the prediction of 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), where θ jointly represents all parameters of the model, i.e., θ = [α , β 0 ,k , y ,k , R 0 ,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 β ,k is distance-independent (β ,k = β 0 ,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 β ,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 β 0 ,k and α with respect to the carbon atom parameters according to α = α − α C , and β 0 ,k = β 0 ,k /β 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 α 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 ,k to compare to, we observe, nevertheless, that the C-C parameter value changes considerably from the initial value of 0.3. For the β 0 ,k parameters, only the values for N-C resemble the literature values. Furthermore, the optimal values of R 0 ,k change the least from the values in Ref. [64].
Using these optimized parameters, we predicted 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 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(R ,k ) in β ,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., , 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 Huckel model as the simulated polarizability only corresponds to the contribution from π-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 π-electrons, the contributions of the σ-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. ,k , b) β lr ,k (Eq. 9), and c) β exp ,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 ( -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. ,k parameters for different Hückel models, a) β 0 ,k , b) β lr ,k (Eq. 9), and c) β exp ,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 ( -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 ,k parameters for different Hückel models, a) β lr ,k (Eq. 9), and b) β exp ,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 ( -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 ,k parameters for different Hückel models, a) β lr ,k (Eq. 9), and b) β exp ,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 ( -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. 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, β 0 ,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. 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, β lr ,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.

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 Figure 13. 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, β exp ,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. α 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, β 0 ,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.
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. Table S1: List of filters employed to create the π-systems subset of GDB-13. In each line, x denotes the value of the corresponding feature. On the structures obtained after filtering, we performed property simulations. To simulate the excited state properties of all the molecules, first, ground state conformational ensembles were generated using CREST [13] (version 2.11) with the iMTD-GC [6,12] workflow (default option) via the GFN2-xTB//GFN-FF [1,14] composite method. The lowest energy conformers were optimized using ORCA [10] (version 5.0.1) at the r2SCAN-3c [7] level of theory. The corresponding geometries were used for subsequent ground-and excited state single-point calculations to assess the properties of vertical excitations. Single points at the TDA-SCS-ωPBEPP86/def2-SVP [3, 16] level  Figure S2: HOMO and LUMO depictions with coefficients for the best optimized molecules using the Van-Catledge parameter set.