A Comparison of Medium-Sized Basis Sets for the Prediction of Geometries, Vibrational Frequencies, Infrared Intensities and Raman Activities for Water

Optimized geometries, vibrational frequencies, as well as infrared intensities and Raman activities were calculated for water (H2O) utilizing popular quantum mechanical approaches. Here, density functional theory (DFT) calculations were performed using the B3LYP (Becke, three-parameter, Lee-Yang-Parr) functional, as well as ab initio calculations using second-order Møller-Plesset (MP2) perturbation theory and coupled-cluster with single, double and perturbative triple excitations [CCSD(T)] levels of theory were used. We assess and benchmark the performance of 69 different atomic orbital basis sets including various popular families of medium-sized basis sets typically of two to four zeta quality and differing levels of augmentation by polar and diffuse functions. The basis sets range from the commonly adopted Pople-style (6-31G & 6-311G), Dunning’s correlation consistent (cc-pV(n+d)Z & aug-cc-pV(n+d)Z, as well as Truhlar’s calendar variations, Jensen’s polarization consistent (pc-n & aug-pc-n), Ahlrichs (def2-…), Sapporo’s and Karlsruches as well as atomic natural orbitals (ANOs) such as NASA Ames (ANOn), Neese-style, and Roos-style. We also compare several basis sets specifically designed to calculate vibrational and electronic properties, including the Sadlej-pVTZ (and LPol-X families), as well as SNS families of Barone. The results are compared to experimental values where available, or calculations performed with 5 or 6 zeta-level (e.g., cc-pV6Z). The performance of each family of basis sets is discussed in terms of their accuracy (and pitfalls), as well as computational resource scaling and efficiency. The Def2 basis family performs very well overall, yielding more accurate results with lower runtimes than traditional basis sets. ‘May’ basis sets also provide accurate predictions of vibrational frequencies at significantly lower costs. Raman activities can be accurately calculated using MP2 under harmonic approximation with several ‘spectroscopic’ families performing well.


Introduction
Predictions of molecular properties, such as geometries, vibrational frequencies, infrared intensities and Raman activities have become powerful tools for chemists and physicists performing spectroscopic analysis of substances. Particularly involving the characterization  [1,2]. Commonly employed quantum mechanical approaches to calculation these properties include density functional theory (DFT) and ab initio methods, of which extensive studies have been done investigating accuracy of their predictions [3][4][5][6][7][8][9][10][11]. However, only a few studies exist systematically comparing the performance of a limited number of differing atomic orbital basis sets to accurately calculate these properties [12]. Moreover, careful selection from the plethora of basis set choices now easily accessible (for example through the EMSL database, https://bse.pnl.gov/bse/portal) can lead to optimization of the approach implemented based upon computational resources and required accuracy levels [12][13][14]. Since computational resources are often constrained, Poplestyle basis sets are often employed since they are often not demanding, and widely available in commercial packages, however, they are often out-performed by other basis sets at similar computational cost. On the other hand, correlation consistent basis sets are widely available and -although computationally heavy -they do systematically converge on accurate results when either 5/6 zeta basis sets are used, or complete basis set extrapolation schemes are incorporated (where each zeta, ζ, typically approximately doubles the number of atomic orbitals). Since computational costs for these calculations scale with N 4 to N 7 , the runtimes for even mediumsized molecules under 5/6 zeta become unfeasible. However, these large calculations don't necessarily produce the most accurate predictions. The results from several basis set families can be extrapolated to complete basis set (CBS) limit via a series of calculations performed at lower zeta level [15][16][17]. This study compares the results of several medium-sized basis sets across three levels of theory, in order to identify which basis sets can be reliably utilized to predict molecular properties (such as geometries, frequencies, IR intensities and Raman activities) maximizing the precision/accuracy for a given set of computational resources. Here, we only focus on preliminary results for water, but the findings are in line with those observed for other molecules we have investigated. Both DFT and ab initio calculations were perfomed, using Becke's three parameter exchange funtional [18] along with Lee-Yang-Parr's correlation functional [19] (B3LYP), secondorder Møller-Plesset (MP2) perturbation theory [20] and coupled-cluster with single, double and perturbative triple excitations [CCSD(T)] [21,22]. For each level of theory 69 different atomic orbital basis sets were used varying from Pople-style (6-31G & 6-311G) [23], Dunning's correlation consistent (cc-pV(n+d)Z & aug-cc-pV(n+d)Z) [24], as well as Truhlar's calendar variations [25][26][27][28][29][30], Jensen's polarization consistent (pc-n & aug-pc-n) [31][32][33], Alrich's (def2-...) [34], Sapporo's and Karlsruhe's, as well as atomic natural orbitals (ANOs) such as NASA Ames (ANOn), Neese-style [35], and Roos-style. We also compare several basis sets specifically designed for calculating vibrational and electronic properties, including the Sadlej-pVTZ (and LPol-X families) [36,37], as well as SNS & NO7 families of Barone [38,39] and a new generation of Thakkar DZ basis sets, NLO-1 [40]. Calculated values are compared to both experimental values [41][42][43] where available, as well as various convergences found either running large zeta calculations or through extrapolation. Though many medium-sized basis sets may converge on the solutions provided by extrapolation schemes and high-zeta results, this does not guarantee, and in fact rarely results in the prediction of molecular properties being accurately calculated due to inherent flaws in the level of theory taken. Therefore, minor improvements at the cost of significantly larger runtimes are rarely worthwhile. It is of note that this study does not account or attempt to correct for errors arising from inadequacies in the levels of theory described here, although it is possible to do so (e.g., by correcting for the truncation of the theoretical approaches, core-correlation effects, or relativistic effects which need to be taken into account to reproduce frequencies to within 1 cm −1 ). However, due to fortuitous cancellation of errors, the CBS limits, or cc-pV6Z results are often within 5-10 cm −1 of the corrected values [44]. Here, we report on which basis sets can be reliably utilized to predict molecular properties (such as geometries, frequencies, IR intensities and Raman activities) maximizing the precision/accuracy IOP Conf. Series: Journal of Physics: Conf. Series 1290 (2019) 012013 IOP Publishing doi:10.1088/1742-6596/1290/1/012013 3 for a given set of computational resources. The need for high zeta basis sets, or extrapolated results can often be avoided, if the basis sets are more carefully considered for the property being evaluated, and more accurate results reliably obtained at a lower computational cost. We briefly address the benefits of using anharmonic corrections to determine vibrational frequencies versus the use of scaling factors for a few select cases.

Methods
Calculations were done using the General Atomic and Molecular Electronic Structure System, GAMESS(US) [45][46][47], using the 14 FEB 2018 (R1) version of Gamess(US). All calculations were run using spherical harmonics; Pople-style basis sets were typically ran using cartesian functionals unless they incorporated f-type polarization functions, where spherical harmonics were used. The segmented versions of polarization consistent basis sets were used, it should also be noted that within the Def2 family QZVPD is equivalent to QVZPPD and that QZVP is equivalent to QVZPP. Tight convergence criteria were used throughout (e.g., integral cutoff 10 −12 , primitive cutoff 10 −25 , 10 −7 SCF convergence, gradient convergence 10 −6 for geometries). Hessian and Raman calculations were run with semi-numerical (B3LYP and MP2) or fullynumerical gradients [CCSD(T)] under the double harmonic approximation with projected frequencies and two displacements in each cartesian direction during force field computations. The frozen-core approximation was used for both MP2 and CCSD(T) calculations. B3LYP calculations were run using an "army grade" pruned grid (JANS=2) with roughly 71,000 grid points and 155 radial points per atom. In general calculations were run on four processors, with only a few exceptions. Larger basis sets such cc-pV5Z, cc-pV6Z and Pcseg-4 were run with fully numerical gradients for all theory methods, and required calculations to be run on a single processor. For a few special cases, the anharmonic frequencies were calculated using Vibrational self-consistent field (VSCF) [48]. These calculations were run with third order mode couplings, ten grid points, with second order degenerate perturbation theory corrections. Extrapolated complete basis set (CBS) limits were obtained using the following two-point linear extrapolation formula generated by Halkier et al [49].
Where X and Y are the basis set's zeta value, with Y = X − 1, CBS limits were found using the highest zeta basis sets calculated per basis "family". Error analysis was done by taking the average error between the calculated values and either experimental data or CBS limits.

Calculations Removed from Study
After careful consideration, certain combinations of basis set and theory methods were dropped from the study due to issues discussed in this section. The version of Gamess(US) used for this study, does not support basis sets with high angular momentum i-functions (5/6 zeta basis sets) to be utilizied for predicting IR intensities or Raman activities. Zeta 5/6 Hessian calculations for CCSD(T) were too computationally intensive and did not provide significant improvements on molecular property predictions to justify their extreme runtimes. Sadlej-LPolX-fl was also

Molecular Geometries
The floating integer values within Figures 1-4 throughout the following sections are the zeta values for the basis sets. As expected, CCSD(T) performs best for calculating molecular geometries yielding consistent results with few errors. B3LYP converges faster than MP2 and CCSD(T), with the latter providing only slightly more accurate bond lengths and angles, for water.

Figure 1: Unsigned Error of CCSD(T) Bond Angle Calculations Compared to Experimental Data
It is notable to mention that the Pople-style basis sets provide reasonably accurate geometries, and it requires considerably high zeta basis sets to converge on the accurate geometry; particularly for the bond angle of water. While bond length converges quickly, for all levels of theory studied here, the bond angle converges more slowly, as shown for CCSD(T) in Figure 1. In general, most basis families show linear trend lines of slower convergence as either diffuseness or zeta values increase; making DFT with pople basis sets particularly well-suited for optimizing the geometries of large molecules. The Def2 family performs well in both runtime and accuracy, similarly both the PCseg and APCseg yield high accuracy predictions with lower runtimes, particularly when compared to their CC counterparts.

Vibrational Frequencies
When comparing vibrational frequencies to experimental data CCSD(T) calculations yielded the most accurate results, as expected. Interestingly, the NO7 and SNS basis set families yield results within ten wavenumbers (here, we consider this a reasonable threshold for the characterization or discrimination of most molecular species) with minimal runtimes [11]. The Def2 basis family consistently outperforms both correlation consistent and polarization consistent basis sets, at significantly reduced computational expense. Both the augmented and non-augmented SPK-basis sets also perform very well in calculating frequencies, producing similar accuracy and runtimes as Def2 sets. While, for B3LYP, the commonly used Pople family, converges above the 10 cm −1 threshold with similar runtimes. It is noteworthy to mention that fortuitous cancellation of errors often plays a large role here, particularly with smaller basis sets. However, if these errors consistently off-set a deficiency in the level of theory encountered at the basis set limit, or that is not overcome unless much larger basis sets are used, this may prove useful. As an example, for the prediction of vibrational frequencies, the 'may' calendar basis sets consistently out-perform even their fully augmented counterparts at only a minor additional computational expense compared to the CCn basis at all levels of theory studied. B3LYP is known to over-estimate bond-lengths, and under-estimate harmonic frequencies. Often, a scaling factor is applied to reduce the errors in predicting anhamonic values from calculations performed under the harmonic approximation [51][52][53]. In Table 4 we can see clearly that the anharmonic frequencies for B3YLP are all below the experimental values. Again in Table 4, the anharmonic frequencies are shown for a few select levels of theory and basis sets. Here, we can see that B3LYP consistently underestimates the experimental frequencies, and that even the may-cc-pVTZ basis set values lie within 5 cm −1 of the cc-pV6Z values. Similarly, we see for MP2 that the def2, may and pc families provide excellent agreement with the cc-pV5Z values. For each level of theory, the def2-QZVPD basis set approaches the accuracy of cc-pV5Z/cc-pV6Z much faster than most other basis sets of similar quality. Since anharmonic calculations are computationally expensive, CCSD(T) calculations were only performed using two of the def2 basis sets, where it could be seen that although an improvement over MP2 is observed, there is still not overall agreement with the experimental harmonic frequencies. This is likely because MP2 and CCSD(T) do not perform as well when used to calculate energetics at geometries far from the equilibrium position. To emphasize this point, calculations were additionally ran using the SCS-MP2 and CCSD(2) T levels of theory, which in each case show a marked improvement in the agreement with experimental anharmonic frequencies. All of the levels of theory considered here, struggle to reproduce the ν 2 /ω 2 frequency. We note that the CCSD(2) T /def2-QVZPD level of theory only has a MUE of ≈13 cm −1 .

Infrared Intensities
From Figure 3, we can see that though the Pople and other small basis sets can outperform larger basis sets at predicting some infrared intensities, however, caution is urged since this is not always the case. Both the B3LYP and MP2 levels of theory have significant errors when calculating infrared intensities under the harmonic approximation, even when extrapolated to their CBS limits. The B3LYP theory level convergence line appears to lie just above the ten percent error margin, some basis sets -particularly small basis sets, such as Pople -deviate from this producing more accurate results. In some cases where small basis sets are utilized, the errors arising actually compensates for the deficiencies in B3LYP, producing seemingly more accurate results. However, these results may occur sporadically, and caution is urged using these results, which could easily be taken out of context [54]. Both the harmonic and anharmonic calculation of the infrared frequencies is subject to the accuracy of the initial Hessian used to project the geometries, as well as the accuracy of the level of theory and basis set used to calculate the energetics and dipole moment. As explained previously, the MP2 and CCSD(T) levels of theory are not as reliable at geometries required for anharmonic, or even harmonic calculations. As a result, infrared intensities are not as dependent upon diffuse and polar functions as Raman activities, but rather how accurately the relative energies are calculated at each geometry, which explains why some basis sets that do not incorporate large polarization and diffuse functions are only necessary to the extent that the energetics and dipole moment are calculated with sufficient accuracy. Table 4 shows that by utilizing SCS-MP2 [55] or CCSD(2) T [56] levels of theory (the former at no additional cost, the latter at twice the cost of regular CCSD(T) calculations), further improvements in the frequencies and IR intensities can be obtained.

Raman Activities
In contrast to the IR intensities, the Raman activities can be very accurately calculated under the harmonic approximation, in particular when the MP2 level of theory is used (in agreement with previous studies [3,4]) and are extremely sensitive to the degree of polarization and diffusivity incorporated into the basis sets. It is clear from Figure 4 that while increasing the zeta-level or overall size of basis set is helpful, those basis sets that incorporate higher polarization levels and more diffuse functions are generally more accurate at reproducing Raman activities. In particular, the spectroscopic basis sets all perform well here (e.g., Sadlej, N07, SNS). The (PPD versions of) Def2 and (aug) Roos basis sets perform very well here also. The NLO-1 basis set also performs fairly well in calculating Raman activities, comparable to many quadruple-zeta basis sets. The dependence on the level of diffusivity is well demonstrated from the correlation consistent and calendar basis sets (cc→ may→ jun→ jul→ aug).

Conclusion
From our investigation we have identified various trends for the predictions of water's molecular properties, of which the most notable is the Def2-n basis set family's overall performance. Poplestyle basis sets tend to be very consistent with their predictions, typically producing minimal changes between each individual basis set's calculation. Pople basis sets are a reasonable choice for exploration, but there is little systematic improvement within their family, and are not recommended when results require high accuracy or quantitative values. Although they are occasionally somewhat accurate for infrared intensities when used with the B3LYP level of theory, they are not recommended to predict Raman activities. The Aug-CCn family is good for predicting Raman activities, particularly with the MP2 level of theory, however, they are computationally heavy, and outperformed by less computationally heavy options, such as aug-PCn, Def2-PD, and SPK augmented versions, as well as the majority of spectroscopic basis sets in this area. These basis sets have been shown to produce results comparable to six zeta basis sets, with minimal runtimes for the calculation of many optical properties. Another notable trend is the NO7 and SNS basis families success with calculating vibrational frequencies with harmonic approximation methods. Regarding the ANO series of basis sets, the NASA and Neese series offer high accuracy for geometries, as well as vibrational frequencies and infrared intensities, but lack sufficient diffuse functions to be able to accurately determine Raman activities (even the Aug-Neese sets). Only the Roos augmented ANO basis sets were able to accurately reproduce Raman activities, and represent an excellent all-round choice. We noticed that for the calculation of Raman activities basis sets which incorporated large numbers of diffuse functions always provided superior accuracy, however our scope was limited to the water molecule and further  11 investigations need to be conducted.We are currently increasing the scope of this study to include various other molecular species in the hopes of validating trends found in this study.