UV-Visible Absorption Spectra of Solvated Molecules by Quantum Chemical Machine Learning

The calculation of the UV-visible absorption spectra follows a multiscale approach which combines FPMD and excited state calculation. This approach accurately predict line shapes as it includes both temperature and solvation effects, within the limits of the method used to compute electronic excitations.Ge et al. 2015; Timrov et al. 2015; Timrov et al. 2016; Zuehlsdorff and Isborn 2019 In this work, density functional theory (DFT) based FPMD simulations are performed using the mixed basis-set Quickstep approach, implemented in the CP2K packageVandeVondele et al. 2005; Hutter et al. 2014. We use the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) for the exchange and correlation functionalPerdew et al. 1996 with D3 van der Waals corrections.Grimme et al. 2010 Valence Kohn-Sham orbitals are expanded in real space on a double-$\zeta$ localized valence polarizable (DZVP) basis setFeller 1996, the electron density in reciprocal space is expanded in plane waves up to a cutoff energy of 300 Ry, and core states are treated implicitly using Geodecker-Teter-Hutter norm-conserving pseudopotentials.Goedecker et al. 1996 We run MD simulations of each molecule in aqueous environment in cubic periodic cells with 128 water molecules. FPMD runs were initialized from classical MD simulations in the constant pressure canonical ensemble at ambient conditions, using the TIP4P/Ice model for water and the generalized Amber force field (GAFF) for the organic molecules.Wang et al. 2004; Abascal et al. 2005 These systems were then equilibrated for 10 ps by FPMD at constant volume at 300 K. The frames used for absorption spectra calculations were extracted from 100 ps long production runs in the constant-volume canonical ensemble, enforced using the stochastic rescaling algorithm with a coupling constant of $\tau=10$ ps.Bussi et al. 2007 The equations of motion were integrated with a timestep of 0.5 fs.

For each molecular geometry sampled from the FPMD simulations, VEEs were computed by “domain-based local pair natural orbital similarity transformed equation of motion coupled cluster singles and doubles” (DLPNO-STEOM-CCSD): an efficient coupled-cluster approach with $O(N{{}^{4}})$ scaling.Berraud-Pache et al. 2019 For the sake of efficiency, the explicit solvent was replaced by a conductor-like polarizable continuum model(CPCM).Pascual-Ahuir and Silla 1990; Silla et al. 1991; Pascual-ahuir et al. 1994 All electron DLPNO-STEOM-CCSD calculations were performed with a triple-$\zeta$ valence polarized Karlsruhe basis set (def2-TZVP)Weigend and Ahlrichs 2005 using ORCA 4.2.1.Neese et al. 2020 For the excited state calculations, the resolution of identity approximation for both Coulomb and exchange integrals was employed to speed up the self-consistent calculation. For each molecule, the lowest-energy absorption band of the UV-visible spectrum is simulated by summing Gaussian functions with a width of 0.027 eV centered on the calculated VEE for 200 statistically independent configurations chosen from the FPMD trajectory. This approach corresponds to the “ensemble method", which provides a reasonable approximation of line widths but does not capture vibronic effects arising from nuclear dynamics.Zuehlsdorff et al. 2019 Furthermore, we do not re-weight the contributions to the spectrum using the oscillator strength obtained from QC calculations, as we have verified that for each molecule this quantity does not vary significantly over the ensemble of structures that we used to obtain the absorption spectra.

The goal of our ML model is to predict VEEs as a function of the molecular coordinates, so to spare the computational burden of the QC calculations. To this scope, we represent the molecular configuration in BC.Bartók et al. 2010 While originally BC was proposed as a descriptor to approximate the Born-Oppenheimer potential energy surface of single element systemsBartók et al. 2010, it has been applied to develop ML potential of multi-component systems and predict material properties such as elastic constants, bulk modulus as well as vibrational free energies and entropies of solidsLegrain et al. 2017; Li et al. 2018; Cusentino et al. 2020.Compared to other atomic environment descriptors such as smooth overlap of atomic positions (SOAP)Bartók et al. 2013; Bartók et al. 2017; Wilkins et al. 2019 and atom-centered symmetry functions (ACSF) Behler 2011; Bartók et al. 2013, BC is more keen at describing the nuances of atomic environments, as it is projected to a more complete set of basis functions with higher dimensions Pozdnyakov et al. 2020; Zuo et al. 2020. Additionally, in the development of spectral neighbor analysis potentials, BC was formulated to retain a linear relation with the target property.Thompson et al. 2015; Cusentino et al. 2020 This development ensures the resulting ML models to achieve robust performance using only a moderate amount of training data. Such a trade-off between model complexity and size of the training data is the most critical factor for us to bridge BC with linear regression.

