Virtual screening of norbornadiene-based molecular solar thermal energy storage systems using a genetic algorithm

We present a computational methodology for the screening of a chemical space of 1025 substituted norbornadiene molecules for promising kinetically stable molecular solar thermal (MOST) energy storage systems with high energy densities that absorb in the visible part of the solar spectrum. We use semiempirical tight-binding methods to construct a dataset of nearly 34,000 molecules and train graph convolutional networks to predict energy densities, kinetic stability, and absorption spectra and then use the models together with a genetic algorithm to search the chemical space for promising MOST energy storage systems. We identify 15 kinetically stable molecules, five of which have energy densities greater than 0.45 MJ/kg and the main conclusion of this study is that the largest energy density that can be obtained for a single norbornadiene moiety with the substituents considered here, while maintaining a long half-life and absorption in the visible spectrum, is around 0.55 MJ/kg.


Introduction
Today, the production of solar and wind energy have become more profitable than nonrenewable alternatives. However, daily and seasonal variations in the renewable energy production as well as large variations in the power demands are two serious challenges for a sustainable energy eco-system. Storing excess power and using this energy in peak times is therefore absolutely crucial, 1,2 but it requires affordable large-scale energy storage systems. One technology that attempts to solve this problem is closed-cycle MOlecular Solar Thermal (MOST) energy storage. [3][4][5][6][7][8][9][10][11][12][13] Such systems relies on photochromic molecules or molecular photoswitches where stable reactants can interconvert to form metastable products using solar irradiation as the driving force. Thus, the solar energy can be stored as chemical energy until a subsequent exothermic reaction releases the captured energy. Depending on the storage lifetime, this energy release can either occur spontaneously or be controlled by e.g. a thermal activation, a heterogeneous catalyst, an electric potential, or light. To paraphrase, MOST systems are able to harvest and store solar energy, which can later be released as clean thermal energy for space heating or heating of domestic water. In fact, the thermal energy can be stored without the need for thermal insulation, which enables prolonged storage times compared to alternative thermal energy storage systems.
One of the most promising MOST systems is the norbornadiene/quadricyclane (NBD/QC) couple (Figure 1), which was introduced as a potential MOST system in 1961 by Dauben and Cargill,14 although the isomerization reaction of a NBD/QC dicarboxylic acid system was already observed in 1954 by Cristol and Snell. 15,16 Since then, the system has been extensively studied due to its high energy density of almost 1 MJ/kg, which is estimated to be the fundamental upper limit of MOST systems and comparable to Li-ion batteries. 7 However, the absorption spectrum of NBD lies in the UV region with absorption onset at 267 nm, and has therefore no overlap with the spectrum of solar radiation. 12 Several studies have shown that it is possible to obtain a large redshift in the absorption spectrum and better quantum yield, but at the expense of a decrease in the energy density due to the added weight of the chromophore. [17][18][19][20][21][22][23][24] For example, in a recent perspective Orrego-Hernández et al. 24 highlight seven examples of NBD systems with absorption onsets of 362-466 nm, but energy densities of 0.10-0.56 MJ/kg -compared to 0.97 MJ/kg for unsubstituted NBD. Furthermore, four of the seven molecules have half-lives of a week or less, which make them unsuitable for long-term energy storage.
In this study we combine quantum chemical calculations, machine learning, and a genetic algorithm to search chemical space of roughly 10 26 NBD/QC derivatives defined in Figure 2 for optimal MOST candidates. Our results suggest that the largest energy storage value that can be obtained with these substituents, while maintaining a long half-life and absorption in the visible spectrum, is around 0.55 MJ/kg.   Figure 2: A representation of the chemical space investigated using machine learning, showing that the NBD/QC motif can be functionalized on four positions with three different groups ("a", "b", and "c"). The group "a" includes electron withdrawing groups (EWGs) and electron donating groups (EDGs) as well as hydrogen, and accounts for a total of 18 different substituents. There are roughly 18 + 2( 1 2 18 5 ) = 1.9M substituents and 1 4 (1.9M) 4 = 10 25 different NBD systems.

