Gaussian process regression adaptive density-guided approach: Toward calculations of potential energy surfaces for larger molecules.

We present a new program implementation of the Gaussian process regression adaptive density-guided approach [Schmitz et al., J. Chem. Phys. 153, 064105 (2020)] for automatic and cost-efficient potential energy surface construction in the MidasCpp program. A number of technical and methodological improvements made allowed us to extend this approach toward calculations of larger molecular systems than those previously accessible and maintain the very high accuracy of constructed potential energy surfaces. On the methodological side, improvements were made by using a Δ-learning approach, predicting the difference against a fully harmonic potential, and employing a computationally more efficient hyperparameter optimization procedure. We demonstrate the performance of this method on a test set of molecules of growing size and show that up to 80% of single point calculations could be avoided, introducing a root mean square deviation in fundamental excitations of about 3 cm-1. A much higher accuracy with errors below 1 cm-1 could be achieved with tighter convergence thresholds still reducing the number of single point computations by up to 68%. We further support our findings with a detailed analysis of wall times measured while employing different electronic structure methods. Our results demonstrate that GPR-ADGA is an effective tool, which could be applied for cost-efficient calculations of potential energy surfaces suitable for highly accurate vibrational spectra simulations.

Adaptive regularized Gaussian process regression for application in the context of hydrogen adsorption on graphene sheets.

We present a Gaussian process regression (GPR) scheme with an adaptive regularization scheme applied to the QM7 and QM9 test set, several protonated water clusters and specifically to the problem of atomic hydrogen adsorption on graphene sheets. For the last system our goal is to achieve good predictive accuracy with only a few training points. Therefore, we assess for these systems a self-correcting multilayer GPR model, in which the prediction is corrected by a chain of additional GPR models. In our adaptive regularization scheme, we impose no noise on the training data, but use an approach based on the data itself to account for its impurity. The strength of this strategy is that the data points are treated differently based on their importance and that the regularization can still be controlled by a single parameter. We assess how the accuracy of the prediction depends on this parameter. We can show that the new regularization scheme as well as the multilayer approach results in more robust predictors. Furthermore, we demonstrate that the predictor can be in good agreement with the density-functional theory results.

Massive Assessment of the Binding Energies of Atmospheric Molecular Clusters.

Quantum chemical studies of the formation and growth of atmospheric molecular clusters are important for understanding aerosol particle formation. However, the search for the lowest free-energy cluster configuration is extremely time consuming. This makes high-level benchmark data sets extremely valuable in the quest for the global minimum as it allows the identification of cost-efficient computational methodologies, as well as the development of high-level machine learning (ML) models. Herein, we present a highly versatile quantum chemical data set comprising a total of 11 749 (acid)1-2(base)1-2 cluster configurations, containing up to 44 atoms. Utilizing the LUMI supercomputer, we calculated highly accurate PNO-CCSD(F12*)(T)/cc-pVDZ-F12 binding energies of the full set of cluster configurations leading to an unprecedented data set both in regard to sheer size and with respect to the level of theory. We employ the constructed benchmark set to assess the performance of various semiempirical and density functional theory methods. In particular, we find that the r2-SCAN-3c method shows excellent performance across the data set related to both accuracy and CPU time, making it a promising method to employ during cluster configurational sampling. Furthermore, applying the data sets, we construct ML models based on Δ-learning and provide recommendations for future application of ML in cluster configurational sampling.

Thermicity of the Decomposition of Oxygen Functional Groups on Cellulose-Derived Chars.