Hereafter, we summarize the key formulations of BC to supplement the following parameterization of our ML model. As reported in the original works,Thompson et al. 2015 the atomic environment is expressed as the weighted atomic density ($\rho_{i}(r)$):

where $r_{ii^{{}^{\prime}}}$ is the interatomic distance between the central atom $i$ and the neighboring atom $i^{{}^{\prime}}$, $f_{cut}$ is the cutoff function to smoothly decay the neighboring atomic density to 0 at the pair-wised cutoff radius ($R_{cut,ii^{{}^{\prime}}}$). $R_{cut,ii^{{}^{\prime}}}$ is computed by summing over cutoff radii pairs between central and neighboring atoms. The dimensionless weighting factor ($\omega_{i^{{}^{\prime}}}$) is used to differentiate the neighboring atoms. After expressed as a sum of $\delta$ functions, $\rho_{i}(r)$ is expanded in hyperspherical harmonics ($U_{j}(\theta,\phi,\theta_{0})$):

In equation (2), $u_{j}$ is the Fourier expansion coefficient given by the inner product between $\rho_{i}(r)$ and $U_{j}(\theta,\phi,\theta_{0})$. In this implementation, $j$, $j_{1}$ and $j_{2}$ are truncated so that $j,j_{1},j_{2}\leq j_{max}$ to ensure a finite spatial resolution of the weighted atomic density. $2j_{max}$, the even integrable of $j_{max}$, represents a hyperparameter to dictate the number of BC used to fit the ML model. With $u_{j}$, the bispectrum components ($B_{j_{1},j_{2},j}$) can be defined as:

where $\otimes_{j_{1}j_{2}j}$ represents a coupling product analogous to angular momentum coupling of spherical harmonics. The $\frac{1}{2j+1}$ prefactor ensures $B_{j,j_{1},j_{2}}$ invariant under permutation of the atom indices. In this work, BC were computed using the FitSNAP package.git

The goal of our ML model is to predict the first VEE from the molecular geometry with quantum chemical accuracy. For this purpose we train a linear model on the DLPNO-STEOM-CCSD/def2-TZVP excited state calculations described above using BC as descriptors of the molecular configurations. As described in the previous section and in former works,Bartók et al. 2010 the BC descriptor consists of projecting the local atomic environment on hyperspherical harmonics for each atomic species. As the BC descriptor is high-dimensional, we apply the least absolute shrinkage and selection operator (LASSO)Tibshirani 2011; Pedregosa et al. 2011, a linear regression model with norm one ($\ell_{1}$) regularization. Training LASSO consists of minimizing the loss function with respect to the set of coefficients $\beta$:

where $N$ is the number of molecular configurations extracted from the FPMD simulations, $\alpha$ is the regularization parameter, and $\varepsilon_{QC}$ are the first excitation energies obtained from the QC calculations. The performance of the ML model is assessed through two statistical metrics: mean absolute error (MAE) and mean signed error (MSE). The latter, defined as:

is useful to spot the occurrence of systematic errors in the predicted VEEs. After training the LASSO model, we compute the UV-visible absorption spectra by estimating $\varepsilon_{ML}$ for 5000 statistically independent configurations from the FPMD trajectories. The final ML spectrum consists of the envelop of Gaussian functions centered on the $\varepsilon_{ML}$ with a width of 0.014 eV. Whereas our ML model is developed in close analogy with the approach used to construct SNAP interatomic potentials,Cusentino et al. 2020 the use of LASSO marks the main difference between SNAP and our excited state model. In a forthcoming section, we discuss extensively the importance of using a norm $\ell_{1}$ selection operator to shrink the space of ML parameters and how it enhances the predictivity of the ML models beyond the set of molecules on which it is trained.

A further advantage of using atomic descriptors and the LASSO model is that we can compute the relative contribution of each atom or group of atoms to $\varepsilon_{ML}$, through a linear decomposition analysis, such that:

where $\varepsilon_{0}$ is the intercept of the LASSO model and $\varepsilon_{\gamma}$ represents prediction contributions from the atom type $\gamma$ for atom $i$. As $\varepsilon_{ML}$ is a scalar quantity and the LASSO model is linear, $\varepsilon_{ML}$ can be partitioned with respect to different functional groups:

