Evaluating Fast Methods for Static Polarizabilities on Extended Conjugated Oligomers

Given the importance of accurate polarizability calculations to many chemical applications, coupled with the need for e ciency when calculating the properties of sets of molecules or large oligomers, we present a benchmark study examining possible calculation methods for polarizable materials. We first investigate the accuracy of highly-e cient semi-empirical tight-binding method GFN2-xTB, and connected D4 dispersion model, comparing its predicted additive polarizabilities to !B97XD results for a subset of PubChemQC and a compiled benchmark set of molecules spanning polarizabilities from approximately 3Å3 to 600 Å3, with a few compounds in the range of approximately 1200 Å3-1400 Å3. Although we find GFN2, and thus D4, to have large errors with polarizability calculations on large oligomers, it would appear a quadratic correction can remedy this. We also compare the accuracy of DFT polarizability calculations run using basis sets of varying size and level of augmentation, determining that a non-augmented basis set may be used for large, highly polarizable species in conjunction with a linear correction factor to achieve accuracy extremely close to that of aug-cc-pVTZ.


Introduction
Polarizability plays a key role in many chemical processes and phenomena, and its accurate calculation is therefore crucial to a variety of applications. Because of its fundamental role in explaining dispersion forces, 1 it is a key component of widely-used dispersion corrections for computational calculations. 2 The importance of accurate electrostatic interactions has led to the development and widespread use of polarizable force field models for studying systems such as biomolecules 3,4 and ionic liquids. 5 Computationally-derived Raman spectra also rely on the calculation of polarizability tensors. 6 Polarizability values are also necessary to calculate the values of more complex material properties such as refractive index 7 and dielectric constant, 8 often in the context of molecular screening.
Because of its wide utility, polarizability has been the topic of a number of recent computational benchmark studies. Hait and Head-Gordon provided a thorough examination of the performance of a large number of density functional theory (DFT) functionals at the complete basis set (CBS) limit for 132 small molecules. 9 Frediani et. al. used a subset of the Head-Gordon study's molecule set to test the veracity of that study's CBS limit claim using alternative multiwavelet bases in order to reduce potential error. 10 Sauer and co-workers used a benchmark set of 14 heteroaromatic molecules to assess the accuracy of various second-order methods for both static and frequency dependent polarizabilities. 11 Afzal and Hachmann tested various DFT methods to determine the best way to balance accuracy and e ciency for high-throughput non-conjugated polymer screening. 12 While all of these studies provide valuable insight into the relative accuracy of various polarizability methods for di↵erent applications, none of them examine such methods for the high polarizability limit. As shown in our previous work using a genetic algorithm (GA) to search for high dielectric oligomers, there is a need for a polarizability method capable of calculating large polarizabilities (on the scale of 10 2Å3 ) while making e cient use of time and computational resources. 8 As a point of reference, all of the molecular species examined by the previously mentioned studies possess isotropic polarizabilities less than 40Å 3 .
The Hachmann study suggests the viability of extrapolating polymer polarizabilities from oligomers for non-conjugated species, but notes that this proves untenable for species with conjugated backbones due to electron correlation e↵ects in the ⇡-system. 12 Wong and coworkers performed a polarizability benchmark on oligomers of polydiacetylene and polybutatriene which included some hexamer polarizabilities greater than 200Å 3 . 13 It is worth noting that despite augmented basis sets being recommended for accurate polarizability calculations, 14 the basis sets used in this study were not augmented due resource constraints. In another benchmark study from Wong, linear polarizabilities (and hyperpolarizabilities) were found for a range of streptocyanine oligomers. 15 By using CCSD(T)-F12 which enhances basis set convergence, they were able to give a basis set extrapolation of computed polarizabilities from a non-augmented triple zeta basis set. Both studies provide useful data, including the e↵ects of short-range exchange on polarizability and CCSD(T) calculations in the former and an assessment of MP2 quality and the importance of potential lower-energy open-shell states in the latter, however the scope of neither study was intended to examine resource e cient polarizability methods for large oligomers. With the realm of computationally-generated novel materials continuing to grow, data is needed on resource-e cient methods and basis sets to find accurate polarizabilities of largely polarizable molecules.
In this work, we analyze the viability of using a semi-empirical method and smaller basis sets to make large polarizability calculations more tractable. We first test the durability of popular semi-empirical tight-binding method GFN2-xTB (GFN2) through the additive polarizabilities of the D4 dispersion model, on both small organic molecules from Pub-ChemQC 16 for which we expect GFN2/D4 to perform well, and highly polarizable species, for which we expect some degree of notable inaccuracy. Testing on a set of small organic molecules gives us a point of comparison which allows us to precisely identify GFN2/D4's strengths and limitations when it comes to polarizability calculations for larger extended systems. While it is understood that GFN2/D4 is empirically tuned to most accurately calculate the polarizabilities of relatively small molecules, testing this method on larger conjugated systems is of special interest since the same basic method used to compute GFN2 polarizabilities is also central to the D4 dispersion correction used widely with DFT methods on species of a wide variety of sizes and chemical structures. 2,17 We also examine four basis sets of varying size and level of augmentation to determine whether smaller basis sets can be used to calculate large polarizabilities since they minimize issues regarding computation time and potential linear dependence.