The evolution of oxygen functional groups (OFGs) and the associated thermic effects upon heat treatment up to 800 °C were investigated experimentally as well as by theoretical calculations. A synthetic carbon with a carbonaceous structure close to that of natural chars, yet mineral-free, was derived from cellulose and oxidized by HNO3 vapor at different temperatures and for varied durations in order to generate char samples with different concentrations and distributions of OFGs. The functionalized samples were subjected to calorimetric temperature-programmed desorption measurements in correlation with an extensive effluent gas analysis, thereby focusing on the specific heat effects of individual OFG evolution. Interpretation of the experimental results was aided by density functional theory (DFT) calculations which allowed one to infer the thermal stability of different OFGs and the reaction energy associated with their evolution upon heating. Results showed that, with increasing temperature, H2O was released due to the loss of physisorbed water, the decomposition of clusters bound to carboxylic acids, and condensation reactions. The associated heat uptake amounted to about 100 kJ mol-1. Contrarily, the release of CO2, attributed to the decomposition of condensed acids, carboxylic acids, anhydrides, and lactones, resulted in a heat release of about 40 kJ mol-1. The most strongly pronounced thermic effects were detected for the release of CO, comprising highly exothermic effects due to the decomposition of condensed acids and carbonyls/quinones as well as endothermic effects attributed to anhydrides and phenols/ethers. Notably, anhydrides can be formed during the oxidative treatment as well as during heating by condensation of adjacent carboxylic acids. In the latter case, the two-step decomposition is overall highly exothermic, indicating the associated occurrence of pronounced carbon matrix rearrangements. DFT investigations suggest that these rearrangements not only affect the immediate OFG proximity but also involve several carbon sheets.

An automatized workflow from molecular dynamic simulation to quantum chemical methods to identify elementary reactions and compute reaction constants.

We present an automatized workflow which, starting from molecular dynamics simulations, identifies reaction events, filters them, and prepares them for accurate quantum chemical calculations using, for example, Density Functional Theory (DFT) or Coupled Cluster methods. The capabilities of the automatized workflow are demonstrated by the example of simulations for the combustion of some polycyclic aromatic hydrocarbons (PAHs). It is shown how key elementary reaction candidates are filtered out of a much larger set of redundant reactions and refined further. The molecular species in question are optimized using DFT and reaction energies, barrier heights, and reaction rates are calculated. The setup is general enough to include at this stage configurational sampling, which can be exploited in the future. Using the introduced machinery, we investigate how the observed reaction types depend on the gas atmosphere used in the molecular dynamics simulation. For the re-optimization on the DFT level, we show how the additional information needed to switch from reactive force-field to electronic structure calculations can be filled in and study how well ReaxFF and DFT agree with each other and shine light on the perspective of using more accurate semi-empirical methods in the MD simulation.

Can Small Polyaromatics Describe Their Larger Counterparts for Local Reactions? A Computational Study on the H-Abstraction Reaction by an H-Atom from Polyaromatics.

Hydrogen abstraction is one of the crucial initial key steps in the combustion of polycyclic aromatic hydrocarbons. For an accurate theoretical prediction of heterogeneous combustion processes, larger systems need to be treated as compared to pure gas phase reactions. We address here the question on how transferable activation and reaction energies computed for small molecular models are to larger polyaromatics. The approximate transferability of energy contributions is a key assumption for multiscale modeling approaches. To identify efficient levels of accuracy, we start with accurate coupled-cluster and density functional theory (DFT) calculations for different sizes of polyaromatics. More approximate methods as the reactive force-field ReaxFF and the extended semi-empirical tight binding (xTB) methods are then benchmarked against these data sets in terms of reaction energies and equilibrium geometries. Furthermore, we analyze the role of bond-breaking and relaxation energies, vibrational contributions, and post-Hartree-Fock correlation corrections on the reaction, and for the activation energies, we analyze the validity of the Bell-Evans-Polanyi and Hammond principles. First, we find good transferability for this process and that the predictivity of small models at high theoretical levels is way superior than any approximate method can deliver. Second, ReaxFF can serve as a qualitative exploration method, whereas GFN2-xTB in combination with GFN1-xTB appears as a favorable tool to bridge between DFT and ReaxFF so that we propose a multimethod scheme with employing ReaxFF, GFN1/GFN2-xTB, DFT, and coupled cluster to cope effectively with such a complex reactive system.

Implementation of the iterative triples model CC3 for excitation energies using pair natural orbitals and Laplace transformation techniques.