where $\varepsilon_{NH_{2}}$, $\varepsilon_{OH}$ and $\varepsilon_{OCH_{3}}$ are the prediction contributions from the amine, hydroxyl and methoxy groups. $\varepsilon_{C_{6}H_{n}}$ is the prediction contributions from carbon and hydrogen atoms within a phenyl ring ($n=5$ for phenol, anisole, aniline and $n=4$ for 1,2-methoxyphenol, 1,4-aminophenol, 1,4-anisdine, 1,4-benzenediol and dimethoxybenzene isomers). $\%_{group}$ is used to express the prediction contributions by percentage.

First, to show the baseline performance of our ML approach, we fitted a model using the full data set of ten molecules with 200 VEE per molecule. Then, to test model generalizability, we developed a 7-molecule model by leaving out of the QC calculations of 1,2-methoxyphenol, 1,4-aminophenol and 1,4-dimethoxybenzene from the training set. Before optimizing the 7- and 10-molecule models, we set the parameters of the BC descriptor, specifically the weights ($\omega_{atom}$) and the cutoff radius for each atomic species ($R_{cut,atom}$). To optimize $R_{cut,atom}$, we fixed $R_{cut,H}=1.2$ Å, and performed a grid search on $R_{cut}$ for each heavy atom between 2.6 and 3.2 Å. We simplified the choice of $R_{cut,atom}$ parameter by setting $R_{cut,C}=R_{cut,N}=R_{cut}$ and $R_{cut,O}=1.05R_{cut}$. Since carbon is the major building block of these aromatic molecules, we set the weighting factor of carbon ($\omega_{C}$) to unity. For the remaining weighting factors ($\omega_{H},\omega_{O}\ \&\ \omega_{N}$), a constraint of $\omega_{H}<\omega_{N}\leq\omega_{O}$ was imposed. After determining $\omega_{atom}$ and $R_{cut,atom}$, we chose $2j_{max}$ based on the number of available training samples. At $2j_{max}=18$, the total number of BC almost equals the number of training samples ($N_{train}$) used to develop the 7-molecule model. Thus, to retain a similar number of input features and training samples, we chose $2j_{max}=18$. The hyperparameters $R_{cut}$ and $\omega_{atom}$ are optimized for the 7-molecule models and their value is used also for the 10-molecule model. In this process no information is used from the three molecules left out of the training set. Table 1 summarizes the hyperparameters used to compute BC.

The $\ell_{1}$ regularization in LASSO is designed to achieve feature eliminationPedregosa et al. 2011, so to prevent overfitting. From Table 2, we see that for both ML models, the number of non-zero features is less than 25 $\%$ of the total number of input features ($N_{features}$). The 7- and 10-molecule models retain similar fractions of non-zero features with respect to the size of the training sets so that the model performance can be compared. A 10-fold cross validation was applied to both models to avoid bias from a specific split of training and testing data sets.

As shown in Table 2, the testing MAE of our 10-molecule model, averaged from the 10-fold cross validation, is 23.05 $\pm$ 3.92 meV. A similar MAE, 23.33 $\pm$ 4.13 meV, is also observed from the 7-molecule model. Besides MAE, we also introduced MSE to gauge if our ML models systematically overestimate or underestimate excitation energies. As shown in Figure 2, most of the molecules have MSE of less than 5 meV, which implies that no noticeable overestimation or underestimation of excitation energies occurs, except for 1,4-anisidine, for which the predicted excitation energies have MAE and MSE of 25.91 and 11.57 meV.Figure 3 illustrates the testing performance of the 7-molecule model. It can be noticed that both 10 and 7-molecule models have similar testing MAE and MSE for molecules in the training set. It is worth highlighting, from the 10-molecule model, that the predicted excitation energies for every single molecule have MAE far below the intrinsic error of the underlying EOM-CCSD method (70 meV)Dutta et al. 2018. Even for the excitation energies interpolated from the 7-molecule model, the MAE are merely around half of the 70 meV. Hence, no significant error is introduced by both ML models.