Computational Methodology
We show that semiempirical tight-binding methods (SQM) can be used used to identify molecules with high energy densities and thermal back reaction (TBR) barriers as well as suitable absorption maxima by benchmarking against DFT calculations and experiment. We use the SQM methods to construct a dataset that provided even coverage of the chemical space of interest and train graph convolutional neural networks to predict energy densities, TBR barriers, and absorption spectra and then use the models together with a genetic algorithm to search chemical space for promising MOST candidates.
Semiempirical tight-binding calculations ∆H • storage and ∆G •, ‡ T BR are approximated as the differences in electronic energy (∆E storage and ∆E ‡ T BR ) computed using GFN2-xTB. 25 For each NBD structure the QC structure is automatically generated using RDKit. 26 5 + 5n rot random conformations (where n rot is the number of rotatable bonds in the molecule) are then generated for each structure using RD-Kit and optimized with GFN2-xTB. Optimizations that result in discrepancies between the input and output connectivity are discarded. The lowest energy conformers of NBD and QC are used to compute ∆E storage .
To compute the energy barrier of the thermal back reaction a concerted adiabatic scan is performed for the two breaking CC single bonds of QC-to-NBD, which are constrained to 20 values from 1.5 Åout to 2.2 Å, starting from the QC structure with the lowest energy. The highest energy structure and energy is used as an estimate for the transition state. The absorption spectra of NBD are computed by the sTDA-xTB method 27 using the lowest energy GFN2-xTB structure.
The entire process is automated and requites only SMILES strings of NBD derivatives as input (see flowchart in supporting information).

The machine learning model
Three undirected graph convolutional networks (GCNs) 28 are trained to reproduce the energy storage, absorption, and TBR barrier. The GCN model is essentially that implemented in DeepChem 29 and written in Python 3.6.8 using PyTorch version 1.2.0 30 and PyTorch Geometric version 1.3.2. 31 The feature vectors of the nodes describe the atom type, number of directly bonded neighbors, number of hydrogen atoms attached to the atom, formal charge, hybridization, and aromaticity, while the feature vectors of the edges include bond orders as well as information about conjugation and presence in a ring system. The GCN uses two graph convolutional layers, with 128 channels in the first layer and 64 in the second. After the final graph convolutional layer, a global max-pooling layer creates a representation with 64 values, which is fed to the feed-forward neural network. The feed-forward network has two hidden layers with 64 nodes in each layer. The final layer outputs a single value i.e. the target value predicted by the network.

The genetic algorithm
In a previous study 32 we have shown that chemical space of MOST candidates can be searched efficiently by genetic algorithms and we apply a similar approach here. The chemical space shown in Figure 2 is encoded by a gene like that shown in Figure 3. Each of the four positions is represented by a label ("A", "B", and "C") followed by five integers ranging from 0 to 17. The labels represent the three types of ligands shown in Figure 2a, b, and c, respectively, while the five integers represent the five different positions on the phenyl rings shown in Figure 2b and c (if the label is "A", only the first digit is read). The GA search generates an initial population of 300 genes corresponding to random singly substituted NBD/QC systems (phenyl and phenylacetylene groups are also singly substituted and random mutations are applied if duplication occurs) after which a new generation of same size is created from crossovers between two parent genes. In the crossover, between one and three ligands from the first parent are randomly selected and combined with the remaining ligands from the second parent to produce a child. Furthermore, random mutations of the children occurs with a mutation rate of 25 %. These mutations allows one of the ligands to be modified by either changing the ligand type ("A", "B", or "C" from Figure 2) or one of the attached "A" ligands. The different parents are selected with a probability that is proportional to their score (roulette selection). This score is computed as sum of three thresholded-linear functions 33 ranging from 0 to 1 with thresholds of 0.6 MJ/kg, {375, 400, 450, 525} nm, and 250 kJ/mol for the energy density, absorption, and TBR barrier, respectively, predicted by the GCN models. To favour a redshift of the absorption spectrum, we multiplied the absorption score by two resulting in a maximum GA score of four. Hereafter, 300 of the worst scoring genes are discarded, which results in a new population of 300 genes. Each GA search runs for 300 generations and the final output is the 50 highest scoring molecules of the final population. We select four different absorption thresholds resulting in four runs of a thousand GA searches to also promote higher energy densities and TBR barriers.

Benchmarking the tight binding calculations
We benchmark the tight binding calculations against DFT on 31 different NBD/QC systems with promising properties collected from the literature. [17][18][19][20]34 Based on a recent benchmark study by some of us, 35 we use M06-2X/6-311+G(d) for the energy density, and the same method is also used for benchmarking the TBR barrier estimates.
Figures S1 and S2 compares the GFN2-xTB storage energy and barriers of the back reaction to the corresponding M06-2X/6-311+G(d) values. While GFN2-xTB significantly underestimates the storage energy and overestimates the barrier heights both properties are strongly correlated with the DFT results, with R 2 values of 0.81 and 0.70 and Pearson's R values of 0.90 and 0.84, respectively. GFN2-xTB is thus capable of identifying molecules with high energy densities and low barrier heights for further study using DFT.
The sTDA-xTB UV-Vis spectrum of a charge-tagged NBD/QC carboxylate is benchmarked against the experimental result by Ugo et al. 34 and several DFT calculations (Figure S3). The relevant peak at 315 nm is blueshifted by just 10 nm by sTDA-xTB. In fact, sTDA-xTB performs better than 13 out of the 16 DFT functionals tested in a benchmark study, 35 while requiring only a few seconds instead of several hours.