We present a pair natural orbital (PNO)-based implementation of CC3 excitation energies, which extends our previously published state-specific PNO ansatz for the solution of the excited state eigenvalue problem to methods including connected triple excitations. A thorough analysis of the equations for the excited state triples amplitudes is presented from which we derive a suitable state-specific triple natural orbital basis for the excited state triples amplitudes, which performs equally well for local and non-local excitations. The accuracy of the implementation is evaluated using a large and diverse test set. We find that for states with small contributions from double excitations, a T0 approximation to PNO-CC3 yields accurate results with a mean absolute error (MAE) for TPNO = 10-7 in the range of 0.02 eV. However, for states with larger double excitation contributions, the T0 approximation is found to yield significantly less accurate results, while the Laplace-transformed variant of PNO-CC3 shows a uniform accuracy for singly and doubly excited states (MAE and maximum error of 0.01 eV and 0.07 eV for TPNO = 10-7, respectively). Finally, we apply PNO-CC3 to the calculation of the first excited state of berenil at a S1 minimum geometry, which is shown to be close to a conical intersection. This calculation in the aug-cc-pVTZ basis set (more than 1300 basis functions) is the largest calculation ever performed with CC3 on excitation energies.

TURBOMOLE: Modular program suite for ab initio quantum-chemical and condensed-matter simulations.

TURBOMOLE is a collaborative, multi-national software development project aiming to provide highly efficient and stable computational tools for quantum chemical simulations of molecules, clusters, periodic systems, and solutions. The TURBOMOLE software suite is optimized for widely available, inexpensive, and resource-efficient hardware such as multi-core workstations and small computer clusters. TURBOMOLE specializes in electronic structure methods with outstanding accuracy-cost ratio, such as density functional theory including local hybrids and the random phase approximation (RPA), GW-Bethe-Salpeter methods, second-order Møller-Plesset theory, and explicitly correlated coupled-cluster methods. TURBOMOLE is based on Gaussian basis sets and has been pivotal for the development of many fast and low-scaling algorithms in the past three decades, such as integral-direct methods, fast multipole methods, the resolution-of-the-identity approximation, imaginary frequency integration, Laplace transform, and pair natural orbital methods. This review focuses on recent additions to TURBOMOLE's functionality, including excited-state methods, RPA and Green's function methods, relativistic approaches, high-order molecular properties, solvation effects, and periodic systems. A variety of illustrative applications along with accuracy and timing data are discussed. Moreover, available interfaces to users as well as other software are summarized. TURBOMOLE's current licensing, distribution, and support model are discussed, and an overview of TURBOMOLE's development workflow is provided. Challenges such as communication and outreach, software infrastructure, and funding are highlighted.

Assessment of the DLPNO Binding Energies of Strongly Noncovalent Bonded Atmospheric Molecular Clusters.

This work assesses the performance of DLPNO-CCSD(T0), DLPNO-MP2, and density functional theory methods in calculating the binding energies of a representative test set of 45 atmospheric acid-acid, acid-base, and acid-water dimer clusters. The performance of the approximate methods is compared to high level explicitly correlated CCSD(F12*)(T)/complete basis set (CBS) reference calculations. Out of the tested density functionals, ωB97X-D3(BJ) shows the best performance with a mean deviation of 0.09 kcal/mol and a maximum deviation of 0.83 kcal/mol. The RI-CC2/aug-cc-pV(T+d)Z level of theory severely overpredicts the cluster binding energies with a mean deviation of -1.31 kcal/mol and a maximum deviation up to -3.00 kcal/mol. Hence, RI-CC2/aug-cc-pV(T+d)Z should not be utilized for studying atmospheric molecular clusters. The DLPNO variants are tested both with and without the inclusion of explicit correlation (F12) in the wavefunction, with different pair natural orbital (PNO) settings (loosePNO, normalPNO, and tightPNO) and using both double and triple zeta basis sets. The performance of the DLPNO-MP2 methods is found to be independent of PNO settings and yield low mean deviations of -0.84 kcal/mol or below. However, DLPNO-MP2 requires explicitly correlated wavefunctions to yield maximum deviations below 1.40 kcal/mol. For obtaining high accuracy, with maximum deviation below ∼1.0 kcal/mol, either DLPNO-CCSD(T0)/aug-cc-pVTZ (normalPNO) calculations or DLPNO-CCSD(T0)-F12/cc-pVTZ-F12 (normalPNO) calculations are required. The most accurate level of theory is found to be DLPNO-CCSD(T0)-F12/cc-pVTZ-F12 using a tightPNO criterion which yields a mean deviation of 0.10 kcal/mol, with a maximum deviation of 0.20 kcal/mol, compared to the CCSD(F12*)(T)/CBS reference.