Figure 4 compares the lowest-energy UV-visible absorption band computed by the multiscale quantum-chemical approach described in the methods section to experimental measurements for each molecule in aqueous solution.Shiobara et al. 2003; Stalin et al. 2005; Malongwe et al. 2016; Corrochano et al. 2017; Hullar et al. 2020; Chen et al. 2020; Hullar et al. 2021 The QC model is in excellent agreement with experiments for 8 out of 10 molecules, with differences in the center of the peak within less than 0.1 eV and nearly overlapping low-energy tails. Experimental bands tend to be broader on the high energy side, possibly because they encompass the tails of higher excited states. 1,4-anisidine and 1,3-dimethoxybenzene make an exception, with differences in peak positions of 0.18 and 0.15 eV. A possible source of discrepancy is that the GGA functional used in FPMD simulations leads to systematic errors in the geometry of the molecules. The differences between the gas-phase excitation energies of 1,3-dimethoxybenzene computed for geometries optimized with the PBE and PBE0Adamo and Barone 1999 functionals support this hypothesis (Table S1). However, the cost of running FPMD with PBE0 would be excessive for the purpose of this study. Choosing a semilocal GGA functional strikes the desired balance between the accuracy of computing spectra and a reasonable computational cost. Additionally, the experimental conditions may be slightly different from those in the models. The absorption spectrum of 1,4-anisidine was measured at pH$=6.0$, while in our FPMD simulations molecules solvated in aqueous solution at neutral pH, and this difference may cause a shift in the UV-visible absorption.Shiobara et al. 2003; Zheng et al. 2018

The QC spectra were obtained by averaging only a few hundred frames for each molecule. Hence, the calculated spectra for some of these molecules, especially 1,4-dimethoxybenzene, exhibit jagged line shapes (Figure S1). As shown in previous works, ML may be employed to obtain smoother and more refined line shapes.Bononi et al. 2020; Xue et al. 2020 The same effect is achieved here, as illustrated in Figure 4. The ML spectra are obtained by averaging over 5,000 frames for each molecule, which guarantees smooth line shapes and converged statistical sampling. At this point, convergence is limited only by the capability of FPMD simulations to sample the configurational space of the molecules in aqueous solution. The main features of the ML spectra obtained with this model are almost indistinguishable from the QC references as ML predictions have the same level of accuracy as the training data. From Figures 2 and 4, we can conclude that a unified ML model is able to accurately predict the excitation energies and thus allowing us to compute the spectra for thousands of molecular configurations at nearly no additional computational cost.

Although useful to save computational time to compute the lowest energy band of the absorption spectra, the application of the ML model described so far is limited to the molecules comprised in the training set. Hereafter, we explore whether our approach can be used to predict the absorption spectra of molecules that were not included in the training set, provided that these molecules are similar to those used to construct the ML model. To this scope, we employ the LASSO model fitted with the excited state calculations of 7 molecules to interpolate the excitation energies and the corresponding UV-Visible absorption spectra for 1,2-methoxyphenol, 1,4-aminophenol and 1,4-dimethoxybenzene. These three molecules consist of a phenyl ring with side groups of the same kind (methoxy, hydroxyl and amine) as the molecules included in the training set, but combined to form different isomers or different molecules altogether. Table 3 and Figure 5 show that the peak positions and the absorption line shapes for the seven molecules in the training set remain the same as those obtained with the 10-molecule ML model. Figure 3 indicates an increase of both MAE and MSE for the three molecules excluded from the training set. MAE, in particular, raises beyond 20 meV for 1,4-aminophenol and 1,2-methoxyphenol. A larger error for these molecules is expected, as the side-group arrangements (e.g., ortho positioning of hydroxyl and methoxyl group, para positioning of amine and hydroxyl group and of two methoxyl groups) are not explicitly known by the ML model and their effects on the excitation energies are inferred from the other seven molecules. Nevertheless, a close comparison between Figures 4 and 5 suggests that the effect of these errors on the predicted UV-visible absorption spectra is small. In fact, the model-generalized spectra remain very similar to those obtained with the 10-molecules ML model, are in excellent agreement with the reference QC spectra. The difference between $\varepsilon_{max,experiment}$ and $\varepsilon_{max,ML}$ is 0.084 eV, 0.077 eV and 0.038 eV for 1,2-methoxyphenol, 1,4-aminophenol and 1,4-dimethoxybenzene. Considering that the 7-molecule model is developed without using any QC calculations of the three excluded molecules, less than 0.1 eV differences are remarkable.

These results suggest that the 7-molecule model can properly interpolate the electronic excitations of the molecules not included in the training set without further tuning the hyperparameters of the BC descriptors and provide UV-visible absorption spectra in good agreement with experiments. However, while the low-energy tail of the absorption band is reproduced very well, our approach fails to reproduce the skewed experimental line shapes. To improve the line shape predictions we have applied the dynamics-based third-order cumulant scheme, using the fluctuations of the VEE estimated by the ML model along the FPMD trajectory (Fig. 5 solid line).Mukamel 1985; Zuehlsdorff et al. 2019; Chen et al. 2020 Whereas we get a minor systematic improvement in the prediction of the line shapes with this approach, the theory still underestimates the intensity of the high-energy part of the absorption band. This may be due to the fact that nuclei are treated as classical particles in FPMD, thus neglecting nuclear quantum effects.