Training the machine learning models
For the dataset, different singly and doubly substituted NBD/QC systems are selected to get an even coverage of chemical space where every small ligand (Figure 2a) is represented at every position on the phenyl rings and NBD scaffold a roughly equal number of times. The singly substituted NBD/QC systems include all directly attached EWGs and EDGs as well as phenyl and phenylacetylene groups with all combinations of up to three EWGs and EDGs in the ortho-, meta-, and para-position (resulting in a total of 18 3 · 2 + 18 = 11, 682 singly substituted systems). The doubly substituted NBD/QC systems (substituted in position 1 and 2, 1 and 3, or 1 and 4, see Figure 2) included all combinations of the EWGs, EDGs, and single ortho-, meta-, or para-substituted phenyl and phenylacetylene groups (resulting in a total of ((17·3+1)·2+17)+2−1 2 · 3 = 22, 143 double substituted systems). Hence, the dataset consisted of 33,825 unique NBD/QC systems. Of these, 119 of the GFN2-xTB calculations failed and are omitted from the dataset. This results in a total of 33,706 NBD/QC systems, which are split 80/20 for training and test set, respectively. A 5-fold cross-validation is used to train the different GCNs for 100 epochs with a batch size of 512 and using the Adam optimizer 36 with a learning rate of 0.01 on a MSE loss (see learning curves in supporting information). To prevent overfitting, the trained GCN with the overall best validation loss for each property are saved.

Identifying promising candidates using a genetic algorithm
The 50 top-scoring molecules from the final populations of 4,000 GA searches are selected, which results in 1,234 unique NBD/QC systems. The energy storage, TBR barrier height, and absorption spectra of these molecules are then recomputed using GFN2-xTB and sTDA-xTB, respectively. Then, the 230 systems with energy densities, absorption spectra, and TBR barriers above 0.2 MJ/kg, 350 nm, and 160 kJ/mol, respectively, are reoptimised at the M06-2X/6-311+G(d) level of theory (using the lowest energy GFN2-xTB structures as starting point) followed by vibrational analysis to ensure correct convergence. TDDFT calculations on all 230 molecules reveal that 116 molecules have absorption maxima above 350 nm. The energy densities for these 116 molecules are then computed and the 38 molecules with energy densities higher than 0.4 MJ/kg are selected for barrier computations. The TS search fail for four of the 38 molecules while the TBR barrier is higher than 150 kJ/mol for 15 of the 34 molecules (shown in Table 1). Five of the molecules have energy densities greater than 0.45 MJ/kg, which is similar to the highest values observed so far for NBD-dimer systems and considerably higher than any previously observed for NBD-monomers. 24 The absorption spectra of these five molecules are shown in Figure 5. However, most of the molecules have one or more amine and hydroxy groups directly attached to the NBD scaffold and may undergo keto-enol or imine-enamine tautomerisation, which may adversely impact the properties. The main conclusion of this study is that the largest energy storage value that can be obtained for a single NBD/QC moiety with these substituents, while maintaining a long half-life and absorption in the visible spectrum, is around 0.55 MJ/kg.  (1)(2)(3)(4)(5). Furthermore, the UV-Vis spectrum of the parent NBD isomer is given as a reference (NBD). All the spectra are obtained using the time-dependent analog of M06-2X/6-311+G(d) and simulated using Eq. S1 from the supporting information.