Computational methods
Two primary data sets were analyzed in this benchmark. The first set is a randomly chosen subset of approximately 8,400 species from PubChemQC's approximately 3.2 million known small (molecular weights less than 500 a.u.) organic molecules. 16 It was chosen to provide a strong basis of comparison as a set of molecules for which we expected GFN2/D4 to perform well. The second "wide [polarizability] range" set is drawn from previous studies and designed to cover a very wide range of polarizabilities. Drawing from the pool of hexamer structures we had created with our GA, we constructed a set of 54 hexamers with GFN2/D4 predicted polarizability values in the approximate range of 80-280Å 3 . In order to balance out our benchmark set, we also added 19 conjugated oligomers and small molecules with GFN2/D4 predicted polarizability values in the "medium polarizability" range of 4-91 A 3 . Hexamer equilibrium geometries were found using preliminary force-field optimization using OpenBabel 18 with MMFF94 19-23 or UFF 24,25 followed by geometry optimization using GFN2. 17 Equilibrium geometries for the medium polarizability molecules were optimized with ORCA 4.0.0.2 26 using DFT with the B3LYP functional 27-30 and 6-31G(d) basis set. 31,32 All polarizabilities reported in this study are isotropic, meaning they are the average of the diagonal elements of the polarizability tensor. We chose to focus on static, isotropic polarizabilities for this study due to their general applicability to a variety of theoretical and computational fields, including method development, dispersion correction, and polarizable electrostatics. For GFN2 calculations, polarizabilities were calculated using xTB which relies on the D4 method in which polarizabilities are calculated using a weighted sum of precomputed atomic polarizabilities. This was done to estimate how well GFN2 compared generally to DFT, without facing potential resource and/or linear dependence issues likely to arise when using augmented basis sets with large molecules. For the comparison of basis set accuracy, aug-cc-pVTZ 34,35 was selected as the standard of comparison. This basis set was chosen as the standard because di↵use functions are necessary to describe the long-range electron behavior and electron correlation important for polarizability calculations, demonstrated in by Rowley and coworkers' assessment that augmenting basis sets with di↵use functions leads to a substantial increase in accuracy, particularly for polarizability calculations. 14 This study also determined that aug-cc-pVTZ performed better than both a similar non-augmented triple-zeta basis set and an augmented double-zeta basis set. Additionally, Sauer and co-workers found using larger augmented basis sets do not yield substantial accuracy increases for polarizability calculations despite the increased time required. 11