A Gaussian process regression adaptive density guided approach for potential energy surface construction.

We present a new iterative scheme for potential energy surface (PES) construction, which relies on both physical information and information obtained through statistical analysis. The adaptive density guided approach (ADGA) is combined with a machine learning technique, namely, the Gaussian process regression (GPR), in order to obtain the iterative GPR-ADGA for PES construction. The ADGA provides an average density of vibrational states as a physically motivated importance-weighting and an algorithm for choosing points for electronic structure computations employing this information. The GPR provides an approximation to the full PES given a set of data points, while the statistical variance associated with the GPR predictions is used to select the most important among the points suggested by the ADGA. The combination of these two methods, resulting in the GPR-ADGA, can thereby iteratively determine the PES. Our implementation, additionally, allows for incorporating derivative information in the GPR. The iterative process commences from an initial Hessian and does not require any presampling of configurations prior to the PES construction. We assess the performance on the basis of a test set of nine small molecules and fundamental frequencies computed at the full vibrational configuration interaction level. The GPR-ADGA, with appropriate settings, is shown to provide fundamental excitation frequencies of an root mean square deviation (RMSD) below 2 cm-1, when compared to those obtained based on a PES constructed with the standard ADGA. This can be achieved with substantial savings of 65%-90% in the number of single point calculations.

Machine learning for potential energy surfaces: An extensive database and assessment of methods.

On the basis of a new extensive database constructed for the purpose, we assess various Machine Learning (ML) algorithms to predict energies in the framework of potential energy surface (PES) construction and discuss black box character, robustness, and efficiency. The database for training ML algorithms in energy predictions based on the molecular structure contains SCF, RI-MP2, RI-MP2-F12, and CCSD(F12*)(T) data for around 10.5 × 106 configurations of 15 small molecules. The electronic energies as function of molecular structure are computed from both static and iteratively refined grids in the context of automized PES construction for anharmonic vibrational computations within the n-mode expansion. We explore the performance of a range of algorithms including Gaussian Process Regression (GPR), Kernel Ridge Regression, Support Vector Regression, and Neural Networks (NNs). We also explore methods related to GPR such as sparse Gaussian Process Regression, Gaussian process Markov Chains, and Sparse Gaussian Process Markov Chains. For NNs, we report some explorations of architecture, activation functions, and numerical settings. Different delta-learning strategies are considered, and the use of delta learning targeting CCSD(F12*)(T) predictions using, for example, RI-MP2 combined with machine learned CCSD(F12*)(T)-RI-MP2 differences is found to be an attractive option.

Approximate high mode coupling potentials using Gaussian process regression and adaptive density guided sampling.

We present a new efficient approach for potential energy surface construction. The algorithm employs the n-mode representation and combines an adaptive density guided approach with Gaussian process regression for constructing approximate higher-order mode potentials. In this scheme, the n-mode potential construction is conventionally done, whereas for higher orders the data collected in the preceding steps are used for training in Gaussian process regression to infer the energy for new single point computations and to construct the potential. We explore different delta-learning schemes which combine electronic structure methods on different levels of theory. Our benchmarks show that for approximate 2-mode potentials the errors can be adjusted to be in the order of 8 cm-1, while for approximate 3-mode and 4-mode potentials the errors fall below 1 cm-1. The observed errors are, therefore, smaller than contributions due to missing higher-order electron excitations or relativistic effects. Most importantly, the approximate potentials are always significantly better than those with neglected higher-order couplings.

Gaussian process regression to accelerate geometry optimizations relying on numerical differentiation.