This result shows that local geometric descriptors are sufficient to predict very accurately excitation energies of a group of molecules with similar features. This is a promising starting point for future work on more complex molecules. However, to extend the transferability of this approach to other classes of molecules, given the non-local nature of electronic excitations in large molecules, it may turn out necessary to supplement our approach by including richer physically-based descriptors, e.g. electronic orbitals, and to compute excitation energies as the eigenvalues of ML effective Hamiltonians.Zhang et al. 2020; Westermayr and Maurer 2021; Nigam et al. 2022

Hereafter, we analyze the features that make the ML model developed in this work accurate and predictive. We first examine the importance of $\ell_{1}$ regularization. To this aim, we built another 7-molecule model using ordinary least square (OLS)Puntanen 2013 as the underlying ML model. BC for this model were computed using the hyperparameters as summarized in Table S2. The $2j_{max}$ was chosen to be 8 so that both 7-molecule models have similar $N_{features\neq 0}$. Table S2 shows that both 7-molecule models have very similar overall MAE. The 7-molecule + OLS model even achieves lower standard deviation (std) in overall MAE than the 7-molecule model. Surprisingly, two 7-molecule models show striking differences in MAE with respect to each molecule. From Figure 3 and S2, one can see that the 7-molecule model has lower MAE for 5 molecules in the training set. The 7-molecule model also records 8.93, 1.37 and 13.63 meV lower in MAE than the 7-molecule + OLS model when interpolating the excitation energies for 1,2-methoxyphenol, 1,4-aminophenol and 1,4-dimethoxybenzene. The std of MAE from the 7-molecule model (0.724 meV) is only about $12.5\%$ of the std from the 7-molecule + OLS model (5.76 meV) for the three model-generalized molecules. The 7-molecule model is fitted against 225 BC selected from 1544 BC generated with $2j_{max}$ up to 18. As shown in Figure S3, $2j_{max}=18$ is a sufficiently high order to generate BC with optimal MAE. Meanwhile, feature elimination schemes, such as $\ell_{1}$ regularization, help capture the essential features of the atomic environments and ensure little to no performance loss as the final model is fitted on a carefully-chosen subset of the initial feature space.Imbalzano et al. 2018. Therefore, these 225 BC are composed of comprehensive descriptions of the atomic environments at resolution as high as $2j_{max}=18$. However, as no feature elimination is imposed in the 7-molecule + OLS model, its feature space can only be confined at lower order of $2j_{max}$ to prevent overfitting. In particular, each of 224 BC generated with $2j_{max}=8$ is used to fit the 7-molecule + OLS model. Thus, despite almost identical $N_{features\neq 0}$, the difference in the detailed description of the atomic environment eventually leads to a drastic difference in the capability for the two 7-molecule models. Whereas a more standard $\ell_{2}$ regularization may be used, the latter is prone to overfitting when the ratio between the number of samples and features approaches unity. It is therefore very advantageous to impose the $\ell_{1}$ regularization as it not only prevents overfitting, but it also enables the ML model to have a sufficiently large initial feature space, which is critical for the model to be generalized to systems not included in the training set.

To trace how the ML models explore the structural similarity among these aromatic molecules, we performed the linear decomposition analysis as formulated above in the methods section. Figures 6 and S4 show that, for $\varepsilon_{ML}$ from both ML models, contributions from the aromatic ring are predominant. This trend is common to all molecules. This result suggests that the observed generalizability across the whole family of aromatic molecules lays primarily on geometric variations of the aromatic ring during the FPMD simulations, and that the primary effects of solvation are local and can be tracked down to configurational changes in the solvated molecule.P. et al. 2017 To interpret this result from a quantum mechanical standpoint, we computed the Natural Transition Orbitals (NTOs)Martin 2003; Badaeva and Tretiak 2008; Lu and Chen 2012 of all the aromatic molecules using their equilibrium structures in aqueous solution. As shown in Figures S5, the dominating NTOs for all these aromatic molecules exhibit the characteristics of a $\pi\rightarrow\pi^{*}$ transition. Such excited state characters are mostly contributed by the aromatic ring, along with a secondary participation of lone pairs on functional groups. Whereas relative contributions of the hydroxyl and amine groups from the NTOs are noticeably larger than those illustrated in the linear decomposition analysis, NTOs confirm the predominant relevance of the carbon atoms in the aromatic ring. Quantitative differences in the relative contributions of the functional groups indicate that the geometrical interpretation of the excitations from ML is qualitative.