Results and discussion
Due to our interest in finding an e cient method for calculating molecular polarizabilities for novel molecular searches, we sought to test the accuracy of GFN2/D4 on common small molecules. As an initial experiment, we calculated the polarizabilities for approximately 8,400 species from the PubChemQC dataset, using both GFN2/D4 and DFT with the !B97XD functional 36 and the cc-pVTZ basis set. Although this does not allow polarizabilites to be as accurate as when calculated with augmented basis sets, it allowed us to perform thousands of calculations quickly and gave us an initial baseline against which we could compare GFN2/D4. As discussed below, such results can be scaled to augmented basis sets.
Due to the presence of a few outliers, robust linear regression was performed using SciKit learn's Huber regressor method 37,38 with default epsilon value of 1.35 to limit outlier e↵ects.
The y-intercept was also forced to zero, representing the physical reality that a completely non-polarizable molecule should be computed to have zero polarizability by any method.
After performing Huber linear regression (Figure 1), two notable observations were apparent.
While the trendline's slope was very close to one, the values calculated with GFN2/D4 were often substantially lower than those calculated with DFT, with di↵erences as great as over 100Å 3 between the two methods. This error appears to be somewhat systematic, as species with lower polarizabilities generally have smaller di↵erences in calculated values ( Figure 1A) whereas those with high polarizabilities generally have larger di↵erences in calculated values ( Figure 1C). A substantial number of species' values appear as outliers from the regression line, suggesting a level of random error in GFN2/D4 calculations. In order to further explore the performance of GFN2/D4 for large polarizability calculations, we pursued testing a smaller group of molecules with a wider range of polarizabilites.

Highly Polarizable Oligomers
While we are aware that because of its minimal basis set approach GFN2 is best geared toward polarizability calculations for small molecules, we believe it is important to test its performance on a range of species, including conjugated oligomers. We are interested in the viability of GFN2/D4 additive polarizabilities for two primary reasons: first, to test them as resource-inexpensive albeit relatively low accuracy calculations that allow for bulk molecular screening. For example, because of its vast speed-up compared to ab initio methods, we previously used GFN2/D4 to calculate the polarizabilities of novel hexamer structures generated by a genetic algorithm (GA)-driven search for high-dielectric organic conjugated oligomers. 8 For that work, we chose to use GFN2/D4 despite suspecting a degree of inaccuracy at high polarizabilities. The speed-up GFN2/D4 provides over DFT methods was necessary to complete the thousands of hexamer polarizability calculations in a reasonable amount of computation time. We qualified our work by noting that while GFN2/D4 appeared to vastly underestimate large polarizabilities, it appeared to do so in a systematic way, such that the order of the polarizabilities' magnitudes relative to one another was preserved (allowing us to accurately rank molecules by polarizability, as was needed for the GA). The second reason we believe it is both relevant and important to test GFN2's ability to calculate polarizabilities for a wide range of chemical species is that shares the additive polarizability model with the D4 disperson-correction method, which is widely used for a variety of species including those not limited to the types of small molecules for which GFN2 has been calibrated.
To test the integrity of the GFN2/D4 polarizability results for the 73 member "wide range" benchmark set, we ran single point DFT polarizability calculations using the !B97X functional 36 and the cc-pVTZ basis set. 34,35,[39][40][41] This functional was chosen because it allowed us to compare the DFT results for the "wide range" set to the previously computed results for the PubChemQC subset, since the former was performed in ORCA 4.0.0.2, which does not have an option for the exact dispersion-corrected functional used in the latter's calculations in Gaussian 09. This basis set was chosen because it was the largest basis set that we were able to reach convergence with for all species in the "wide range" set in a reasonable amount of computation time.
Performing Huber regression with a fixed intercept at the origin on the GFN2/D4 and !B97X polarizabilities from our "wide range" set ( Figure 2), we observed trends similar to those seen in the PubChemQC study. We again observed smaller di↵erences in polarizabilities less than 50Å 3 , with a mean absolute error of 4.72Å 3 (Figure 2A), and increasingly larger di↵erences as polarizabilities increased, with an MAE of 37.31Å 3 for polarizabilities less than 100Å 3 ( Figure 2B), and an MAE of 145.02Å 3 for the entire set ( Figure 2C). We similarly observed greater variation in polarizability di↵erences as their magnitudes grew, shown by the drastically smaller R 2 value for the full set as compared to the subset with values less than 50Å 3 . Three extreme outliers appeared ( Figure 2C), in which the percent error of the GFN2/D4 calculated value was in excess of 80%. These outliers were run with the same DFT method and basis set using Gaussian (Table S1)

