Fast and accurate quantum mechanical modeling of large molecular systems using small basis set Hartree-Fock methods corrected with atomcentered potentials
There has been significant interest in developing fast and accurate quantum mechanical methods for modeling large molecular systems. In this work, by utilizing a machine-learning regression technique, we have developed new low-cost quantum mechanical approaches to model large molecular systems. The developed approaches rely on using one-electron Gaussian-type functions called atom-centered potentials (ACPs) to correct for the basis set incompleteness and the lack of correlation effects in the underlying minimal or small basis set Hartree-Fock (HF) methods. In particular, ACPs are proposed for ten elements common in organic and bio-organic chemistry (H, B, C, N, O, F, Si, P, S, and Cl) and four different base methods: two minimal basis sets (MINIs and MINIX) plus a double-ζ basis set (6-31G*) in combination with dispersion-corrected HF (HF-D3/MINIs, HF-D3/MINIX, HF-D3/6-31G*), and the HF-3c method.The new ACPs are trained on a very large set (73832 data points) of non-covalent properties (interaction and conformational energies) and validated additionally on a set of 32048 data points. All reference data is of complete basis set coupled-cluster quality, mostly CCSD(T)/CBS. The proposed ACP-corrected methods are shown to give errors in the tenths of a kcal/mol range for non-covalent interaction energies and up to 2 kcal/mol for molecular conformational energies. More importantly, the average errors are similar in the training and validation sets, confirming the robustness and applicability of these methods outside the boundaries of the training set. In addition, the performance of the new ACP-corrected methods is similar to complete basis set DFT but at a cost that is orders of magnitude lower, and the proposed ACPs can be used in any computational chemistry program that supports effective-core potentials without modification. It is also shown that ACPs improve the description of covalent and non-covalent bond geometries of the underlying methods and that the improvement brought about by the application of the ACPs is directly related to the number of atoms to which they are applied, allowing the treatment of systems containing some atoms for which ACPs are not available. Overall, the ACP-corrected methods proposed in this work constitute an alternative accurate, economical, and reliable quantum mechanical approach to describe the geometries, interaction energies, and conformational energies of systems with hundreds to thousands of atoms.
Introduction<!>Theoretical background<!>ACP development process<!>Training and validation data sets<!>Technical details<!>Performance of ACPs for the training set<!>(i) Non-covalent interaction energies<!>(ii) Molecular conformational energies<!>(iii) Molecular deformation energies<!>Performance of ACPs for the validation set<!>Performance of ACPs for molecular geometries<!>Applications of ACPs developed for HF-3c<!>Summary and Outlook
<p>Quantum mechanical (QM) methods are an indispensable tool for understanding chemical phenomena. When combined with a nearly complete basis set, high-level wavefunction theory methods can predict various thermochemical and structural properties with an accuracy comparable to, or even better than, experiments. However, such approaches have limited applicability because their computational cost increases steeply with the size of the system. [1][2][3][4][5] This precludes high-level wavefunction methods from being applied to study chemical and biological processes involving large molecular systems, such as enzymatic catalysis, protein folding, supra-molecular host-guest complexation, and many others. [6][7][8][9][10][11][12][13][14] In the past few decades, a significant amount of effort has been devoted to developing efficient and accurate QM methodologies that can be applied to large molecular systems. [15][16][17][18][19][20][21][22][23][24][25][26][27][28] The application of QM modeling begins by selecting a set of approximations to solve the Schrödinger equation. One of the simplest QM approaches with a low computational expense is the Hartree-Fock (HF) method. However, HF has a major shortcoming in that, by definition, it does not calculate any correlation energy, which results in overly repulsive dispersion interactions, bond lengths that are too short, and the poor prediction of various other molecular properties. In addition, the cost and accuracy of any QM method is strongly dependent on the choice of basis set, the set of functions used to describe the system's molecular orbitals.</p><p>In HF, the computational cost scales roughly as the fourth power of the number of basis functions, with more sophisticated methods presenting an even steeper scaling. Calculations using either minimal or double-ζ basis sets are relatively inexpensive, but the use of these small basis sets introduces an additional error due to the insufficient number of basis functions. This basis set incompleteness error is severely detrimental to the method's accuracy. Therefore, even though minimal and double-ζ basis set HF offers a computationally inexpensive approach for modeling large molecular systems, a way needs to be devised to effectively mitigate the deleterious effect of missing electron correlation and basis set incompleteness error. 29,30 Many existing semi-empirical QM methods are based on approximations to minimal basis set HF. 31-integrals from the underlying minimal basis set HF approach while incorporating empirical parameters obtained by fitting to experimental or high-level theoretical reference data. These approximations substantially limit the accuracy of semi-empirical QM methods but in exchange reduce the computational cost below that of minimal basis set HF. Due to their reduced computational cost, semi-empirical QM methods have found extensive use in modeling large molecular systems. An example of a popular and more recent semi-empirical QM approach is the PM7 method of Stewart, which was also modified by Throssel and Frisch. 35,36 Another approach that is similar in spirit to conventional semi-empirical QM methods is the HF-3c 37 method proposed by Sure and Grimme. The HF-3c method uses three separate geometry-dependent formulas 38,39 to add energy corrections ("3c") for the various deficiencies of minimal basis set HF: one to account for some of the missing dispersion interactions, and two to mitigate the effects of basis set incompleteness errors. Several other techniques have been proposed in the literature [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54] , reflecting the interest in developing computationally inexpensive methods for large systems.</p><p>Finding a good compromise between cost and accuracy is critical when modeling large molecular systems. HF-3c is an example of a QM method that strikes a good balance between these two desirable characteristics. Even though the cost of HF-3c is higher than most semi-empirical QM methods, it is still orders of magnitude cheaper than nearly complete basis set wavefunction theory or density-functional theory (DFT) based methods. On the other hand, the accuracy of HF-3c in describing molecular structures and non-covalent interaction strengths is similar to large basis set DFT. 55 These features allow HF-3c to be applied for fast geometry optimizations, conformer exploration, and prediction of non-covalent interaction energies in fairly large systems, with sizes between many hundreds and a few thousand atoms. This allows the QM description of (small) biological systems (proteins, nucleic acids, carbohydrates, lipids) as well as supramolecular host-guest complexes. The downside of HF-3c is that it is unable to accurately describe thermochemical quantities such as bond breaking and formation energies. 56 Grimme and co-workers have applied the 3c correction to a few density functional approximations to address this problem. [57][58][59][60][61][62] Our previous works have shown that atom-centered potentials 63 (ACPs) offer a convenient means of improving the accuracy of HF and DFT based methods. [64][65][66][67][68][69][70][71][72][73][74][75][76] ACPs are one-electron potentials that share the same mathematical form as effective-core potentials 77,78 (ECPs) but do not replace any electrons. This allows ACPs to be used in most computational chemistry software packages without modifying the code.</p><p>Additionally, ACPs are an economical way of mitigating the errors in the underlying methodology, since using them incurs only a small additional cost. In a previous proof-of-concept work, we developed a single set of ACPs for the H, C, N, and O elements to mitigate the shortcomings of dispersion-corrected minimal basis set HF. 64 The parameters for the ACPs were obtained by fitting to a set of 9814 data points of noncovalent properties (interaction, conformational, and molecular deformation energies). In that work, we demonstrated the feasibility of the ACP correction approach by showing that ACPs developed for dispersion-corrected minimal basis set HF were able to accurately predict the mentioned non-covalent properties.</p><p>In this work, we build upon our previous study 64 and develop four ACP-corrected small basis set HF based methods. In all cases, the target applications are similar to those of HF-3c and our previous work, 64 namely structures and non-covalent interaction strengths. However, ACPs are developed for a larger set of atoms (H, B, C, N, O, F, Si, P, S, and Cl) than in our previous work, greatly increasing the applicability of the proposed methods. In addition, the use of the LASSO (Least Absolute Shrinkage and Selection Operator) regression [79][80][81] for fitting of ACP parameters greatly simplifies ACP development and allows using a training set about eight times as large and much more diverse than in earlier works, 64,[68][69][70]73 resulting in more robust and more widely applicable ACPs. Three of the four new ACP-corrected methods are based on HF with minimal (MINIs 82 and MINIX 37 ) or small double-ζ basis sets (6-31G* 83,84 ) and use Grimme's D3 38,85,86 correction to account for the missing dispersion in HF. In addition, we also present a set of ACPs designed to improve the performance of the HF-3c method. In addition, our intention with the HF-3c-ACP method is to overcome the limitation imposed by the fact that ACPs are available only for the ten elements mentioned above. Since HF-3c parameters are available for most elements in the periodic table, we expect HF-3c-ACP to reduce to HF-3c performance for the atoms for which ACPs have not been developed, which in general should be in the minority. The newly developed ACP-corrected methods are assessed using an extensive validation set, demonstrating their performance and robustness.</p><!><p>The procedure employed to develop the ACPs proposed in this work is similar to our earlier proofof-concept study 64 . The mathematical form of an ACP is:</p><p>where (2)</p><p>where the coefficients (𝑐) and exponents (𝜉) are adjustable parameters that are determined via a regularized least-squares fitting to reference data during ACP development (Section 2.2). The sum in Equations 2 and 3 runs over the total number (𝑁) of Gaussian-type functions defined for atom 𝛼 for the local and non-local potential terms. For ease of notation, we will represent the 𝑉 𝑙𝑜𝑐𝑎𝑙 α (𝑟) and 𝛿𝑉 𝑙 𝛼 (𝑟) together as:</p><p>If the ACP operator (Equation 1) with the functional form of Equation 4 is treated as a perturbative correction to any Hamiltonian then the first-order perturbation energy correction induced by the ACPs is:</p><p>where the sum in Equations 5 runs over the occupied molecular orbitals. Substituting the expressions from Equations 1 and 4 into Equation 5gives:</p><p>The</p><p>where 𝚫𝐄(𝝃) 𝑇 is the vector of ACP energy terms. It should be noted that Equation 7is only correct to first order in the ACP perturbation as the coefficients 𝑐 𝑙𝑛 𝛼 have an influence on the underlying wavefunction.</p><p>Equation 7 becomes exact only in the limit of 𝒄 → 𝟎, and it is approximately correct if the ACP coefficients are small in magnitude. The deviation between the 𝐸 𝐴𝐶𝑃 obtained using a self-consistent calculation with the corresponding ACP and the linear estimate in Equation 7, which assumes the coefficients have no influence on the underlying wavefunction, is called the non-linearity error.</p><!><p>ACPs have features that make them useful to develop energy corrections for a QM method (see Reference 69 for more details): (i) ACPs generate energy correction terms based on the molecular orbitals (Equation 5) and are wavefunction-dependent, which means they include information from the chemical environment and the electronic wavefunction, (ii) the angular projection operators in the potential (Equation 1) produces energy correction terms that are dependent upon the local anisotropic environment of a given atom, (iii) the exponential form of the ACPs (Equation 4) ensures that a given ACP produces a correction that decays exponentially with interatomic distances, (iv) ACPs can be used with any software that uses Gaussian-type basis sets and ECPs for geometry optimizations as well as calculation of energies and energy derived properties, and (v) the use of ACPs incurs only in a small computational cost (see Next, the exponents (𝜉 𝑙𝑛 𝛼 ), atoms (𝛼), and angular momentum channels (𝑙) are chosen, and the corresponding ACP energy terms (Δ𝐸 𝑙𝑛 α (𝜉 𝑙𝑛 α )) are calculated. The ACP energy term evaluation process is carried out by first obtaining SCF energy for every training set entry and then evaluating the ACP terms post-SCF. This approach speeds up the ACP energy term evaluation process which is important given the size of the training set and the number of ACP terms. Once the ACP energy terms have been computed for each target method/basis-set, exponents, angular momenta, and systems in the training set, the optimal ACP coefficients 𝑐 𝑙𝑛 𝛼 are determined using a regularized least-squares fit subject to a constraint on the sum of the absolute values of the coefficients. This constraint limits the magnitude of the ACP contributions and ensures that the correction arising from the ACP does not lead to significant non-linearity error, which would lead to disagreements between the predictions of our linear model (Equation 7) and the results obtained from the application of the ACP in an actual self-consistent calculation.</p><p>The training set is organized into subsets, corresponding to different molecular properties and data sources from the literature. Each subset of the training set is assigned a weight in the fitting procedure.</p><p>The weight of subset 𝑖 is calculated in the same way as in our previous work 64 :</p><p>where 𝑀 𝑖 is the average of the absolute value of the reference energies and 𝑁 𝑖 is the number of data points in subset 𝑖. These weights account for differences in reference data magnitude and number of points in the subsets. The error function minimized in the fit is the weighted root-mean-square-error (𝑤𝑅𝑀𝑆𝐸):</p><p>𝑁 𝑖 𝑗 𝑖 ∑ 𝑁 𝑖 𝑖 (9)   where 𝑗 are the data points in the 𝑖 th subset, 𝑦 𝑟𝑒𝑓,𝑗 𝑖 are the high-level reference energies for system 𝑗 in the 𝑖 th subset, and 𝑦 𝑚𝑒𝑡ℎ𝑜𝑑,𝑗 𝑖 are the energies of the underlying method for which the ACPs are being developed.</p><p>The LASSO (Least Absolute Shrinkage and Selection Operator) regression [79][80][81] , commonly employed in statistics and machine learning, is used to carry out the regularized least-squares fit. In LASSO, the 𝑤𝑅𝑀𝑆𝐸 in Equation 9 is minimized subject to the condition that 𝑙 1 -norm of the ACP coefficients does not exceed a certain bound chosen beforehand:</p><p>The LASSO method is used to limit the magnitude of the ACP coefficients. In addition, for a given constraint, LASSO automatically selects the best subset of ACP terms and discards the others, resulting in ACPs with fewer terms, which is beneficial because it curbs the computational cost of applying the ACPs.</p><!><p>The training set (Table 1) comprises non-covalent interaction energies, molecular conformational energies, and molecular deformation energies. This choice of training set properties is justified by the potential target applications of small basis set HF based methods, namely fast geometry optimizations and non-covalent interaction strengths in large systems as well as high-throughput 87 screening of conformers in combination with conformer search techniques [88][89][90] . These applications are useful, for instance, when performing exhaustive conformational searches of macrocyclic drugs [91][92][93][94][95][96] and other pharmaceutical candidates 97 , and studying biochemical processes like protein folding [98][99][100] and puckering of nucleotides 101,102 .</p><p>A successful method for non-covalent interactions must be able to accurately describe diverse noncovalent interaction motifs, which means that the training set must contain some of this diversity. For instance, the importance of π-π interactions is well-known in medicinal chemistry 103 , structural biology 104,105 , and organic electronics 106 . Such interactions also contribute to the stabilization of DNA 107,108 and proteins 109 , control the strength and specificity of drug-protein interactions 110 , and help in the rational design of supramolecules [111][112][113] . Other types of non-covalent interactions are also important in practice.</p><p>For example, Berka et al. reported that aliphatic-aliphatic (or hydrophobic) interactions between amino acid backbone chains are the most abundant in proteins, particularly in the hydrophobic active site. 114 Hydrophobic interactions also control the structure and properties of lipid bilayers 115,116 and selfassembled supramolecules 117 . On the other hand, hydrogen bonding is probably the most studied noncovalent interaction. 118,119 Several other stabilizing non-covalent interactions with potential chemical and biological applications have also been studied, including halogen bonding [216][217][218][219] , pnicogen bonding [220][221][222] , and anionic [120][121][122][123] interactions. Therefore, we designed our training set to contain representative candidates from the interaction types mentioned above. This also allows us to assess the strengths and weaknesses of the developed ACPs regarding each interaction type.</p><p>In most cases, the subsets of the training set, and the molecular geometries and reference data in them, were adopted from the literature. Occasionally, the reference energies were re-calculated at a higher-level to improve their quality. In each subset of the training set, data points involving molecules containing atoms other than H, B, C, N, O, F, Si, P, S, and Cl were excluded. A detailed list of all the subsets used for our ACP training set is given in Table S1 of the SI. The number of data points in each subset varies depending on the availability of benchmark data sets.</p><p>The training set used here is almost eight times larger than in our previous work 64 . In total, the training set comprises 73832 data points (167275 molecular geometries) calculated mainly with complete basis set wavefunction theory methods (Table S1 in SI). The total number of non-covalent interaction energies, molecular conformational energies, and molecular deformation energies are 19439, 44105, and 10288, respectively. The most abundant type of non-covalent interaction in the training set is hydrogen bonding. The mixed non-covalent interactions subset is second in abundance and features a mix of all common interactions found in large molecular systems. The large size of the training set ensures that no overfitting occurs when the least-squares fit is carried out.</p><p>In order to evaluate the performance and robustness of the new ACPs, we also assembled a validation set (Table 2), different from the training set, by compiling additional data sets from the literature. A detailed list of all the data sets included in the validation set is provided in Table S2 of the SI. In total, the validation set consists of 32047 high-level data points (92161 molecular geometries) calculated mainly with complete basis set wavefunction theory. The validation set contains 27811 non-covalent interaction energies and 4237 molecular deformation energies.</p><p>The structures and reference energies of all data points in the training and validation sets are given in the SI. In addition, the subsets that comprise the training and validation sets, grouped into categories to facilitate the analysis of the results, are listed in Tables 1 and 2. It should be noted that this subset categorization is in no way an exhaustive representation of the various types of systems in the training or validation set.  a) defined as the difference between the energy of the complex and the sum of the monomer energies. A negative interaction energy indicates the complex is more stable than the separated monomers. b) the reference data was recalculated in this work at the DLPNO-CCSD(T)/CBS level of theory (see SI for more details), at geometries reported in the literature. c) comprises non-covalently bound dimers where at least one of the monomers is negatively charged. d) the reference data was calculated in this work at CCSD(T)/CBS level using the same extrapolation method as in Reference 177.</p><p>e) defined as the difference between the energy of a particular conformer and a lower-energy conformer of the same molecule.</p><p>f) contains mostly conformational energies but also some molecular deformation energies. g) comprises negatively charged conformers.</p><p>h) defined as the difference between the energy of a molecule deformed along a particular normal mode and the energy of the same molecule at equilibrium. HW6Cl-anionic 215,216 , HW6Fanionic 215,216 , FmH2O10anionic 215,216 , SW49Bind345anionic 217 , SW49Bind6anionic 217 , Anionpi-anionic 218 , IL236-anionic 219 , DES15Kanionic 206 , NENCI-2021anionic 207 , CHAL336anionic 212 , XB45-anionic 211   +0.17 to +33.79 negatively charged molecules a) defined as the difference between the energy of the complex and the sum of energy of the monomers. A negative interaction energy indicates the complex is more stable than the separated monomers. b) comprises non-covalently bound complexes with at least one negatively charged monomer. c) defined as the difference between the energy of a particular conformer and a lower-energy conformer of the same molecule. d) only 30 systems used; we could not find the rest systems mentioned in Reference 228 in the supporting information e) only macrocyclic peptides considered.</p><p>f) only a subset from the PEPCONF 188 database for which reference data was recalculated at the DLPNO-CCSD(T)/CBS level of theory (see SI for more details). g) available reference data was calculated at LC-ωPBE-XDM/aug-cc-pVTZ level of theory. h) comprises negatively charged conformers.</p><!><p>Three sets of ACPs were developed for HF-D3 in combination with the minimal basis set MINIs, MINIX, and the double-ζ basis set 6-31G*. An additional set of ACPs was developed for HF-3c, which uses the MINIX basis set. The MINIX basis set was proposed at the same time as HF-3c 37  Twenty-nine ACP exponents (𝜉 𝑙𝑛 𝛼 ) were chosen, with values: 0.12 to 0.30 in 0.02 steps, 0.40 to 2.00 in 0.10 steps, and 2.50 to 3.00 in 0.50 steps. It should be noted that the choice of exponents is different than in our previous work 64 . A careful evaluation of the computational cost associated with the calculation of ACP energy term integrals (Equation 7) with low exponent functions (0.01 < 𝜉 𝑙𝑛 α < 0.11) suggested that ACPs containing exponents lower than 0.12 can lead to a significant increase in calculation time (especially for large molecules) compared to methods where ACPs are not applied. Therefore, choosing exponents higher than or equal to 0.12 ensures that the computational overhead is limited to a 10-30% increase relative to the uncorrected method.</p><p>The combination of ten atoms and twenty-nine exponents along with the various angular momentum channels results in 841 ACP terms for MINIs, 957 ACP terms for MINIX, and 1102 ACP terms for 6-31G*. For our training set of 73832 data points (167275 molecular geometries), ACP development required a total of 140678275 (HF-D3/MINIs), 160082175 (HF-D3/MINIX), and 184337050 (HF-D3/6-31G*) single-point energies. Combined with the self-consistent calculations used to evaluate the impact of non-linearity error (1338200) and the calculations on the validation set to evaluate the performance of the ACPs (737288), the total number of calculations for this project is 487172988. This complexity required the development of specialized software which we briefly describe next.</p><p>The parameters for the D3 dispersion correction used in this work correspond to those for the HF/aug-cc-pVTZ method with Becke-Johnson damping: 𝑠 6 = 1.0, 𝑠 8 = 0.9171, 𝑎 1 (𝐵𝐽) = 0.3385, and 𝑎 2 (𝐵𝐽) = 2.8830 Å. It should be noted that these D3 parameters are very close to those used in HF-3c 37 .</p><p>The ACP energy term evaluation and fitting processes were carried out using the dcp 230 and acpfit 231 packages available in our GitHub repository 232 . These programs automatize and collate all the data required for ACP development. For the LASSO regression, we used the local linearization plus active set method proposed by Osborne et al. 233 and implemented in octave/MATLAB by Schmidt 234,235  Once all the ACP energy terms for a particular target method were successfully computed, they were passed to the LASSO fit, resulting in an optimal set of ACPs for that method with minimum 𝑤𝑅𝑀𝑆𝐸 for a constraint of 25.0 au on the 𝑙 1 -norm of coefficients. The ACPs proposed in this work contain approximately 6-19 terms per atom, and they are designed to be paired with the specific method for which they were developed (i.e., HF-D3/MINIs, HF-D3/MINIX, HF-D3/6-31G*, or HF-3c), and are not transferable to other methods. The ACP coefficients and exponents for each method are provided in the SI. An example of the usage of ACPs in the Gaussian software is also given in the SI.</p><!><p>The optimal ACPs paired with their respective methods (HF-D3/MINIs, HF-D3/MINIX, HF-D3/6-31G*, and HF-3c) were applied self-consistently on the entire training set. The resulting non-covalent interaction energies, molecular conformational energies, and molecular deformation energies are compared to the corresponding reference data in Figures 1 and 2. The strip charts represent the error distribution as vertical lines for each method. The mean signed errors (MSEs) (open circle) and the standard deviations (SDs) of the errors (horizontal black lines) are also represented. The mean absolute errors (MAEs) and percentage change in the MAEs upon the application of ACPs (%∆MAE) for each method are listed on the right. A detailed breakdown of the errors for each method and subset can be found in Table S3 of the SI.</p><p>Table S3 of the SI also lists the deviation between the prediction of the ACP performance from our linear model (the LASSO fitting procedure) against the actual results from using the ACP in self-consistent calculations. This comparison, which measures the extent of non-linearity error, shows that the deviation between the MAEs predicted by the linear model and the self-consistent calculations is under 10% for most of the subsets of the training set, with only a few exceptions. This indicates that the 𝒍 𝟏 -norm constraint imposed in the LASSO fit was effective in preventing excessive non-linearity error and that the linear model used to develop ACPs is a faithful representation of their eventual performance as a correction method. In the following, the results obtained from the self-consistent application of ACPs are discussed for the different molecular properties in the training set.</p><!><p>The ACPs developed in this work have been trained on a wide range of non-covalent interaction types, including stacked and non-stacked π-π interactions, hydrophobic interactions, pnicogen bonding, halogen bonding, hydrogen bonding, and interactions of mixed and anionic nature. The proper description of each of these interactions is important for modeling large molecular systems, like proteins, where they operate co-operatively. 238    These two examples suggest that ACPs reduce the MAEs of the underlying methods when they are high and have a lesser impact on those subsets where the MAEs of the underlying method are already low.</p><p>Consequently, in a few rare instances the performance of an underlying method with low initial MAE can be negatively, but only slightly, impacted by the application of ACPs. This is the case for HF-3c-ACP applied to π-π interactions.</p><p>A particular limitation of all methods in this work is the performance for the anionic systems in the Anionic subset. This subset is challenging for basis sets like MINIs, MINIX, and 6-31G* because of the lack of diffuse functions required to properly describe negatively charged species. Figure 1 shows that the error spread and the SDs of the Anionic subset are larger than other interaction types. Because the Anionic subset was used in the training set, ACPs improve the performance of all four methods for anionic interaction energies, with MAE reductions of 19-49%. However, there is obvious room for improvement, and it is likely that it can only be achieved by the inclusion of diffuse basis functions, which would incur in an additional computational cost.</p><p>An interesting observation from Figure 1 (also  37,39 and SRB 37 corrections of HF-3c. Also, since the ACPs developed for HF-3c in most cases improve on HF-3c, ACPs provide additional error mitigation beyond that offered by gCP 37,39 and SRB 37 . On the other hand, ACPs developed for double-ζ basis set HF-D3 result in lower MAEs than those used in combination with minimal basis set HF-D3 for each interaction type, indicating that systematic improvement can be obtained by using ACPs with larger basis sets. However, going beyond double-ζ would lead to a significant increase in computational cost, and would result in methods with limited applicability for large molecular systems. 239 In this regard, a better alternative would be the development of ACPs for use with double-ζ DFT methods, an idea that is currently being explored in our group. 240 The performance of the proposed ACP based methods for the different interaction types and especially for the Mixed NCIs subset makes them promising for various applications. Keeping in mind our goal of designing low-cost approaches for modeling supramolecular and biological systems, we also assembled subsets and generated reference interaction energy data for prototypical non-covalently bound complexes relevant in biochemistry. Non-covalent interactions present in such systems are covered by the Biomolecule-Biomolecule subset. The Biomolecule-Biomolecule subset contains model systems representative of nucleotide-nucleotide interactions as well as protein fragments interacting with carbohydrates, nucleotides, drugs, water, and with other proteins. Such interactions are relevant in applications like protein folding 241,242 , protein structure refinement 243,244 , protein-ligand binding [245][246][247] , intercalation 248 , nucleobase stacking 249 , and protein hydration 250 , to name a few. Uncorrected minimal or double-ζ basis set HF-D3 in general overestimate the interaction energies in this subset, and application of the ACPs reduces this overestimation and decreases the error spread and SDs. Specifically, ACPs reduce the MAEs by about 74% (minimal basis set HF-D3) and 84% (double-ζ basis set HF-D3). HF-3c errors are centered around the zero-error average line with a relatively small spread, indicating that HF-3c is well suited for the complexes present in the Biomolecule-Biomolecule subset. Application of ACPs to HF-3c further reduces the MAE (by about 16%) except for a few systems: some nucleotide trimers and some nucleotide-amino acid complexes.</p><p>The Gas-Ligand subset comprises small molecules like CO2, CH4, and N2 interacting with benzene, coronene, polycyclic aromatic hydrocarbons, polyheterocyclic aromatic compounds, and other functionalized organic molecules. The complexes present in the Gas-Ligand subset are representative of potential applications in the areas of chemical sensing, gas storage, and gas separation. [251][252][253] By training our ACPs to this subset we expect to extend their applicability to the modeling of gas adsorption on various porous materials. 254 The application of ACPs decreases the MAEs of all considered methods for the Gas-Ligand by about 77% (minimal basis set HF-D3), 76% (double-ζ basis set HF-D3), and 18% (HF-3c). The bias of minimal or double-ζ basis set HF-D3 and HF-3c towards over-estimating the interaction energies are also reduced with the application of ACPs, resulting in lower error spread and SDs.</p><p>The Water-Water subset contains interaction energies of water dimers at various intermolecular separations as well as small water clusters (H2O)n with n=3-10. Potential target applications of ACPs trained against this subset are the modeling of aqueous environments, the study of surfaces of astrochemical interest [255][256][257][258][259][260] , as well as performing ab initio molecular dynamics simulations of water [261][262][263][264] . The ACPs improve the MAEs of minimal or double-ζ basis set HF-D3 and HF-3c methods for Water-Water by about 68% (minimal basis set HF-D3), 78% (double-ζ basis set HF-D3), and 19% (HF-3c).</p><p>Furthermore, Figure 1 shows that the large error spreads obtained with all underlying methods are significantly reduced by the ACPs, which bring the MSEs close to zero.</p><p>The last non-covalent interaction energy subset in the training set is BFSiPSCl, which contains complexes of monomers containing B, F, Si, P, S, and Cl. This subset extends the applicability of ACPs to systems like disulfide-linked proteins, covalent organic frameworks, functionalized silicon surfaces, and others. The application of ACPs results in a decrease in the MAEs of the BFSiPSCl subset by about 67%</p><p>(minimal basis set HF-D3), 70% (double-ζ basis set HF-D3), and 52% (HF-3c), with a decrease in the error spread and SDs in all cases. The drop in MAE observed for the HF-3c-ACP is more significant for this subset than all other interaction energy subsets, and the reduction is also close to that observed in the Pnicogen-bonding subset, indicating that perhaps the HF-3c parametrization is not as good for these systems as for the more "usual" non-covalent interactions in the previous sets.</p><p>Finally, we consider a few illustrative examples for which we compare the performance of our ACPcorrected methods with some commonly used DFT methods in combination with large basis sets. For this purpose, we use representative data sets from Mardirossian and Head-Gordon's benchmarking work 265 , for which nearly complete basis set DFT results have been reported in the literature. Specifically, we use the following sets: BzDC215 127 for π-π stacking interactions, HC12 131 for aliphatic-aliphatic interactions, S66x8 [144][145][146] and 3B-69-DIM 152 for interactions of mixed nature, SSI 154 and HSG 136,156 for biomoleculebiomolecule interactions, Water38 177 for water-water interactions, and Sulfurx8 181 for interactions involving S atoms. The reported MAEs (in kcal/mol) of the DFT methods with the very large def2-QZVPDD basis set as well as the MAEs of the ACP-corrected methods are shown in Table 3. The table</p><p>shows that the ACPs reduce the MAEs of minimal or double-ζ basis set HF-D3 and HF-3c methods in all the selected data sets and brings their MAE to a value close to or even lower than the almost complete basis set DFT methods. Therefore, Table 3 demonstrates that ACP-corrected methods have a performance similar to almost complete basis set DFT, but naturally at a cost that is reduced by orders of magnitude.</p><!><p>The purpose of the molecular conformational energy subsets of our training set is to inform the ACPs regarding how the potential energy surfaces of various molecules depend on the changes in rotatable bonds and torsional angles due to effects like π-conjugation, steric interactions, intramolecular hydrogenbonding, and electron repulsion. Our Small molecule conformational energy subset is a good representative of such interactions that can be used to assess the performance of ACPs for conformational energies. The Therefore, the ACPs developed for HF-3c offer a significant improvement beyond gCP 37,39 and SRB 37 for molecular conformation energies.</p><p>Finally, we compare the performance of our ACP-corrected methods relative to nearly complete basis set DFT results from the literature. For this, we consider a few representative data sets such as YMPJ 191 for amino acid conformers, SCONF 129,195 for carbohydrate-like conformers, BCONF 200 for butane-1,2-diol conformers, and PentCONF 201 for pentane conformers. The MAEs (in kcal/mol) of the DFT/def2-QZVPDD and the ACP-corrected methods for these data sets are shown in Table 3. Similar to non-covalent interaction energies, the application of ACPs reduces the MAEs of minimal or double-ζ basis set HF-D3 and HF-3c methods in all the selected data sets and are close to or even lower than the MAEs reported for the various functionals. Table 3 demonstrates that the proposed ACP-corrected methods are able to predict the conformational energies with an accuracy similar to large basis set DFT methods at a significantly lower computational cost.</p><!><p>The Deformation subset of the ACP training set contains energy differences between a molecule at its equilibrium geometry and the same molecule deformed along its various normal modes. Our intention with this subset is to improve the description of the molecular potential energy surfaces around the equilibrium geometries, and consequently improve the prediction of bond lengths and molecular geometries in general.</p><p>The fact that small basis set HF methods predict erroneous geometries is important for the study of large molecules like proteins, as discussed by Kulik et al. 269 and Schmitz et al. 270 Their findings suggest that small basis set HF methods without any correction give, in general, quite inaccurate protein structures. This is likely the reason why the HF-3c 37 method employs the semi-empirical SRB correction. In fact, the SRB correction itself was parametrized by fitting to the geometries of 107 small organic molecules computed at a higher level of theory.</p><p>The performance of our ACP-corrected methods for actual geometry optimizations is discussed in Section 3.4. The results the Deformation subset in Figure 2 already suggest that ACPs improve the prediction of molecular geometries substantially. On application of ACPs, the MAEs of all four methods for this subset are reduced by about 35-51% (minimal and double-ζ HF-D3) and 31% (HF-3c). Figure 2 shows that even though the decrease in the spread of errors using ACPs is modest, the under-estimation in the prediction of molecular deformation energies of the underlying methods is greatly corrected by the ACPs, and the MSEs as well as the SDs decrease. Molecular deformations that are farthest from equilibrium have relatively high reference energies and result in errors higher than ±5 kcal/mol.</p><p>Nevertheless, the application of ACPs predict individual errors that are lower than ±5 kcal/mol for 85% (or more) of the data points out of a total of 10,288.</p><!><p>The results regarding the application of ACPs to the systems in the validation set (Table 2) are presented in Figure 3. The figure includes the signed error distribution, MSEs, MAEs, and SDs of minimal or double-ζ basis set HF-D3 and HF-3c methods with and without ACPs, as well as the percentage change in MAEs upon application of ACPs (%∆MAE) for each method. A detailed breakdown of the errors by method and subset can be found in Table S4 of the SI. The results show that the when the ACPs are applied to the Mixed NCIs validation subset ("Mixed NCIs (val)"), the MAEs of all methods decrease, by 53% (minimal basis set HF-D3), 70% (double-ζ basis set HF-D3), and 37% (HF-3c). This reduction in MAEs is similar to what was observed for the mixed character interactions in the training set, confirming the robustness of the ACPs for non-covalent interactions when applied to systems outside the training set. One particular data set present in the Mixed NCIs (val) subset is BlindNCI 205 . Taylor et al. 205 reported an MAE of 0.34 kcal/mol with the M11/aug-cc-pVTZ method for the BlindNCI data set, which is almost equivalent to the MAE obtained with HF-D3/MINIs-ACP, HF-D3/MINIX-ACP, and HF-3c-ACP, and almost 28% higher than HF-D3/6-31G*-ACP (see Table 4). As in the training set, the performance of ACPs in the description of non-covalent interaction energies in the validation set is similar in quality to large basis set DFT methods. Two data sets (CE20 208,209 , WaterOrg  4), the MAEs of the CE20 data set with ACPs are close to or lower than most of the benchmarked DFT methods with a 6-311+G(3df,2p) basis set in the work of Chan et al. 208 Also, the B3LYP-D3(BJ)/6-31+G**-BSIP method that yields results that are close to B3LYP-D3(BJ)/aug-cc-pVQZ has an MAE of 0.44 kcal/mol for WaterOrg data set, which is close that predicted via HF-D3/6-31G*-ACP approach. 210 The assessment of the ACPs on the validation set also helps us understand what types of interactions are poorly represented in the training set. Based on the analysis performed with the validation set, these interactions are halogen bonding, chalcogen bonding, and close contact repulsions, as discussed below. It should be noted that an assessment of ACPs on interaction types such as π-π stacking, pnicogen bonding, and hydrophobic interactions, discussed earlier for the training set, was not possible in the validation stage because of the scarcity of high-level reference data in the literature.</p><p>For the halogen bonding interactions in the validation set ("Halogen-bonding (val)" subset), all methods in absence of ACPs show a large over-estimation of the interaction energies. The application of ACPs correct for this over-estimation and lead to a decrease in the MAEs by about 60-76%. Figure 3 shows that the MAEs of the minimal basis set HF-D3 and HF-3c methods without ACPs are almost three times higher than HF-D3/6-31G* (2.01 kcal/mol). Figure 3  bonding interactions are a difficult test not only for the methods considered in this work but also for many other electronic structure methods. For example, several dispersion-corrected DFT methods tested with ma-def2-QZVPP basis set have MAEs above 1 kcal/mol for the entire 212 data set. 212 In this context, Figure 3 suggests that the HF-3c method is the best suited among the minimal basis set HF methods for modeling chalcogen bonding interactions.</p><p>Steric repulsive interactions, even though found in some molecules that are forced to be in close contact due to the presence of other attractive interactions or external pressure, seldom occur naturally. 271 Repulsive interactions ("Repulsive contacts" subset) are captured well by the minimal or double-ζ basis set HF-D3 and HF-3c methods (MAEs of 0.59-0.89 kcal/mol) and only small improvements are seen with the application of ACPs. Specific subsets for repulsive interactions were missing from our training set.</p><p>Still, the slight reduction in the MAEs of minimal or double-ζ basis set HF-D3 and HF-3c methods observed with ACPs probably comes from some of the data sets in our training set that contain some data points with repulsive character (e.g. S22x5 136,142,143 , S66x8 [144][145][146] , S66a8 145 , A21x12 3,147,148 , and NBC10ext 128,136,[149][150][151] ). It should be noted that the MAEs of minimal or double-ζ basis set HF-D3 and HF-3c methods with and without ACPs are lower than the newly reparametrized PM6-D3H4R, DFTB3-D3H4R, PM6-D3H4X, and DFTB3-D3H4X methods (MAEs of 0.94-1.48 kcal/mol). These methods attempt to capture the repulsive interactions via the use of a repulsive energy correction term (parametrized against R160x6 213 and R739x5 214 data sets that constitute the validation Repulsive contacts subset) specifically designed for PM6 and DFTB3 methods with the D3H4 correction for dispersion and hydrogenbonding. 214,272 Next, we turn our attention to the (H2O)20 cluster (H2O20Bind10 216 data set), C60 dimer (C60dimer 221 data set), and Large molecule (L7 13,222,223 , S12L 9,11,223 , S30L 220 , Ni2021 224 data sets) subsets.</p><p>These subsets are a good test for ACPs as they contain non-covalently bound complexes that are relatively ACPs for HF-3c also reduce the MAEs by about 7-76%. The systems with large errors, higher than ±8 kcal/mol even after the use of ACPs, are, as expected, those that have very large reference energies. Table 4 shows a comparison of MAEs of various HF based approaches with large basis set DFT methods for the H2OBind10, C60dimer, L7, S12L, S30L, and Ni2021 data sets. It can be seen that, with the exception of Ni2021 data set, the MAEs in Table 4 for all other data sets shows that ACPs have a performance similar to large basis set DFT, making the approach particularly promising for modeling interaction energies in large molecular systems.</p><p>We now consider the results for the conformational energies evaluated with the subsets "Small molecule (val)" and "Proteinogenic". The Small molecule (val) subset contains conformational energies of various small organic and biaryl drug-like molecules. On the other hand, the Proteinogenic subset contains a collection of polypeptide conformers like tripeptides, peptides with disulfide linkages, macrocyclic peptides, and peptide sequences with associated bio-functionality. For the Proteinogenic subset, the improvement in MAEs seen on the application of ACPs is about 40%</p><p>(minimal basis set HF-D3), 49% (double-ζ basis set HF-D3), and 17% (HF-3c). Figure 3 shows that the spread of errors is more or less symmetric about the zero-error average line, except for some systems with errors higher than ±8 kcal/mol, corresponding to the disulfide linkages which were not present in the training set.</p><p>For non-covalent interaction energies of anionic interactions ("Anionic (val)" subset), although the</p><!><p>One attractive possible use of ACPs is for fast geometry optimizations. To gauge the performance of the various HF based methods for this task, we compared and analyzed the structures obtained after energy relaxation with those obtained using dispersion-corrected DFT methods with large basis sets. The 296 structures used for the test contain both non-covalently bound complexes and single molecules, ranging in size between 2 and 205 atoms. For single molecule structures, we used the equilibrium geometries of small organic molecules taken from our MOLdef data set, a variety of organic molecules taken from Reference 274, selected structures from LB12 57 and CLB18 275 data sets, and polypeptide structures from Reference 270. For non-covalently bound structures, we selected the equilibrium complex structures from the A21 148 , S66 144 , L7 222 , and S30L 220 data sets. Wherever high-level geometries were not available, we obtained reference geometries using dispersion-corrected DFT and a reasonably large basis set (CAM-B3LYP-D3(BJ)/6-311++G**) with the "tight" convergence criteria in Gaussian-16. All the optimized and reference geometries used for testing are provided in the SI. We used Kabsch's algorithm 276 to compare the optimized structures with the reference.</p><p>The results are summarized in To examine the source of the improvement in the molecular geometries upon application of ACPs, we calculated the average error in the intermolecular separation distances for the dimer complexes. We also compared the average error in a few selected bond lengths and angles for the single molecules. Table 5 shows that, on average, the intermolecular separation distances are improved, and the under-estimation in the separation distances yielded by the small basis set HF based methods is corrected by the ACPs, leading to better geometries for non-covalently bound complexes. Furthermore, the average deviations in bond lengths presented in Table 5 indicate that the average deviation of the small basis set HF based methods is between 0.002 Å to 0.092 Å for the selected bonds. It can be seen that the inclusion of ACPs with small basis set HF based methods leads to a better description of bond lengths as the polarity of a bond increases, leading to lower errors in bond lengths (except for HF-3c-ACP). A general improvement in the prediction of bond angles is also observed on application of ACPs. The combination of low average errors in bond lengths and angles leads to better overall geometries of single molecule structures. Despite overall good geometries for both single molecules and complexes, upon application of ACPs some tested geometries tend to have slightly higher RMSDs (greater than 0.7 Å) than the uncorrected methods due to slight deterioration in the bond lengths of C-H and C-C bonds. Such deviations relative to the reference geometry are visible in Table S5 of SI for some purely planar systems and peptides with highly flexible backbones. a) calculated as the distance between the centers of mass of each monomer.</p><p>b) excluding the geometries of the non-dimer complexes from the L7 data set for simplicity.</p><!><p>This section explores the use of HF-3c-ACP for modeling systems where most but not all atoms in the system have an associated ACP. For the atoms for which ACPs are not available, our intention is that HF-3c, which is the overall best of the underlying methods in this work will still give a reasonable description of the system. A particular example of an application where HF-3c-ACP could be used is in modeling metalloproteins where the atoms for which ACPs are unavailable are the metal ion(s) in the active site.</p><p>We calculated the interaction energies of two systems from Reference 224 using HF-3c and HF-3c-ACP to demonstrate the above idea. These two systems represent the adsorption of ethanol and benzene with different-sized cluster models of zeolite ZSM-5. We now explore the same idea by taking some subsets of the validation set and purposefully applying only part of the available ACPs so that not all atoms in the system have an associated correction. Table 6 presents a summary of the MAEs using various methods for two data sets from the validation set: the DES15K 206 set of non-covalent interaction energies (11474 data points) and the Torsion30 228 set of molecular conformational energies (2107 data points). The table shows that the MAEs of the HF/MINIX method are 4.83 and 0.92 kcal/mol for the DES15K and Torsion30 data sets, respectively. Using the HF-3c method, the MAE decreases for DES15K (2.14 kcal/mol) and increases to 1.18 kcal/mol for Torsion30.</p><p>Table 6 shows that using HF-3c-ACP but applying ACPs only to hydrogen and one of the non-hydrogen atoms indicates that ACPs improve the results progressively as the correction is applied to more atoms in the system. If only hydrogen and one of the non-hydrogen atoms are corrected, the performance of HF-3c-ACP is similar to HF-3c. If the correction is applied to hydrogen and two non-hydrogen atoms, most atoms are corrected and the MAEs decrease substantially and resemble HF-3c-ACP, in which all atoms receive an ACP. Therefore, we conclude that the application of ACPs is, in general, beneficial, and greater performance is obtained as more atoms receive an ACP, so the use of ACPs is recommended even in systems containing atoms for which ACPs are not available.</p><!><p>An important field of research in modern computational chemistry is the development of new quantum mechanical methods that are accurate and can be applied to model large molecular systems. Small basis set Hartree-Fock (HF) methods are orders of magnitude less expensive than more accurate nearly complete basis set wavefunction theory or DFT methods, but suffer from basis set incompleteness error and lack of electronic correlation. Provided these shortcomings can be addressed, such methods could be applied for modeling large molecular systems as well as for routine applications of fast geometry optimizations, conformational exploration, and prediction of non-covalent interaction strengths.</p><p>In this work, we show that HF with small and minimal basis sets can be effectively corrected by applying atom-centered potentials (ACPs, one-electron potentials similar to effective-core potentials) that are designed to correct for the inaccuracies in the underlying method. Four new sets of ACPs were developed for use with HF-D3 and small basis sets (MINIs, MINIX, 6-31G*) and HF-3c. The advantages of ACPs include that they can be used in most computational chemistry software packages without changes to the code and that they incur only a modest computational cost. The ACPs developed in this work apply to ten elements (H, B, C, N, O, F, Si, P, S, Cl), and our purpose is that the presented ACPs serve to address problems in organic chemistry and biochemistry. For the occasional system containing atoms for which no ACPs are available, we have shown that the improvement of the performance of the underlying method is progressive with the number of atoms where ACPs have been applied. Therefore, the use of ACPs is beneficial even if some atoms are not corrected, and in this case, we recommend the use of HF-3c-ACP.</p><p>We anticipate that the ACP based approaches developed in this work will allow efficient and accurate modeling of biomolecules such as proteins, nucleic acids, carbohydrates, lipids, and other molecules containing B, Si, and halogen atoms such as covalent organic frameworks, functionalized polyaromatic hydrocarbons, functionalized silicon surfaces, and more.</p><p>The ACPs were developed by using a large training set of 73832 data points calculated at a very high level of theory (CCSD(T)/CBS, in general). The training set contains a mixture of non-covalent interaction energies, molecular conformational energies, and molecular deformation energies. We expected that the size of the training set ensures the robustness and applicability of the ACPs. To test this, we validated the new ACPs on a validation set with 32047 data points. The assessment of minimal and double-ζ basis set HF-D3 and HF-3c methods, before and after the application of their corresponding ACPs showed that ACPs lower the MAEs of most subsets in the training set and that this good performance is carried over to the validation set with approximately the same performance in terms of average error. Relative to the uncorrected methods, ACP-corrected approaches improve the prediction of non-covalent interaction energies and molecular conformational energies. Furthermore, the addition of molecular deformation energies to the training set results in an improvement of the equilibrium molecular structures upon application of the ACPs. However, ACP-corrected methods showed relatively poor performance for some interaction types that were not part of the training set, such as chalcogen bonding or repulsive contacts, indicating that more diverse systems need to be included in the training set for greater robustness of the resulting methods.</p><p>Our analysis of representative data sets indicates that our ACP-corrected methods yield results similar to almost complete basis set DFT methods, naturally at a much lower computational cost.</p><p>Nonetheless, there remains a limitation regarding the description of negatively charged systems probably caused by the lack of diffuse functions in the basis sets employed. In spite of this, our ACPs offer a improvement even in this case. ACPs for small basis sets that include some diffuse functions are currently under development. We are also currently working on expanding the set of ACPs to DFT-D3 methods with small basis sets for prediction of accurate thermochemical properties along with non-covalent properties. Despite the limitation, we have shown that ACPs provide a way of developing methods that combine low-cost with robustness and wide applicability. We anticipate that ACPs will be useful to practitioners interested in modeling large systems or other time-intensive applications.</p>