Combining quantum mechanics and machine-learning calculations for anharmonic corrections to vibrational frequencies
CCombining quantum mechanics and machine-learning calculations for anharmoniccorrections to vibrational frequencies
Julien Lam, Saleh Abdul-Al, and Abdul-Rahman Allouche Center for Nonlinear Phenomena and Complex Systems, Code Postal 231,Université Libre de Bruxelles, Boulevard du Triomphe, 1050 Brussels, Belgium ∗ Lebanese International University, Bekaa, Lebanon and International University of Beirut, Beirut, Lebanon Institut Lumière Matière, UMR5306 Université Lyon 1-CNRS,Université de Lyon, 69622 Villeurbanne Cedex, France † Several methods are available to compute the anharmonicity in semi-rigid molecules. However,such methods are not routinely employed yet because of their large computational cost, especiallyfor large molecules. The potential energy surface is required and generally approximated by a quar-tic force field potential based on ab initio calculation, thus limiting this approach to medium-sizedmolecules. We developed a new, fast and accurate hybrid Quantum Mechanic/Machine learning(QM//ML) approach to reduce the computational time for large systems. With this novel ap-proach, we evaluated anharmonic frequencies of 37 molecules thus covering a broad range of vibra-tional modes and chemical environments. The obtained fundamental frequencies reproduce resultsobtained using B2PLYP/def2tzvpp with a root-mean-square deviation (RMSD) of 21 cm − and ex-perimental results with a RMSD of 23 cm − . Along with this very good accuracy, the computationaltime with our hybrid QM//ML approach scales linearly with N while the traditional full ab initiomethod scales as N , where N is the number of atoms. I. INTRODUCTION
Recent advances in vibrational spectroscopy originatefrom progress at both experimental and theoretical lev-els. On the one hand, the development of ultrafastlasers has allowed for probing molecular dynamics on ex-tremely short time scales[1, 2]. On the other hand, com-putational spectroscopy has benefited from novel quan-tum mechanics calculations whose accuracy is compa-rable to experimental measurements[3–8]. In this re-search field, the combining role of experiments and sim-ulations is two-fold: Experiments help to assess the ac-curacy of the numerical calculations, which in turn, en-able for disentangling the complexity of the experimen-tal measurements[9–13]. To obtain additional progressin computational spectroscopy, the current challenge isto carry out calculations for large molecules while main-taining the accuracy obtained for the smaller ones[14–19].Such glass ceiling proves to be even more difficult to breakwhen it comes to computing the anharmonic correctionsto vibrational frequencies[14–19].Indeed, anharmonicity of the electronic potential isoften neglected and taken into account a posteriori us-ing corrective factors that are empirically determined tomatch experimental measurements. More rigorously, ex-plicit modeling of anharmonicity is also carried out usingseveral methods including VSCF[20], VSCF-PT2[21], cc-VSCF[22], VCI-P[23–25], VT2[26], GVPT2[27–31]. Be-cause of the required high-level of ab-initio calculations,the current calculations are computationally demand-ing for large molecules, in particular with double-hybrid ∗ [email protected] † [email protected] functionals which consist on mixing hybrid functionalswith an additional second-order perturbation theory cor-relation.Alternatively, instead of such a full ab-initio calcula-tion, several studies have been suggested. For solvatedmolecules, the computing time can be reduced using ahybrid QM//MM potential[18]. Other methods removemany derivatives terms, using only the 2-mode couplingquartic force constants[19]. But, the most popular ap-proach is a hybrid one where the geometric structureoptimization and harmonic calculations are computedat a high level of theory as the DFT theory and an-harmonic corrections are obtained with semi-empiricalmethods.[17]. More recently, we proposed to use molec-ular mechanic (MM) potential and in particular, theMMFF94[32] method which represents an even lowerlevel of theory[16]. The success of such a hybrid approachis two-fold. On the one hand, the accuracy is compara-ble to full ab initio calculations, and on the other hand,it allows for a considerable reduction of computationalcosts. However, MMFF94 method can only be employedfor organic molecules and an alternative approach mustbe derived.Machine-learning methods have been recently em-ployed to bridge the gap between molecular mechanic po-tentials and ab initio calculations by providing a frame-work to reach the accuracy of the latter while maintaininga low computational cost[33–35]. Four classes of meth-ods are usually considered: (1) Gaussian approxima-tion potentials[36, 37], (2) Kernel Ridge Regression[38–40], (3) Sparse Linear Regression[41–44] and (4) Artifi-cial Neural Network (ANN)[45–47]. Recent comparisonsbetween these different methods have also been carriedout[48, 49].In this context, we propose the use of machine-learningtechniques and in particular ANN potential. This way, a r X i v : . [ phy s i c s . c h e m - ph ] F e b our novel approach overcomes the limitations of the pre-vious one and is now versatile enough to compute anhar-monic corrections to vibrational frequencies in organic,as well as inorganic molecules. Such hybrid approachhas already been investigated in the literature and testedon few molecules but required the use of molecular dy-namics simulations[50]. In this paper, we wish to extendthe methodology by using only static calculations. Wewill first describe the methodology with particular detailson the choices that were made for the ANN parametersand the benchmarking molecules. Then, we will comparethe obtained results with experiments and also more ad-vanced ab initio calculations. II. METHODSA. The hybrid Quantum mechanics/Machinelearning approach
Anharmonic corrections to vibrational frequencies arecomputed within the explicit framework of the gen-eralized second-order vibrational perturbation theory(GVPT2)[28, 29]. We will only briefly outline key fea-tures of the approach before describing more thoroughlythe choices that were made for our own implementation.For a given molecule of N atoms, the potential energysurface is a function of the normal coordinates denoted Q . There are f normal coordinates with f = 3 N − fora non-linear molecule and f = 3 N − otherwise. WithinGVPT2, the potential energy surface (PES) is approxi-mated by a quartic force field using: V ( Q ) = V + V ( Q ) + V ( Q ) + V ( Q ) V ( Q ) = f (cid:88) i =1 h i Q i + 16 t iii Q i + 124 u iiii Q i V ( Q ) = f (cid:88) ij,i (cid:54) = j t ijj Q i Q j + 16 u ijjj Q i Q j + f (cid:88) ij,i Compute energy Using the calculated energies, compute cubic and quarticenergy derivatives Using these derivatives, compute anharmonic frequenciesby GVPT2 methodFIG. 1. Hybrid B2PLYP//ANN approach-pseudocode algorithmic representation of the entire hybrid method.Altogether in our approach, only N DFT calcula-tions of energies and forces per molecules are requiredto construct the quartic force field while for a full DFTapproach, it is necessary to perform (3 N − N − calculations of energies and forces, as implemented inthe Crystal program[60, 61] and N − Hessian calcula-tions as implemented in Gaussian[28]. However, the timeneeded to perform a Hessian calculation, even analyti-cally, is more expensive than a single-point energy+forcescalculation. If the Hessian matrix is computed numer-ically from forces, N − energies and forces cal-culations are needed for each Hessian and in overall N − N − energies and forces calculations arethen required to build the quartic force field.For simplicity, in the next sections, we refer to :• Full B2PLYP as a GVPT2 calculation where all thevalues are calculated using the B2PLYP/def2tzvppmethod via Gaussian software.• Hybrid B2PLYP/ANN as our hybrid method wherethe B2PLYP/def2tzvpp was used to compute theharmonic modes while the cubic and fourth onesare computed using the ANN potential.• Full ANN as calculation where all derivatives arecalculated using ANN potential via a custom-madeinterface to n2p2 library. The fundamental frequen-cies are calculated using GVPT2 method. B. Reference data The hybrid Quantum mechanics/Machine learning ap-proach is tested on 37 molecules: Water (H O), Kryp-ton difluoride (KrF ), Carbonyl selenide (OCSe), Sul-fur dioxide (SO ), selenium dioxide (SeO ), Carbondiselenide (CSe ), Titanium dioxide (TiO ), Ammo-nia (NH ), chlorine trifluoride (ClF ), Formaldehyde(H CO), Hydrogen peroxide (H O ), Gallium trichlo-ride (GaCl ), Aluminum chloride difluoride (AlF Cl),aluminum dichloride fluoride (AlFCl ), Zinc methy-lene (ZnCH ), Methane (CH ), Bromoform (CHBr ),silicon tetrabromide (SiBr ), Titanium tetrachloride(TiCl ), Carbon tetrabromide (CBr ), Formic acid(CH O ), Methanimine (CH NH), Methyl Alcohol(CH OH), Methane selenol (CH SeH), Vinyl bromide(C H Br), Chlorethene (CH CHCl), Selenium hexaflu-oride (SeF ), 1,2-Dibromoethane (C H Br ), methylgermane (GeH CH ), Ethyl bromide (C H Br), Cy-clopropane (C H ), Propylene oxide (C H O), pyruvicacid (C H O ), benzene (C H ), Thiophene (C H S),Dimethyl sulfate (C H O S) and Naphthalene (C H ).All of these molecules were specially chosen becausethey span over a large chemical space and becausetheir vibrational frequencies were experimentally mea-sured (taken from NIST Database)[62]. Altogether,our benchmarking is made on 407 experimental funda-mental frequencies available to the 37 molecules. Formolecules where it is possible to compute the fullB2PLYP/def2tzvpp anharmonic corrections to vibra-tional frequencies, we used them as a second way to assessthe accuracy of our approach. For this second compari-son, we worked with 371 frequencies available to the 34molecules. Three molecules could not be investigatedwith the full B2PLYP/def2tzvpp calculations because oflarge computational cost (C H ) or because of high sym-metry of the molecule whose anharmonic calculations cannot be performed with the employed version of Gaussian(SeF and C H ). III. BENCHMARKS For each molecule of our data set, after optimizationof the geometry using B2PLYP/def2tzvpp method, the N × geometries (where N is the number of atoms)are generated. Their B2PLYP/def2tzvpp energies andforces are calculated and used to build the correspondingANN potential. Fitting results are given in Table S3 ofsupporting information. The RMSD for forces used intraining varies from 5.0 to 23.0 meV/Å with an averagedvalue of 14.8 meV/Å. That of energies varies from 0.1to 0.6 meV with an averaged value of 0.21 meV. Thisexcellent result proves that our symmetry functions withonly 10 neurons in 2 hidden layers can reproduce, withexcellent accuracy the potential energy surface aroundthe equilibrium geometry.To validate the performance of our hybridB2PLYP/ANN GVPT2 approach, we carried outthe calculation of the fundamental frequencies on allmolecules of our data set. To compare the accuracyof our Hybrid method to that of the standard GVPT2approach (full B2PLYP), we calculated, when possiblethe fundamental frequencies of molecules using thismethod. Full B2PLYP is known as a very accuratemethod to study small- and medium-sized moleculesbut the computational cost rapidly becomes prohibitivefor large molecules[58]. All calculated and experimen-tal frequencies are given in Table S4 in supportinginformation. A. Assessment of the hybrid approach bycomparison with full B2PLYP one Frequencies obtained from our hybrid model are com-pared to results obtained from the full B2PLYP approachin Figure 2.a. The correlation coefficients R between ourhybrid ANN and full B2PLYP calculations are equal to0.9996 for fundamental frequencies and 0.934 for anhar-monic corrections. Averaging through our entire bench-mark set (Table I), we obtained a root mean squareddeviation (RMSD) of cm − for all the frequencies andit becomes cm − and cm − when respectively re-stricted to only the low and high frequencies. The un-signed absolute error (AUE) are about 13 cm − , 11 cm − and 19 cm − for all, low and high frequencies respec-tively. In addition, we also plotted in Figure 2.b thecorresponding deviation as a function of the full QM fre-quencies. Altogether, these results reflect the overall verygood quality of the agreement between the hybrid and thefull B2PLYP approaches.More precisely, Fig.3 shows the distribution of devia-tion between our hybrid calculations and full B2PLYPones. Most frequencies are obtained with a deviationsmaller than cm − . The distribution is peaked around0, where 59% of frequencies are predicted with a devia-tion between − to cm − , 77% between − and cm − and 81% between − and cm − .Following Fig.4, it appears that our results dependon the considered molecule. In particular, the largestmolecules including Cyclopropane(C H ) and MethylGerman GeH CH exhibit the largest error which is dueto the higher difficulty in obtaining a good descriptionof anharmonic effect on bend scissor HCH modes. Thiserror is not due to the number of neurons or to the GroupSymmetry Function. Indeed, the forces are reproducedwith a small RMS (11 meV/Å for C H and 14 meV/Åfor GeH CH ). Despite this specific case, the averagedRMSD on over all modes of all molecules remains verysmall.Altogether, our results are consistent with previousprecision obtained using MMFF94[16]. However, the ad-vantage of this novel approach is the large versatility asvirtually any molecules can be considered. TABLE I. Statistical errors in cm − using all frequencies, lowfrequencies ( < = − ) and High ones( > − ). Inthe column denoted "Reference", the line "B2PLYP" and"Experiment" correspond to cases where statistical errorsare computed using respectively B2PLYP and Experimentalfrequencies[62] as reference. All FrequenciesMethod Type Reference RMSD AUE ASE UMAXHybrid Fundamental B2PLYP 20.54 13.11 3.22 96.74Hybrid Fundamental Experiment 22.93 16.19 0.56 117.91B2PLYP Fundamental Experiment 19.80 13.07 -3.82 98.51ANN Fundamental B2PLYP 21.99 13.97 2.42 174.10ANN Fundamental Experiment 23.11 16.32 -1.00 86.93ANN Harmonic B2PLYP 18.15 11.88 -2.47 98.47Low FrequenciesMethod Type Reference RMSD AUE ASE UMAXHybrid Fundamental B2PLYP 19.60 11.62 5.08 98.51Hybrid Fundamental Experiment 20.80 14.26 3.69 117.91B2PLYP Fundamental Experiment 16.18 10.35 -1.15 85.86ANN Fundamental B2PLYP 23.11 14.03 3.15 174.10ANN Fundamental Experiment 21.79 15.16 1.27 85.55ANN Harmonic B2PLYP 19.10 12.21 -3.72 98.47High FrequenciesMethod Type Reference RMSD AUE ASE UMAXHybrid Fundamental B2PLYP 23.66 18.54 -3.53 62.88Hybrid Fundamental Experiment 29.40 23.23 -10.85 84.38B2PLYP Fundamental Experiment 28.60 21.95 -12.52 96.74ANN Fundamental B2PLYP 17.29 13.77 -0.27 47.62ANN Fundamental Experiment 27.40 20.55 -9.25 86.93ANN Harmonic B2PLYP 13.82 10.60 2.38 37.14 B. Assessment of the hybrid approach bycomparison with experiment Fig.5 shows the deviation between our hybrid calcu-lations and experimental measurements. Most frequen-cies are obtained with a deviation smaller than cm − with a distribution peaked around 0. Averaging throughour entire benchmark set, we obtained a root meansquared deviation (RMSD) of cm − for all the frequen-cies(Table I) and it becomes cm − and cm − whenrespectively restricted to only the low and high frequen-cies. The mean absolute errors are cm − , cm − and cm − for respectively all the frequencies, only the lowand only the high frequencies. This error can originatefrom two sources: (1) Dataset used to fit the ANN po-tential to the ab initio method (B2PLYP) and (2) Themethod (GVPT2) used to compute the fundamental fre-quencies. To distinguish between these two, the frequen-cies obtained with the full B2PLYP approach are com-pared with experimental ones which leads to a RMSDequal to 20, 16 and cm − for respectively all, low andhigh frequencies. The full ab initio calculations lead to anagreement with experiments that is similar to the hybridapproach. Altogether, this demonstrates that the ANNpotential does not really provide any additional error.Fig.6 shows that the errors depend on the consideredmolecules. The large error obtained for Cyclopropane is FIG. 2. (a) Scatter plots for the fundamentalfrequencies/cm − and the anharmonic corrections in in-set. Black points are for present data while the red linecorresponds to the ideal case. (b) Deviation between theHybrid QM/ML frequencies and the full QM frequencies.Blue dotted lines correspond to deviations of ± cm − . < - - t o - - t o - - t o - - t o - - t o 00 t o 1 01 0 t o 2 02 0 t o 3 03 0 t o 4 04 0 t o 5 0 > D i s t r i bu t i o n ( % ) Deviation (cm − )All Low High FIG. 3. Distribution of the vibrational frequency deviationbetween hybrid calculations and full ab initio calculations. due to the frequency mode discussed in the previous sec-tion. In addition, we observe a large error for ZnCH .This error is due to B2PLYP/GVPT2 approach not toANN potential. Indeed the RMSD obtained with fullB2PLYP/GVPT2 is about cm − . The CH2 waggingmode cannot be reproduced correctly with this approach.We note that there is a Fermi resonance between CH2wagging and CH2 scissor modes, treated, in GVPT2 ap-proach, variationally after removal of resonant term usingMartin et al. criteria[27]. A VCI calculation should be FIG. 4. Error in vibrational frequency obtained by hybridcalculations relative to full ab initio calculations for eachmolecule. more appropriate to study the anharmonic effect on thismolecule. An error in the experimental value cannot beexcluded either. Indeed, the . cm − of CH2 waggingmode in CH2 seems to be too small compared to the samemode in other molecules. < - - t o - - t o - - t o - - t o - - t o 00 t o 1 01 0 t o 2 02 0 t o 3 03 0 t o 4 04 0 t o 5 0 > D i s t r i bu t i o n ( % ) Deviation (cm − )All Low High FIG. 5. Distribution of the vibrational frequency deviationbetween hybrid calculations and experimental measurements. C. Assessment of full ANN approach We have seen above that hybrid B2PLYP/ANN ap-proach gave very good results as compared with the fullB2PLYP approach. Therefore, we further considered analternative approach, where equilibrium geometry, sec-ond, cubic and quartic derivatives are all calculated us-ing the PES obtained with the ANN potential. Whencomparing with the full B2PLYP approach, RMSD forfundamental frequencies is equal to , and cm − for respectively all, low and high frequencies, and theRMSD for harmonic frequencies is equal to , and cm − for respectively all, low and high frequencies (Ta-ble 1). This shows that our full ANN approach also leadsto very good results as compared with B2PLYP(See Fig-ures 7 and 8 ). The comparison between Figures 7 and 2shows very similar deviations in the vibrational frequencywith slightly better results for Hybrid B2PLYP/ANN ap-proach. Similarly, the comparison of our full ANN ap-proach with the experimental results shows very goodagreement. While the RMSD of full ANN approachis similar to that of our hybrid B2PLYP/ANN one, thisfull ANN approach leads, in general, to slightly less ac-curate results which arises from the errors already ob-tained in the ANN calculation of few very-low harmonicfrequencies. Such error could be coped with by increas-ing the number of data in the training set which natu-rally increases the computational cost. Altogether, thefull ANN approach remains a good alternative to hybridB2PLYP/ANN one, especially for very large moleculeswhere the harmonic B2PLYP calculation can take verylarge computational time. D. Timings Regarding the computational timing, we take C H S,relatively one of the large molecules in our benchmark set, FIG. 6. Error in vibrational frequency obtained by hybrid cal-culations and experimental measurements for each molecule. as an example. As illustrated in Figure 9, when using thefull B2PLYP treatment, it took 2240 hours in CPU time < - - t o - - t o - - t o - - t o - - t o 00 t o 1 01 0 t o 2 02 0 t o 3 03 0 t o 4 04 0 t o 5 0 > D i s t r i bu t i o n ( % ) Deviation (cm − )All Low High FIG. 7. Distribution of the vibrational frequency deviationbetween full ANN calculations and full ab initio calculations. and 72 hours in real-time on 40 cores. Calculations usingthe hybrid B2PLYP/ANN was considerably reduced as ittook only 400 hours in CPU time and 7 hours in real-time(using 100 cores for Data and only 8 cores for the trainingpart). The main advantage here is that for each molecule,calculations of the distorted structures were run in par-allel. It is evident that the computational time shouldgrow with the number of atoms. In the hybrid calcula-tion, the time grows as N × while it grows as ≈ (3 N ) for a full ab initio calculation.[16]. Hence, beside its verygood accuracy, the hybrid calculation approaches is tentimes faster than full ab initio methods. This significantreduction in the computational time enables us to carryout the anharmonic corrections to vibrational frequenciesof large molecules, which was not feasible before. E. Effect of data sets The main goal in this work was to construct an accu-rate PES through a fast method. This requires the useof data sets made of ab initio energies and forces. Wechoose a simple and general method to generate thesedata, making several small (important for modes withhigh frequencies) and large (important for modes withlow frequencies) displacement of each atom from the equi-librium position in the three Cartesian directions. With alarge number of displacements, the PES is well described.However, the computing time increases with the num-ber of displacements. It will therefore be necessary tofind a good compromise between precision and numberof geometries. So far, our results were obtained with N × displacements. To study the effect of the numberof data points, we calculated the fundamental frequenciesusing a ANN potential, fitted on N × displacements ( ± . , ± . , ± . Å), by removing the very large ones.First, for the neural network optimization, we obtaineda RMSD about 9.2 meV/Å for forces and 0.11 meV for FIG. 8. Error in vibrational frequency obtained by full ANNcalculations and full ab initio ones. energies, which is smaller than what was obtained witha larger data set (14.8 meV/Å for forces and 0.21 meVfor energies). This is simply because there are fewer datapoints to fit. Then, the RMSD between the full B2PLYPfundamental frequencies and the Hybrid ones is about24.5 cm − to be compared to 21.0 cm − observed using FIG. 9. CPU time/hours for full B2PLYP and hybridB2PLYP/ANN calculations. Time for each step of hybridB2PLYP/ANN ones. the N × geometries. The result with the smaller datasets is slightly less accurate than that of the larger onesbut the effect is not significantly important. By remov-ing geometries with very large displacements ( . Å ), theerror for modes with lowest frequenies grows, especiallyfor CH NH and C H O molecules, which explains theincrease in the value of RMSD. To explain this overall in-fluence of the dataset size, the correlation between the fit-ting error and the number of data points depends on thenumber of free parameters in the descriptor space. In thecase of the implemented ANN, the number of free param-eters is not much lower than the number of data points.As such, fitting with too few points leads to over fittingissues and we obtain better accuracy within the learn-ing set but worse outside of it. In addition, going from N × to N × data points was achieved by removingthe largest displacements which turns out to be crucial.Indeed, low-frequency modes require a large displacementand it has been mentioned in the literature that trunca-ture to the fourth order of the polynomial representationof the PES is not sufficient for a correct representationof the vibrational motion of theses modes[63]. We notethat by adding large displacements, the CPU comput-ing time is increased by a factor of 8/6, but it does notnecessarily increase the real computing time because thecalculations of energies and forces are independent andcan be done in parallel. IV. CONCLUSION To summarize, anharmonic contributions to vibra-tional frequencies are usually either computed withempirical correction factors or with a full quantummechanics calculations. We proposed an alternativeapproach which combines quantum chemistry calcula-tions with machine-learning potential. In particular,B2PLYP/def2tzvpp is employed to compute equilibriumgeometries and harmonic contributions while the an-harmonicity is obtained through the GVPT2 frameworkwhere derivatives of the PES are computed with a neu-ral network force-field potential. With this approach,we managed to reach a RMSD equal to 23 cm − whencompared to experimental results and to 21 cm − whencompared to full B2PLYP/def2tzvpp calculations. Yet,we note that in few cases the deviation can reach up to cm − . Further refinements regarding the size of thedatabase or the neural network methodology should helpfixing those outlying results. In terms of computationaltiming, we gained a factor of almost ten with the lattermethod which allows us to study large molecules. More-over, the use of a neural network force-field makes ourapproach transferable to molecules made of any types ofatoms. Furthermore, we demonstrate that the PES ob-tained with machine-learning is accurate enough to com-pute anharmonic frequencies. We implemented the ap-proach within the GVPT2 framework, but the obtainedPES could also be employed with other vibrational fre-quency methods such as VSCF, cc-VSCF, VCI-P. Fi-nally, beyond the calculation of vibrational frequencies,our work participates in the common effort to speed upab initio calculations and is therefore an additional evi-dence that machine-learning technique is a very effectivetool for material science. V. ACKNOWLEDGEMENT This work was only possible thanks to the generousgrants of computer time by P2CHPD (Université Lyon1) and the "Centre de calcul CC-IN2P3" at Villeurbanne,France. JL acknowledges financial support of the Fondsde la Recherche Scientifique - FNRS. VI. SUPPORTING INFORMATION AVAILABLE File available free of charge :• SuppInfo.pdf: List of all calculated frequencies,Group symmetry functions and the statestical pa-rameters of machine learning study. [1] L. Whaley-Mayda, S. B. Penwell, and A. Tokmakoff, J.Phys. Chem. Lett. , 1967 (2019). [2] M. Kondo, J. Nappa, K. L. Ronayne, A. L. Stelling, P. J. Tonge, and S. R. Meech, J. Phys. Chem. B , 20107(2006).[3] O. Christiansen, Phys. Chem. Chem. Phys. , 6672(2012).[4] S. M. Gruenbaum, C. J. Tainter, L. Shi, Y. Ni, and J. L.Skinner, J. Chem. Theory Comput. , 3109 (2013).[5] T. Joutsuka and A. Morita, J. Chem. Theory Comput. , 5026 (2016).[6] C. Bistafa, Y. Kitamura, M. T. C. Martins-Costa, M. Na-gaoka, and M. F. Ruiz-López, J. Chem. Theory Comput. , 4615 (2019).[7] B. Yang, S. Liu, and Z. Lin, Sci. Rep. , 1 (2017).[8] G. R. Medders and F. Paesani, J. Chem. Theory Comput. , 1145 (2015).[9] B. Schindler, L. Barnes, G. Renois, C. Gray, S. Cham-bert, S. Fort, S. Flitsch, C. Loison, A.-R. Allouche, andI. Compagnon, Nat. Commun. , 1 (2017).[10] B. Schindler, G. Laloy-Borgna, L. Barnes, A.-R. Al-louche, E. Bouju, V. Dugas, C. Demesmay, and I. Com-pagnon, Anal. Chem. , 11741 (2018).[11] J. Gaynor, A. Petrone, X. Li, and M. Khalil, J. Phys.Chem. Lett. , 6289 (2018).[12] J. Kozubal, T. R. Heck, and R. B. Metz, J. Phys. Chem.A , 23 (2019).[13] S. Brünken, F. Lipparini, A. Stoffels, P. Jusko,B. Redlich, J. Gauss, and S. Schlemmer, J. Phys. Chem.A , 8053 (2019).[14] B. Hirshberg, L. Sagiv, and R. B. Gerber, J. Chem.Theory Comput. , 982 (2017).[15] I. Ivani, V. Baumruk, and P. Bouř, J. Chem. TheoryComput. , 2095 (2010).[16] L. Barnes, B. Schindler, I. Compagnon, and A.-R. Al-louche, Journal of Molecular Modeling , 285 (2016).[17] T. K. Roy, R. Sharma, and R. B. Gerber, Phys. Chem.Chem. Phys. , 1607 (2016).[18] K. Yagi, K. Yamada, C. Kobayashi, and Y. Sugita, J.Chem. Theory Comput. , 1924 (2019).[19] M. Keçeli and S. Hirata, The Journal of Chemical Physics , 134108 (2011).[20] J. M. Bowman, J. Chem. Phys. , 608 (1978).[21] K. Yagi, T. Taketsugu, K. Hirao, and M. S. Gordon, J.Chem. Phys. , 1005 (2000).[22] J. O. Jung and R. B. Gerber, J. Chem. Phys. , 10332(1996).[23] P. Carbonnière, A. Dargelos, and C. Pouchan, Theoret-ical Chemistry Accounts , 543 (2010).[24] A. Erba, J. Maul, M. Ferrabone, P. Carbonnière,M. Rérat, and R. Dovesi, J. Chem. Theory Comput. , 3755 (2019).[25] A. Erba, J. Maul, M. Ferrabone, R. Dovesi, M. Rérat,and P. Carbonnière, J. Chem. Theory Comput. , 3766(2019).[26] H. H. Nielsen, Rev. Mod. Phys. , 90 (1951).[27] J. M. L. Martin, T. J. Lee, P. R. Taylor, and J. François,J. Chem. Phys. , 2589 (1995).[28] V. Barone, J. Chem. Phys. , 014108 (2005).[29] V. Barone, M. Biczysko, and J. Bloino, Phys. Chem.Chem. Phys. , 1759 (2014).[30] Y. Cornaton, M. Ringholm, O. Louant, and K. Ruud,Phys. Chem. Chem. Phys. , 4201 (2016).[31] A. P. Charmet and Y. Cornaton, Journal of MolecularStructure , 455 (2018).[32] T. A. Halgren, J. Comput. Chem. , 490 (1996).[33] J. Behler, J. Chem. Phys. , 170901 (2016). [34] J. Schmidt, M. R. G. Marques, S. Botti, and M. A. L.Marques, npj Comput. Mater. , 1 (2019).[35] J. Hill, G. Mulholland, K. Persson, R. Seshadri,C. Wolverton, and B. Meredig, MRS Bull. , 399(2016).[36] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi,Phys. Rev. Lett. , 136403 (2010).[37] A. P. Bartók and G. Csányi, Int. J. Quantum Chem. ,1051 (2015).[38] D. Hu, Y. Xie, X. Li, L. Li, and Z. Lan, J. Phys. Chem.Lett. , 2725 (2018).[39] K. Vu, J. C. Snyder, L. Li, M. Rupp, B. F. Chen, T. Khe-lif, K.-R. Müller, and K. Burke, Int. J. Quantum Chem. , 1115 (2015).[40] M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. vonLilienfeld, Phys. Rev. Lett. , 058301 (2012).[41] A. Seko, A. Takahashi, and I. Tanaka, Phys. Rev. B ,054113 (2015).[42] A. Seko, A. Takahashi, and I. Tanaka, Phys. Rev. B ,024101 (2014).[43] A. Takahashi, A. Seko, and I. Tanaka, J. Chem. Phys. , 234106 (2018).[44] A. Takahashi, A. Seko, and I. Tanaka, Phys. Rev. Ma-terials , 063801 (2017).[45] J. Behler and M. Parrinello, Phys. Rev. Lett. , 146401(2007).[46] J. Behler, Phys. Chem. Chem. Phys , 17930 (2011).[47] N. Artrith and A. Urban, Comput. Mater. Sci , 135(2016).[48] W. Li and Y. Ando, Phys. Chem. Chem. Phys , 30006(2018).[49] Y. Zuo, C. Chen, X. Li, Z. Deng, Y. Chen, J. Behler,G. Csányi, A. V. Shapeev, A. P. Thompson, M. A.Wood, and S. P. Ong, “A performance and cost as-sessment of machine learning interatomic potentials,”(2019), arXiv:1906.08888.[50] M. Gastegger, J. Behler, and P. Marquetand, Chem. Sci. , 6924 (2017).[51] K. Yagi, K. Hirao, T. Taketsugu, M. W. Schmidt, andM. S. Gordon, J. Chem. Phys. , 1383 (2004).[52] J. Behler, J. Chem. Phys. , 074106 (2011).[53] A. Singraber, T. Morawietz, J. Behler, and C. Dellago,J. Chem. Theory Comput. , 3075 (2019).[54] A. Singraber, J. Behler, and C. Dellago, J. Chem. TheoryComput. , 1827 (2019).[55] S. Grimme, J. Chem. Phys. , 034108 (2006).[56] F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. , 3297 (2005).[57] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria,M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone,B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Cari-cato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino,G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toy-ota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima,Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Mont-gomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark,J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov,R. Kobayashi, J. Normand, K. Raghavachari, A. Ren-dell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi,N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B.Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts,R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi,C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma,V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannen- berg, S. Dapprich, A. D. Daniels, d. Farkas, J. B. Fores-man, J. V. Ortiz, J. Cioslowski, and D. J. Fox, “Gaus-sian 09 Revision D.01,” (2009), gaussian Inc. WallingfordCT 2009.[58] M. Biczysko, P. Panek, G. Scalmani, J. Bloino, andV. Barone, J. Chem. Theory Comput. , 2115 (2010).[59] A.-R. Allouche, “igvpt2 is a program for computing an-harmonic corrections to vibration frequencies, based onforce field expansion of the potential energy surface innormal mode coordinates.” (2016).[60] A. Erba, J. Maul, M. Ferrabone, P. Carbonnière,M. Rérat, and R. Dovesi, J. Chem. Theory Comput. , 3755 (2019).[61] A. Erba, J. Maul, M. Ferrabone, R. Dovesi, M. Rérat,and P. Carbonnière, J. Chem. Theory Comput. , 3766(2019).[62] M. Jacox, NIST Chemistry WebBook, NIST StandardReference Database Number 69, National Institute ofStandards and Technology , edited by P. J. Linstromand W. G. Mallard (National Institute of Standardsand Technology(retrieved August 23, 2019), Gaithers-burg MD, 20899, 2005).[63] N. Inostroza, X. Huang, and T. J. Lee, TheJournal of Chemical Physics135