We study how with means of Gaussian Process Regression (GPR) geometry optimizations, which rely on numerical gradients, can be accelerated. The GPR interpolates a local potential energy surface on which the structure is optimized. It is found to be efficient to combine results on a low computational level (HF or MP2) with the GPR-calculated gradient of the difference between the low level method and the target method, which is a variant of explicitly correlated Coupled Cluster Singles and Doubles with perturbative Triples correction CCSD(F12*)(T) in this study. Overall convergence is achieved if both the potential and the geometry are converged. Compared to numerical gradient-based algorithms, the number of required single point calculations is reduced. Although introducing an error due to the interpolation, the optimized structures are sufficiently close to the minimum of the target level of theory meaning that the reference and predicted minimum only vary energetically in the µEh regime.

Combining Accuracy and Efficiency: An Incremental Focal-Point Method Based on Pair Natural Orbitals.

In this work, we present a new pair natural orbitals (PNO)-based incremental scheme to calculate CCSD(T) and CCSD(T0) reaction, interaction, and binding energies. We perform an extensive analysis, which shows small incremental errors similar to previous non-PNO calculations. Furthermore, slight PNO errors are obtained by using TPNO = TTNO with appropriate values of 10-7 to 10-8 for reactions and 10-8 for interaction or binding energies. The combination with the efficient MP2 focal-point approach yields chemical accuracy relative to the complete basis-set (CBS) limit. In this method, small basis sets (cc-pVDZ, def2-TZVP) for the CCSD(T) part are sufficient in case of reactions or interactions, while some larger ones (e.g., (aug)-cc-pVTZ) are necessary for molecular clusters. For these larger basis sets, we show the very high efficiency of our scheme. We obtain not only tremendous decreases of the wall times (i.e., factors >102) due to the parallelization of the increment calculations as well as of the total times due to the application of PNOs (i.e., compared to the normal incremental scheme) but also smaller total times with respect to the standard PNO method. That way, our new method features a perfect applicability by combining an excellent accuracy with a very high efficiency as well as the accessibility to larger systems due to the separation of the full computation into several small increments.

Accuracy of Frequencies Obtained with the Aid of Explicitly Correlated Wave Function Based Methods.

We asses the basis set convergence of harmonic frequencies using different explicitly correlated wave function based methods. All commonly available CCSD(T) variants as well as MP2-F12 and MP4(F12*) are considered, and a hierarchy of the different approaches is established. As for reaction and atomization energies, CCSD(F12*)(T*) is a close approximation to CCSD(F12)(T*) and clearly superior to the other tested approximations. The used scaling for the triples correction enhances the accuracy relative to CCSD(F12*)(T) especially for small basis sets and is very attractive since no additional computational costs are added. However, this scaling slightly breaks size consistency, and therefore we additionally study the accuracy of CCSD(F12*)(T*) and CCSD(F12*)(T) in the context of calculating anharmonic frequencies to check if this causes problems in the generation of the potential energy surface (PES). We find a fast basis set convergence for harmonic and anharmonic frequencies. Already in the cc-pVDZ-F12 basis, the RMSD to the CBS limit is only around 4-5 cm-1.

Accuracy of Explicitly Correlated Local PNO-CCSD(T).

In recent years PNO-based local correlation methods have gained popularity since they allow Coupled Cluster (CC) calculations with reduced computational costs, yet only a few systematic studies concerning their accuracy are available, in particular for the explicitly correlated versions. In this work we take a deeper look at the explicitly correlated local PNO-CCSD(F12*)(T0) and PNO-CCSD(F12*)(T) methods. The first variant uses the so-called semicanonical triples correction (T0) which neglects off-diagonal elements in the occupied block of the Fock matrix. In PNO-CCSD(F12*)(T) this approximation is avoided by means of Laplace transformation techniques and convergence to the canonical results in the limit of no PNO truncation is restored. We assess the accuracy of both methods using well established benchmark sets for reaction energies and weak molecular interactions and take a look at a system with strong cooperative many-body effects. For reaction energies a close agreement with canonical methods is observed, and chemical accuracy can be reached. Also for weak intermolecular interactions the accuracy is easily controlled, and the methods even allow for improving existing benchmark data.

Atomic-batched tensor decomposed two-electron repulsion integrals.