Investigation of Potential GFN2/D4 Improvement Strategies
As shown in the above assessment, GFN2/D4 performs well when calculating relatively small isotropic polarizabilities, in the range <50Å 3 common for most small molecules. For species with larger polarizabilities, especially for long conjugated systems like the hexamers in the "wide range" set, GFN2/D4 appears to systematically underestimate the polarizability. This derives from its use of an atom-additive polarizability model, neglecting the nonlocal polarizability-enhancing e↵ects of electron delocalization. Our results suggest a correction  Figure 2C.
is needed to allow GFN2/D4 to more accurately predict larger polarizabilities.
We began by constructing an additive polarizability model as a baseline comparison for GFN2/D4 performance. For this model, each atom in a molecule was assigned its GAFF atom type, 42 which was used to further assign it to an Alexandria polarizability type and then finally a corresponding atomic polarizability. 4 The molecular polarizability was calculated as the simple sum of these atomic polarizabilities. We found that the additive model performed extremely similarly to GFN2/D4 ( Figure S1), providing computed polarizabilities with high correlation with GFN2/D4 ( Figure S2).
Given GFN2/D4's systematically increasing inaccuracy for large polarizability calcula-tions, we considered additional chemical properties related to electron delocalization that could potentially correct the additive GFN2/D4 polarizability.
Since the polarizability is connected to chemical hardness ⌘ in conceptual DFT, 43,44 which is defined by the HOMO-LUMO gaps, we first examined the GFN2-computed HOMO-LUMO gaps for both the PubChemQC subset and the hexamers from the wide range polarizability set. In principal, the highly polarizable ⇡-conjugated species should have smaller HOMO-LUMO gaps. Unfortunately, there was not a useful correlation between the calculated polarizability and molecular HOMO-LUMO gap ( Figure 4A), likely because GFN2 is not parameterized for HOMO or LUMO eigenvalues to connect with ionization potential or electron a nity.
We then used an empirical descriptor of the geometric size of the largest conjugated ⇡-system. 45,46 While better correlated to GFN2/D4 polarizability than HOMO-LUMO gap, this information was not enough to meaningfully correct large polarizabilities ( Figure 4B). Plotting DFT polarizabilities against GFN2/D4 polarizabilities for the PubChemQC subset and the "wide range" set, we examine the e↵ects of using a polynomial fit. As an aside, because the PubChemQC subset and "wide range" sets were computed at di↵erent times using slightly di↵erent methods, with Gaussian with a dispersion correction and with ORCA without a dispersion correction, respectively, the functional has been labeled !b97X(D) here to indicate that for part of the data set a dispersion correction was used. We do not believe the dispersion correction or program makes a meaningful di↵erence in this case, as shown by the low MAE demonstrated for a sample of PubChemQC species in Table S2, and therefore the PubChemQC and "wide range" results may be grouped together and treated as one large dataset. We note that a quadratic fit provides a better correlation description than a linear fit, where the former has a MAE of 2.47Å 3 compared to !B97X(D), while the latter has an MAE of 7.94Å 3 compared to !bB7X(D) ( Figure 5). As an aside, for comparison and the sake of completeness, we also tested the correlation between DFT and sTDA-xTB, a simplified time-dependent DFT procedure with a larger inherent basis set. 47 This yielded worse results, with a linear fit MAE of 15.85Å( Figure S3). We therefore conclude that although not as physically meaningful as an adjustment based on a related molecular property, we find a quadratic fit with zero intercept provides the best correction to GFN2/D4 for large polarizability calculations, and note the linear coe cient remains close to unity.  Comparing each smaller basis set to the largest set considered, aug-cc-pVTZ, we see similar results to the increasing pairwise comparison (Table 1). Again, simple linear regression analysis reveals slopes close to one and R 2 values of or nearly 1.00, even when comparing the largest basis set to a non-augmented double zeta basis set. The speed-ups are also worth noting, since even using a partially-augmented basis set (jun-cc-pVTZ) provides over a 2x speedup over the traditionally augmented set. Timings are shown in more detail in the box plot in Figure 7, where the range of calculation times for each basis set is shown to decrease substantially as the sets become smaller.  The large speed-ups provided by non or partially augmented basis sets, combined with lower risk of linear dependence issues, make them better, albeit less accurate, choices for polarizability calculations for large systems. The linearity between systematically larger ba-sis sets suggests that for species with large polarizabilities, the increase in accuracy of the magnitude of the polarizability is not substantial, and that a simple linear correlation coecient could be used to correct large polarizabilities found with smaller basis sets. Given the increased RMSE observed with cc-pVDZ, we suggest that for routine use on large molecules, non-augmented triple zeta basis sets, as used here, are an e cient balance of time and accuracy.