Conclusions
We present a computational methodology for the screening of a chemical space of 10 25 substitutet NBD molecules for promising MOST systems with high energy density and TBR barrier that absorb in the visible part of the solar spectrum. We show that semiempirical tight-binding methods (SQM) can be used used to identify molecules with high energy densities and thermal back reaction (TBR) barriers as well as suitable absorption maxima by benchmarking against DFT calculations and experiment. We use the SQM methods to construct a dataset of nearly 34,000 molecules that provided even coverage of the chemical space of interest and train graph convolutional networks to predict energy densities, TBR barriers, and absorption spectra and then use the models together with a genetic algorithm to search chemical space for promising MOST candidates.
The 50 top-scoring molecules from the final populations of 4,000 GA searches are selected, which results in 1,234 unique NBD/QC systems. Then M06-2X/6-311+G(d) TDDFT caculations are performed on the 230 systems with SQM-predicted energy densities, absorption spectra, and TBR barriers above 0.2 MJ/kg, 350 nm, and 160 kJ/mol, respectively. TDDFT calculations on all 230 molecules reveal that 116 molecules have absorption maxima above 350 nm. The energy densities for the 116 with absorption maxima above 350 nm are then computed and the 38 molecules with energy densities higher than 0.4 MJ/kg are selected for barrier computations. The final results are 15 molecules with a TBR barrier higher than 150 kJ/mol.
Five of the molecules have energy densities greater than 0.45 MJ/kg and the main conclusion of this study is that the largest energy density that can be obtained for a single NBD/QC moiety with these substituents, while maintaining a long half-life and absorption in the visible spectrum, is around 0.55 MJ/kg.

Data Availability
Additional figures and tables can be found in supporting information. The code and data are available at https://sid.erda.dk/sharelink/DeRV97z1Nz

Estimation of chemical space size
There are 18 different small substituents (including H) as shown in Figure 2a. There are five possible substitution sites on the phenyl ring (Figure 2b-c), so there is on the order of 18 5 different substituted phenyl rings but since the ortho-and para-positions are symmetrically equivalent this number should be reduced by roughly a factor of two. Thus there are roughly 2( 1 2 18 5 ) phenyl substituents (with and without the ethyl linker) plus 18 non-phenyl substituents (including H). These 1.9M substituents can be placed on four different, but symmetrically equivalent, sites in the NBD scaffold, resulting in roughly 1 2 (1.9M ) 4 different NBD systems. Table S1: Information about the systems used in Figure S1 and Figure S2 including references to experimental studies examining the systems.  Procedure for obtaining the learning curves

Reactant
• Perform a random shuffle of the training set containing 26,965 NBD/QC systems.
• Use 5-fold cross-validation to train the GCN on the selected subset for 100 epochs with a batch size of 512 and using the Adam optimizer with a learning rate of 0.01 on a MSE loss.
• Test the performance on the validation set for every epoch and save the best model.
• Use the best model to obtain the mean absolute error (MAE) on the training, validation and test sets.
• Continue to the next fold and report the mean and standard deviation of the 5 models obtained in the cross-validation loop.
• Continue to the next subset.
The code for training the GCNs and creating the learning curves are available at https://sid.erda.dk/sharelink/DeRV97z1Nz   Figure S11: Comparison of the results for the investigated properties obtained by graph convolutional network (GCN), GFN2-xTB, and M06-2X/6-311+G(d). The top panels (a, b, and c) show the GCN vs xTB results for results of the 1,234 unique NBD/QC systems obtained from saving the 50 top-scoring molecules from the final populations of 4,000 GA searches. The bottom panels (d, e, and f) show the results of the NBD/QC systems with SQM-predicted energy densities, absorption spectra, and TBR barriers above 0.2 MJ/kg, 350 nm, and 160 kJ/mol, respectively. However, systems without a confirmed transition state structure at the DFT level of theory are omitted in the bottom panels, which limits the number of presented NBD/QC systems from 230 to 147. The green dots corresponds to the 15 NBD/QC systems shown in Table 1 and the orange line represents the line of equality.

Simulating UV-Vis Spectra
In this section, we present the expression used to simulate the UV-Vis spectra. The expression is derived from the assumption of Gaussian band shapes by applying Gaussian functions to convolute the calculated oscillator strengths. Thus, the UV-Vis spectra can be plotted as the extinction coefficient, , vs. the wavelength, λ, using the following equation.
where f i is the calculated oscillator strength, λ i is the corresponding wavelength in nm, λ is an independent variable defining the simulated spectrum, and σ is the standard deviation also known as the full width at half maximum of the Gaussian band (in these simulations σ = 0.4 eV = 0.4 · 8065.54 cm −1 = 3226.22 cm −1 ). Furthermore, the constant k is given by k = N A · e 2 2 · m e · c 2 · 0 · ln(10) · ln (2) π = 2.1751 · 10 8 L mol · cm 2 (S2) where N A is Avogadro's constant, c is the speed of light, e is the elementary charge, m e is the mass of an electron, and 0 is the vacuum permittivity.