We present a new integral format for 4-index electron repulsion integrals, in which several strategies like the Resolution-of-the-Identity (RI) approximation and other more general tensor-decomposition techniques are combined with an atomic batching scheme. The 3-index RI integral tensor is divided into sub-tensors defined by atom pairs on which we perform an accelerated decomposition to the canonical product (CP) format. In a first step, the RI integrals are decomposed to a high-rank CP-like format by repeated singular value decompositions followed by a rank reduction, which uses a Tucker decomposition as an intermediate step to lower the prefactor of the algorithm. After decomposing the RI sub-tensors (within the Coulomb metric), they can be reassembled to the full decomposed tensor (RC approach) or the atomic batched format can be maintained (ABC approach). In the first case, the integrals are very similar to the well-known tensor hypercontraction integral format, which gained some attraction in recent years since it allows for quartic scaling implementations of MP2 and some coupled cluster methods. On the MP2 level, the RC and ABC approaches are compared concerning efficiency and storage requirements. Furthermore, the overall accuracy of this approach is assessed. Initial test calculations show a good accuracy and that it is not limited to small systems.

Perturbative triples correction for local pair natural orbital based explicitly correlated CCSD(F12*) using Laplace transformation techniques.

We present an implementation of pair natural orbital coupled cluster singles and doubles with perturbative triples, PNO-CCSD(T), which avoids the quasi-canonical triples approximation (T0) where couplings due to off-diagonal Fock matrix elements are neglected. A numerical Laplace transformation of the canonical expression for the perturbative (T) triples correction is used to avoid an I/O and storage bottleneck for the triples amplitudes. Results for a test set of reaction energies show that only very few Laplace grid points are needed to obtain converged energy differences and that PNO-CCSD(T) is a more robust approximation than PNO-CCSD(T0) with a reduced mean absolute deviation from canonical CCSD(T) results. We combine the PNO-based (T) triples correction with the explicitly correlated PNO-CCSD(F12*) method and investigate the use of specialized F12-PNOs in the conventional triples correction. We find that no significant additional errors are introduced and that PNO-CCSD(F12*)(T) can be applied in a black box manner.

Explicitly correlated PNO-MP2 and PNO-CCSD and their application to the S66 set and large molecular systems.

We present our current progress on the combination of explicit electron correlation with the pair natural orbital (PNO) representation. In particular we show cubic scaling PNO-MP2-F12, and PNO-CCSD[F12] implementations. The PNOs are constructed using a hybrid scheme, where the PNOs are generated in a truncated doubles space, spanned by orbital specific virtuals obtained using an iterative eigenvector algorithm. We demonstrate the performance of our implementation through calculations on a series of glycine chains. The accuracy of the local approximations is assessed using the S66 benchmark set, and we report for the first time explicitly correlated CCSD results for the whole set and improved estimates for the CCSD/CBS limits. For several dimers the PNO-CCSD[F12] calculations are more accurate than the current reference values. Additionally, we present pilot applications of our PNO-CCSD[F12] code to host-guest interactions in a cluster model for zeolite H-ZSM-5 and in a calix[4]arene-water complex.

Auxiliary basis sets for density-fitted correlated wavefunction calculations: weighted core-valence and ECP basis sets for post-d elements.

We report optimised auxiliary basis sets for the resolution-of-the-identity (or density-fitting) approximation of two-electron integrals in second-order Møller-Plesset perturbation theory (MP2) and similar electronic structure calculations with correlation-consistent basis sets for the post-d elements Ga-Kr, In-Xe, and Tl-Rn. The auxiliary basis sets are optimised such that the density-fitting error is negligible compared to the one-electron basis set error. To check to which extent this criterion is fulfilled we estimated for a test set of 80 molecules the basis set limit of the correlation energy at the MP2 level and evaluated the remaining density-fitting and the one-electron basis set errors. The resulting auxiliary basis sets are only 2-6 times larger than the corresponding one-electron basis sets and lead in MP2 calculations to speed-ups of the integral evaluation by one to three orders of magnitude. The density-fitting errors in the correlation energy are at least hundred times smaller than the one-electron basis set error, i.e. in the order of only 1-100 µH per atom.