Conclusion
Based on our studies, GFN2/D4 appears best parameterized for species with polarizabilities less than 50Å 3 . In its current implementation GFN2/D4 does not compare favorably to DFT-computed polarizabilities in highly polarizable oligomers. We note that in addition to limiting GFN2's usefulness in its present form as a polarizability method for larger and/or conjugated systems, this raises important concerns about the D4 model's accuracy for similar systems, since the methods share the additive polarizability model.
The method's underestimation of polarizability values, which systematically grows as polarizability increases, indicates that the application of a quadratic scaling correction factor could provide a relatively simple solution to drastically improve the accuracy of large polarizability calculations. The presence of three significant outliers in the GFN2/D4 comparison data, all containing similar sulfur motifs, suggests the need to examine GFN2/D4 parameterization for such chemical structures. Beyond increasing GFN2/D4 usefulness for e cient polarizability calculations for large polarizability molecular screening applications, these improvements would notably also improve the accuracy of the D4 dispersion method for large ⇡-conjugated species. We also note that the semi-empirical additive model for GFN2/D4 remains more accurate than calculating polarizabilities using the related sTDA-xTB model.
With regard to the accuracy and e ciency of basis sets, our study suggests that using smaller, even non-augmented basis sets to save time and resources is appropriate for large po-larizability calculations. The substantial linear correlations seen between methods of varying size and levels of augmentation suggests that using a basis set such as cc-pVTZ with a linear scaling factor is appropriate for large polarizability molecules. By calculating polarizability in this manner for molecules with polarizabilities over 200Å 3 , we believe accuracies near the level of those achieved with an aug-cc-pVTZ basis set are attainable at nearly a six-fold speed-up and without the convergence issues we faced when attempting to use this basis set to calculate polarizabilities in this range. Using a basis set with a linear correction factor opens up the possibility of calculating highly accurate polarizabilities for increasingly large molecules using conventional DFT methods. Additional work will need to be done to test the limits of the polarizability magnitudes that can be accurately calculated in this manner.
We hope that the results of this study aid work where the calculation of large polarizability values is crucial. We note that for some applications, such as dielectric device development, the frequency dependence of polarizability should be considered, but is outside the scope of this work. While GFN2/D4 is not currently fit to provide accurate calculations for large polarizability values, we hope that after some minor corrections it will be a viable method for such applications and improve the accuracy of future dispersion correction methods. Considering the incredible e ciency, this would provide a valuable tool for future molecular screening studies of highly polarizable materials. Meanwhile, using lower-cost non-augmented basis sets with a correction factor vastly increases the number of potential molecular species for which highly accurate polarizabilities can be now obtained.