Uncertainty Quantification in First-Principles Predictions of Harmonic Vibrational Frequencies of Molecules and Molecular Complexes

TL;DR

This paper introduces a computationally efficient method to quantify uncertainty in density functional theory predictions of molecular vibrational frequencies, improving confidence in spectroscopic and thermodynamic analyses.

Contribution

The authors develop a novel, low-cost uncertainty estimation technique for DFT vibrational frequencies using BEEF-vdW functional's error bounds, applicable to various molecular systems.

Findings

Uncertainty estimates reliably bound different exchange-correlation functionals.

Modes with bending, torsional, or non-covalent interactions show higher uncertainty.

The method is simple, efficient, and suitable for routine vibrational frequency analysis.

Abstract

Accurate prediction of molecular vibrational frequencies is important to identify spectroscopic signatures and reaction thermodynamics. In this work, we develop a method to quantify uncertainty associated with density functional theory predicted harmonic vibration frequencies utilizing the built-in error estimation capabilities of the BEEF-vdW exchange-correlation functional. The method is computationally efficiency by estimating the uncertainty at nearly the same computational cost as a single vibrational frequency calculation. We demonstrate the utility and robustness of the method by showing that the uncertainty estimates bounds the self-consistent calculations of six exchange correlation functionals for small molecules, rare gas dimers, and molecular complexes from the S22 dataset. Ten rare-gas dimers and the S22 dataset of molecular complexes provide a rigorous test as they are systems with complicated vibrational motion and non-covalent interactions. Using coefficient of variation as a uncertainty metric, we find that modes involving bending or torsional motion and those dominated by non-covalent interactions are found to have higher uncertainty in their predicted frequencies than covalent stretching modes. Given the simplicity of the method, we believe that this method can be easily adopted and should form a routine part of DFT-predicted harmonic frequency analysis.

Tables

Table 1. Table 1 : Harmonic vibrational frequencies for eight small molecules in meV as predicted with the PBE, RPBE, PBEsol, PW91, optPBE-vdW, and BEEF-vdW XC functional. Experimental values and the BEEF-vdW ensemble statistics are also listed. Standard deviation values are in meV. COV, skew, and kurtosis are dimensionless.

	Stretch	Bend	Tors.	Expt.	PBE	RPBE	PBEsol	PW91	optPBE	BEEF	$\mu$	$\sigma$	COV	Skew	Kurt.	JB
									-vdW	-vdW
H2	1	0	0	545.553	536.1	540.4	527.7	538.6	543.5	555.0	555.3	7.5	0.01	-0.03	-0.07	0.7
CH4	0	1	0	168.2 54	159.0	160.1	156.2	159.6	164.6	164.5	163.7	10.1	0.06	-0.17	0.20	13.0
	0	1	0	194.3 54	186.8	187.5	184.5	187.4	191.6	191.6	191.1	7.6	0.04	-0.10	0.10	4.2
	1	0	0	376.5 54	368.5	367.4	367.0	369.4	372.0	372.0	372.2	6.2	0.02	-0.06	0.10	2.0
	1	0	0	391.5 54	382.6	381.1	382.2	383.2	385.1	385.1	385.4	7.9	0.02	-0.09	0.08	3.2
NH3	0	1	0	126.7 54	125.8	128.7	122.4	125.4	127.9	130.6	129.5	13.2	0.10	-0.38	0.72	91.3
	0	1	0	209.6 54	200.3	201.5	197.7	200.8	202.7	205.8	205.0	8.7	0.04	-0.13	0.17	8.0
	1	0	0	434.7 54	420.5	418.9	420.4	421.5	418.6	425.2	425.4	7.4	0.02	-0.06	0.04	1.3
	1	0	0	443.5 54	435.7	433.8	436.3	436.6	433.0	440.1	440.4	8.5	0.02	-0.08	0.06	2.4
H2O	0	1	0	204.3 54	197.1	198.7	194.5	197.1	199.1	202.3	201.8	9.4	0.09	-0.32	0.47	52.5
	1	0	0	475.1 54	461.1	460.0	461.7	462.2	458.5	467.1	467.4	8.2	0.02	-0.08	0	2.1
	1	0	0	488.8 54	473.7	472.5	474.7	474.9	470.9	479.7	480.0	8.6	0.03	-0.11	0.04	4.2
CO	1	0	0	269.053	265.3	262.3	267.2	265.8	264.3	265.4	265.5	5.7	0.02	-0.08	0.05	2.3
N2	1	0	0	292.353	292.4	289.4	294.1	292.9	291.0	293.6	293.8	6.0	0.02	-0.08	0.01	2.1
H2CO	0.04	0.96	0	147.7 54	142.4	142.3	141.4	142.8	142.8	144.8	144.1	9.4	0.07	-0.17	0.22	13.7
	0	1	0	159.6 54	151.3	151.2	150.3	151.8	152.3	154.1	153.7	7.3	0.05	-0.09	0.09	3.4
	0.10	0.90	0	193.8 54	183.5	183.8	181.7	184.1	185.0	187.3	186.3	5.5	0.03	-0.70	0.53	186.7
	0.90	0.10	0	218.7 54	219.0	217.0	221.1	219.2	217.5	219.4	220.0	3.4	0.02	0.55	0.28	107.4
	1	0	0	365.0 54	346.5	345.9	344.6	347.5	346.1	351.6	351.8	7.3	0.02	-0.08	0.07	2.5
	1	0	0	373.0 54	352.3	351.4	350.6	353.3	351.7	357.5	357.8	7.5	0.02	-0.09	0.05	2.9
CO2	0	1	0	82.7 55	79.7	79.1	80.2	79.7	78.7	79.6	79.0	7.3	0.09	-0.32	0.47	52.4
	1	0	0	165.3 55	164.4	162.6	165.9	164.5	163.3	164.5	164.6	3.6	0.02	-0.08	0	2.1
	1	0	0	291.255	292.7	289.5	296.0	292.8	290.0	292.4	292.6	7.7	0.03	-0.11	0.04	4.2
ME					7.8	8.3	8.6	7.2	7.3	3.8
MAE					8.0	8.5	9.5	7.4	7.4	5.3

Table 2. Table 2 : Harmonic vibrational frequencies in meV as predicted with the PBE, RPBE, PBEsol, PW91, optPBE-vdW, and BEEF-vdW XC functionals for C 2 H 6 and HCOOH. Experimental values and the BEEF-vdW ensemble statistics are also listed. Standard deviation values are in meV. COV, skew, and kurtosis are dimensionless.

	Stretch	Bend	Tors.	Expt.	PBE	RPBE	PBEsol	PW91	optPBE	BEEF	$\mu$	$\sigma$	COV	Skew	Kurt.	JB
									-vdW	-vdW
C2H6	0	0	1	37.6 54	36.7	36.5	37.1	36.7	37.6	37.8	34.6	23.8	0.69	-0.14	-1.12	111.1
	0	0.73	0.27	101.9 54	99.1	99.4	98.0	99.3	100.1	101.3	99.5	11.5	0.12	-0.72	1.03	261.2
	0.97	0.03	0	126.0 54	122.8	120.6	125.0	122.6	120.7	122.0	122.3	2.7	0.02	0.16	-0.45	25.4
	0	0.64	0	154.5 54	145.9	146.1	144.3	146.3	147.6	149.4	148.8	8.1	0.05	-0.11	0.12	5.2
	0	1	0	178.3 54	168.2	168.8	165.9	168.9	170.8	172.9	172.3	9.3	0.05	-0.15	0.16	9.6
	0.03	0.97	0	179.6 54	168.7	170.1	167.9	170.1	171.6	173.9	173.4	8.8	0.05	-0.07	0.01	1.6
	0	0.71	0.29	189.2 54	180.0	180.7	177.6	180.6	182.2	184.8	183.6	7.9	0.04	-0.12	0.17	7.2
	0	0.62	0.38	192.4 54	180.0	180.8	177.6	180.7	182.3	184.8	184.2	8.1	0.04	-0.09	0.07	3.1
	1	0	0	377.3 54	367.7	366.8	366.4	368.6	366.1	371.2	371.4	6.7	0.02	-0.08	0.09	2.8
	1	0	0	379.5 54	368.0	367.0	366.7	368.9	366.4	371.4	371.6	6.6	0.02	-0.07	0.10	2.5
	1	0	0	389.3 54	377.1	375.9	376.5	377.8	374.6	380.3	377.6	8.0	0.02	-0.09	0.07	2.5
	1	0	0	393.6 54	377.5	376.3	376.8	378.1	374.9	380.4	380.5	7.9	0.02	-0.09	0.08	3.2
HCOOH	0.10	0.90	0	77.5 55	74.8	74.6	74.8	75.0	75.1	76.1	73.0	9.1	0.12	-2.74	10.91	1.2e4
	0	0	1	79.1 55	83.5	82.4	84.6	83.4	82.2	83.2	84.0	12.2	0.14	0.48	0.41	90.8
	0	0	1	128.1 55	124.2	123.9	123.9	124.5	123.9	125.7	124.0	8.0	0.06	-0.83	0.47	248.0
	0.73	0.27	0	137.0 55	133.4	131.6	135.4	133.6	132.0	133.6	133.8	3.5	0.03	0.99	5.63	2968.1
	0.20	0.80	0	152.4 55	155.5	155.3	155.8	155.8	155.9	158.2	158.7	5.6	0.04	0.67	0.17	152.0
	0.07	0.93	0	172.0 55	167.4	167.9	165.6	167.9	168.7	171.3	170.8	6.5	0.04	-0.20	-0.23	17.7
	0.91	0.09	0	219.4 55	218.2	216.2	220.5	218.4	216.6	218.9	219.3	4.7	0.02	0.21	-0.08	15.2
	1	0	0	364.9 55	367.6	366.7	365.7	368.6	367.3	373.0	373.2	7.2	0.02	-0.09	0.08	3.2
	1	0	0	442.6 55	448.1	448.6	447.3	449.1	446.4	455.4	455.7	8.4	0.02	-0.07	0.06	1.9
ME					5.1	5.5	5.7	4.6	5.2	2.2
MAE					6.6	6.9	7.1	6.3	6.4	5.2

Table 3. Table 3 : Harmonic vibrational frequencies for rare-gas dimers in meV as predicted with the PBE, RPBE, optPBE-vdW, vdW-DF2, and BEEF-vdW XC functionals. BEEF-vdW ensemble statistics are also presented. Standard deviation values are in meV. COV, skew, and kurtosis are dimensionless. Experimental values are from Ogilvie and Wang 58 , 59 . The PBE and TPSSh values from Tao and Perdew 13 include BSSE corrections.

	Expt.	PBE13	PW9112	TPSSh13	PBE	RPBE	optPBE	vdW	BEEF	$\mu$	$\sigma$	COV	Skew	Kurt.	JB
							-vdW	-DF2	-vdW
He2	4.116	8.936	12.4	4.342	8.4	10.0	10.3	9.1	15.6	15.1	4.7	0.31	0.76	0.20	195.9
HeNe	4.334	6.539	10.0	6.609	6.9	8.3	9.2	9.5	13.1	13.1	5.9	0.45	0	-0.68	38.5
HeAr	4.317	6.220	8.1	4.896	7.6	6.7	8.5	7.5	11.0	10.9	5.6	0.52	-0.03	-0.86	61.9
HeKr	3.972	4.842	7.6	5.606	7.0	6.1	8.1	6.9	11.2	10.9	5.5	0.51	-0.13	-0.86	67.3
Ne2	3.536	3.942	6.2	4.066	5.6	7.0	6.3	5.1	8.7	8.5	4.4	0.52	-0.07	-0.86	63.3
NeAr	3.494	3.681	4.6	2.311	3.9	5.2	4.5	3.6	4.2	5.0	3.1	0.61	0.32	-0.77	83.5
NeKr	3.243	3.335	4.1	2.834	4.9	4.1	5.5	4.9	5.5	5.4	3.0	0.55	0	-0.89	66.0
Ar2	3.837	3.086	3.2	2.534	2.8	4.5	5.8	5.3	6.4	6.0	2.3	0.38	-0.73	0.16	179.8
ArKr	3.465	2.420	2.9	2.289	3.3	2.7	3.9	3.8	4.6	4.3	1.9	0.43	-0.40	-0.47	71.7
Kr2	2.923	1.808	2.1	2.120	2.2	2.2	3.7	3.5	4.2	4.0	1.4	0.36	-0.71	0.29	175.0
ME		-0.757	-2.4	-0.037	-1.5	-2.0	-2.9	-2.2	-4.7
MAE		1.339	2.8	1.012	1.9	2.3	2.9	2.2	4.7

Table 4. Table 4 : Water dimer vibrational frequencies (meV) predicted with various XC functionals. Modes 1 through 6 are intermolecular, while modes 7 through 12 are intramolecular. BLYP and B3LYP values were predicted with aug-cc-pVTZ(-f) basis sets 62 , and MP2 values with the aug-cc-pVDZ basis set 67 . MP2 values were multiplied by a scale factor by Dunn et al. and are rescaled here. All experimental intermolecular modes except mode 4 are not corrected for anharmonic effects. ME and MAE are calculated for intramolecular frequencies only.

Mode	Stretch	Bend	Tors	Expt.	PBE	RPBE	PBEsol	PW91	optPBE	BLYP62	B3LYP62	MP267	BEEF	$\mu$	$\sigma$	COV	Skew	Kurt.	JB
									-vdW				-vdW
1				10.971	17.0	13.2	17.2	14.2	13.4	16.1	16.5	16.0	13.5	26.3	25.2	1.22	0.49	-1	164.7
2				12.871	19.7	15.3	21.3	18.8	17.4	19.3	19.7	18.5	25.3	23.3	25.3	1.10	0.73	-0.77	228.9
3				13.471	21.1	18.5	22.3	20.0	17.8	20.2	19.8	18.8	35.7	26.1	25.6	0.98	0.63	-0.90	199.7
4				18.672	23.4	19.2	27.0	23.4	20.9	23.8	23.9	22.9	39.8	31.1	23.6	0.76	0.64	-0.79	188.6
5				38.673	45.9	38.4	49.0	44.0	42.1	46.1	46.5	44.6	47.7	40.3	22.6	0.56	0.51	-0.90	154.6
6				64.873	79.5	69.1	88.2	80.2	72.2	77.0	78.6	79.8	74.6	64.8	21.1	0.33	-0.32	-0.47	52.5
7	0	1	0	204.969	196.9	198.7	193.9	196.8	199.0	198.9	203.0	201.5	204.6	201.6	9.2	0.05	-0.15	0.19	10.5
8	0	1	0	206.969	200.0	201.1	198.2	200.4	201.0	201.0	205.3	203.8	206.9	203.7	9.1	0.04	-0.16	0.17	10.9
9	1	0	0	460.969	440.4	447.8	431.6	440.3	444.6	437.1	456.0	459.3	456.5	457.0	9.2	0.02	-0.06	-0.03	1.4
10	1	0	0	470.769	460.3	459.9	460.7	461.7	457.9	452.9	470.4	470.9	465.6	466.6	8.2	0.02	-0.06	0.03	1.2
11	1	0	0	481.269	470.0	469.6	470.6	471.3	467.8	462.3	480.2	484.3	475.3	476.9	8.8	0.02	-0.06	0.03	1.3
12	1	0	0	483.469	472.5	471.9	473.2	473.8	469.8	464.7	482.1	486.8	476.6	478.7	8.6	0.02	-0.07	0.04	1.6
ME					11.3	9.8	13.3	10.6	11.3	15.2	1.8	0.2	3.8
MAE					11.3	9.8	13.3	10.6	11.3	15.2	1.8	2.5	3.8

Equations

$$\begin{array}{rl}{m}_i{\ddot{x}}_i^{\alpha }& =F_i^{\alpha }\\ & =-\sum _{\beta }\sum _j\left(\frac{{\partial }^{2}E}{\partial x_i^{\alpha }\partial x_j^{\beta }}\right)x_j^{\beta }\\ & =-\sum _{\beta }\sum _j{\Phi }_{ij}^{\alpha \beta }{x}_j^{\beta },\end{array}$$

$${\omega }^2\mathbf{e}=\mathbf{He}.$$

$${\text{H}}_{3\ast (i-1)+\alpha ,3\ast (j-1)+\beta }=\frac{1}{\sqrt{m_i{m}_j}}{\Phi }_{ij}^{\alpha \beta }.$$

$$\Phi_{ij}^{\alpha\beta}\approx\begin{cases}\frac{1}{4h^{2}}\bigg{(}E_{i+h,j+h}+E_{i-h,j-h}-E_{i+h,j-h}-E_{i-h,j+h}\bigg{)},&\text{if}\ i\neq j\\ -\sum_{k\neq i}\Phi_{ik}^{\alpha\beta},&\text{if}\ i=j.\\ \end{cases}$$

$$E=E_{KE}+E_{ions}+E_{e-e}+E_{xc},$$

$$E_{xc}=\sum _{m=0}^{29}{a}_m{E}_m^{\text{GGA}-\text{x}}+{\alpha }_c{E}^{\text{LDA}-\text{c}}+(1-{\alpha }_c)E^{\text{PBE}-\text{c}}+E^{\text{nl}-\text{c}},$$

$$E_m^{\text{GGA}-\text{x}}=\int {ϵ}_{\text{x}}^{\text{UEG}}[n(\mathbf{r})]B_m\left\{s[n(\mathbf{r}),\nabla n(\mathbf{r})]\right\}d\mathbf{r}.$$

Full text

\SectionNumbersOn

CMU]Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA

\abbreviationsIR,NMR,UV

1 Introduction

Accurate prediction of molecular vibrational frequencies by ab initio methods is important in many areas of chemistry and physics1, 2. Calculating the enthalpy and free energy of a reaction, for example, requires zero-point energy (ZPE) and finite temperature contributions, both of which depend on vibrational frequencies. The accuracy of the frequencies depends on the level of theory used to model the system 3, convergence criteria 4, and if anharmonic corrections have been included5. A widely-used method to predict vibrational frequencies is density functional theory (DFT) 6. The accuracy of a DFT calculation depends strongly on the choice of the exchange-correlation (XC) functional. There has been much work devoted to identifying the best XC functional and DFT calculation parameters for accurately predicting the frequencies of small-to medium-sized molecules. 3, 7, 8, 9, 10, 11 Effort has also been extended to predicting frequencies of weakly-bonded systems.12, 13 Patton and Pederson 12 and Tao and Perdew 13 showed that DFT calculations of rare-gas dimers at the generalized-gradient (GGA) 14 level correct the overbinding introduced by the local density approximation (LDA) functional. This body of results has shown that DFT-predicted frequencies at the GGA level may be accurate to within tens of meV of experimental frequencies for many small molecules. However, for complex vibration modes, the sensitivity of the predictions to the choice of XC functionals remains unclear.

One naïve way to estimate the uncertainty associated with DFT-predicted frequencies is to perform the calculation with multiple XC functionals. The selection of functionals is somewhat arbitrary, however, and the calculation must be performed multiple times making it computationally inefficient. The Bayesian error estimation functional with van der Waals correlation (BEEF-vdW) 15 is a GGA-level XC functional that can systematically estimate uncertainty in DFT predictions. It possesses built-in uncertainty estimation capabilities in the form of an ensemble of GGA XC functionals that are calibrated to reproduce the error observed between experimental measurements and DFT predictions. Using the BEEF-vdW ensemble is computationally efficient compared to performing many calculations using different XC functionals, as results for thousands of XC functionals are obtained non-self-consistently by using one self-consistent calculation. BEEF-vdW has been applied to quantify uncertainty in magnetic ground states 16 heterogeneous catalysis 17, 18, 19, electrocatalysis 20, 21, 22, and mechanical properties of solid electrolytes 23. The BEEF-vdW model space includes contributions from the non-local vdW-DF2 XC functional 24. It therefore offers potential improvements compared to other GGA-level functionals in describing van der Waals forces, which are traditionally not well described by DFT 25.

In this work, we present a computationally efficient method to estimate the uncertainty of harmonic vibrational frequencies of non-periodic systems such as molecules and molecular complexes. The harmonic approximation is valid at temperatures well below the dissociation limit 26, which are on the order of $1,000$ K for most systems considered herein but can be as low as $370$ K for the rare-gas dimers 27, 28. Frequency ensembles are calculated by solving an eigenvalue problem for an ensemble of Hessian matrices that depend on the second derivatives of the energy of the system with respect to its nuclear coordinates. The formulation of the eigenvalue problem is presented in Section 2.1. The ensemble of Hessians is calculated from an ensemble of energies determined with BEEF-vdW, which is described in Section 2.2. Computational details are provided in Section 2.3.

We first apply our method to a set of small benchmark molecules in Section 3.1 to determine if it can capture the uncertainty associated with choosing an XC functional. We also make comparisons to experimental measurements. We then consider a case study of the ten rare-gas dimers studied by Patton and Pederson 12 and Tao and Perdew 13 in Section 3.2. Rare-gas dimers offer a simple test of the ability of the BEEF-vdW XC functional and the BEEF-vdW ensemble to describe weakly-bonded systems. In Section 3.3, we study the non-covalently bonded complexes in the S22 dataset 29. Because some XC functionals are known to perform well on rare-gas dimers but poorly on larger complexes (and vice versa) 30, 29, the S22 dataset offers an extensive test of the ability of BEEF-vdW to describe complicated vibrations in non-covalently bonded systems. We compare the predictions for all the molecules and molecular complexes in Section 3.4. In considering all the results, we utilize coefficient of variation as a metric to quantify uncertainty and should that the coefficient of variation is largest for modes that involve bending or torsion or whose bonds are non-covalent. Localized stretching modes tend to have lower uncertainty in their predicted frequencies. The spread of the frequencies can also be used to quantify the uncertainty associated with quantities derived from the frequencies (e.g., Gibbs free energy or ZPE). Given the simplicity of the developed method, we believe that this should form a routine part of DFT-predicted harmonic frequency analysis.

2 Methods

2.1 Harmonic vibrational frequencies

In this section, we discuss a method for predicting harmonic vibrational frequencies by computing energies using DFT. Let $N$ be the number of atoms in a molecule, and let $i$ and $j$ index the atoms so that $1\leq i,j\leq N$. $\alpha$ and $\beta$ denote the Cartesian directions [i.e., $\alpha,\beta=x(1),y(2),z(3)$]. We make the harmonic approximation, in which a Taylor series expansion of the system’s potential energy, $E$, is truncated after the second-order term 31, 32, 33. In this approximation, the force on atom $i$ in the $\alpha$-direction, $F_{i}^{\alpha}$, is proportional to the displacement of every other atom in the system from its equilibrium position. Newton’s 2nd law for atom $i$ in the $\alpha$-direction is thus

[TABLE]

where $m_{i}$ is the mass of atom $i$. The second derivatives in eq 1 are called the harmonic force constants, $\Phi_{ij}^{\alpha\beta}$. By assuming temporally periodic atomic motions, the set of equations represented by eq 1 for all atoms in all directions can be converted into the eigenvalue problem

[TABLE]

Here, $\omega$ is a vibrational frequency, $\mathbf{e}$ is a vector of length $3N$ describing the motion of the atoms, i.e. the mode shape, and $\mathbf{H}$ is the 3$N\times 3N$ Hessian matrix. The elements of the Hessian are related to the harmonic force constants by

[TABLE]

The harmonic force constants are calculated by numerically approximating the second derivative in eq 1 using a central finite difference of the energies with respect to perturbations of the equilibrium structure. Using the shorthand $E(x_{i}^{\alpha}\pm h,x_{j}^{\beta}\pm h)=E_{i\pm h,j\pm h}$, where $h$ is the perturbation magnitude of the atomic displacement, the central difference formulas are 34

[TABLE]

Note that we do not directly calculate the harmonic force constants where $i=j$. This practice reduces numerical error associated with the eggbox effect by enforcing conservation of momentum (i.e., translational invariance) 35. The harmonic force constants are symmetrical with respect to permutation of the atomic and direction indices 36 (i.e., $\Phi_{ij}^{\alpha\beta}=\Phi_{ji}^{\beta\alpha}$), so that we only need to calculate the upper triangular portion of the Hessian. Other molecule-specific symmetries can further reduce the number of force constants to be calculated, but we do not consider these here.

The harmonic force constants can also be calculated with a finite difference of the force on atom $i$ 32. We use the energies because, as will be explained in Section 2.2, BEEF-vdW provides uncertainty estimation in the system energy and not in the atomic forces.

2.2 Bayesian error estimation

The energy of a system can be predicted using DFT and takes the form 6

[TABLE]

where $E$ is the electronic energy of the system used in eq 4. $E_{KE}$ is the kinetic energy, $E_{ions}$ is the potential energy of the electrons due to Coulombic interactions with the ions, and $E_{e-e}$ is the potential energy of the electrons due to Coulombic interactions with other electrons, all of which can be calculated exactly. $E_{xc}$ is the exchange-correlation energy, whose value is approximated by the chosen XC functional.

BEEF-vdW is an XC functional at the GGA level15. The BEEF-vdW model space takes the form

[TABLE]

where $a_{m}$ and $\alpha_{c}$ are multiplicative factors, $E^{\text{LDA}-\text{c}}$ is a correlation contribution from the local Perdew-Wang LDA correlation 37, $E^{\text{PBE}-\text{c}}$ is a correlation contribution from the PBE semi-local correlation 38, $E^{\text{nl}-\text{c}}$ is a correlation contribution from the vdW-DF2 non-local correlation 24, and $E_{m}^{\text{GGA}-\text{x}}$ is the contribution to the exchange energy given by

[TABLE]

In eq 7, $n(\mathbf{r})$ and $\nabla n(\mathbf{r})$ are the electron density and its gradient, $s$ is a function that parameterizes $n$ and $\nabla n$, $\epsilon_{\text{x}}^{\text{UEG}}$ is the exchange energy density of the uniform electron gas, and $B_{m}$ is the $m$th Legendre polynomial. To determine the optimal BEEF-vdW XC functional, Wellendorff et al. fit $a_{m}$ and $\alpha_{c}$ to energetic and structural data describing bonding in a variety of chemical and condensed matter systems. These parameters were regularized to prevent overfitting to the training data.

BEEF-vdW provides a systematic approach to estimating uncertainty in a DFT energy calculation by using an ensemble of XC functionals around the optimal BEEF-vdW XC functional. A self-consistent DFT calculation is first performed using the optimal parameters, yielding a converged electron density. This density is then used with distributions of $a_{m}$ and $\alpha_{c}$ to non-self-consistently generate an ensemble of energies using eq 6. The distributions of $a_{m}$ and $\alpha_{c}$ are tuned such that the spread of the ensemble energies reproduces the errors observed when comparing the experimental training data to BEEF-vdW self-consistent predictions using the optimal XC functional.

BEEF-vdW provides an ensemble of energies rather than a single energy for a DFT calculation. This ensemble can be propagated in eq 4 to obtain an ensemble of numerical derivatives $\partial^{2}E/\partial x_{i}^{\alpha}x_{j}^{\beta}$ and Hessians. The eigenvalue problem can be solved for each Hessian in the ensemble to determine the ensemble of frequencies.

2.3 Computational details

Self-consistent DFT calculations were performed with the real-space projector-augmented wave method 39, 40 as implemented in GPAW 41, 42. The BEEF-vdW XC functional was used with 2000 ensemble functionals for each calculation. Using more than 2000 functionals has been found to have little effect on the standard deviation of the ensemble energy values 15, 23. We used a real-space grid spacing of $0.12$ Å for the rare-gas dimers and $0.18$ Å otherwise. Molecules were surrounded by vacuum in cubic boxes. Box lengths were determined so that adding eight additional real-space grid points along each axis changed the energy of the relaxed structure by no more than $0.01$ eV. This criterion resulted in a box length of at least $10$ Å for each system. Equilibrium geometries were determined by relaxing the structure so that each atom experienced a force of less than $0.01$ eV$/\mbox{\AA}$. Starting geometries for the S22 dataset were obtained from the Benchmark Energy and Geometry DataBase 43. All single-point calculations were converged so that the energy variation between the final three iterations was less than $10^{-8}$ eV.

To obtain perturbed energies for numerical estimation of the second derivative in eq 4, atomic displacements ($h$) of at most $\pm 0.01\text{ }\mbox{\AA}$ were applied. The $0.01\text{ }\mbox{\AA}$ displacement size led to variations in the harmonic force constants on the order of $0.001$ eV$/$Å2 and at most $0.0625$ eV$/$Å2 with respect to force constants calculated with displacements as small as $0.001$ Å. A smaller displacement leads to smaller energy variations, which require more stringent convergence criteria and longer computation time, but can yield more accurate numerical derivatives. Numerical variation in the force constants will impact the calculated frequencies, but this effect is suppressed by two factors: (i) the requirement that the force constants satisfy conservation of momentum (eq 4), and (ii) the low sensitivity of the Hessian to the eigenvalue problem because it is Hermitian44. A comparison of the effect of numerical uncertainty in the force constants to the BEEF-vdW uncertainty for two vibrational frequencies (one high and one low) of the benzene-ammonia complex, a member of the S22 dataset, is shown in Figure 1. The numerical uncertainty histograms were generated by adding draws from a Gaussian distribution $N(0,0.0625)$ eV$/$Å2 to the force constants calculated using the BEEF-vdW XC functional. This process was repeated $10,000$ times to ensure converged error estimates. For both frequencies, the spread due to the numerical uncertainty in the force constants is smaller than the ensemble spread, indicating that in considering the latter we may ignore the effects of the former.

For larger molecules and molecular complexes, the atomic motions described by the eigenvectors in eq 2 are often delocalized. To analyze the atomic motions in terms of the movement of local groups of atoms, they were transformed from Cartesian to internal coordinates. The internal coordinates were determined using the MolMod Python package 45. Bond length, bend angle, and dihedral angle internal coordinates, corresponding to stretching, bending, and torsional modes, were considered. The atomic motions were then mapped onto the internal coordinates by following the decomposition method outlined by Boatz and Gordon46. The components of stretch, bend, and torsion internal coordinate motion for each mode sums to unity within a numerical error of at most $0.01$. As such, the relative contribution of each type of internal coordinate motion can be compared across different modes. Intermolecular modes are excluded from this mapping procedure, as their motion does not correspond to motion of the internal coordinates.

3 Results

3.1 Benchmark molecules

To test the effectiveness of the proposed method, we analyze the estimated uncertainty due to choice of XC functional in predictions of vibrational frequencies of small molecules. We first examine a set of eight small molecules: H2, N2, CO, CO2, H2O, NH3, H2CO, each of which has only stretching and bending modes. We then compare the results to those for two larger molecules, HCOOH and C2H6, that have more complicated modes (e.g., a torsional mode in the case of C2H6). We designate these ten molecules as our “benchmark” set. Frequencies were predicted using the PBE 38, RPBE 47, PBEsol 48, PW91 49, and optPBE-vdW 50, and BEEF-vdW functionals self-consistently.

Results for the eight small molecules are presented in Table 1 and for HCOOH and C2H6 in Table 2. Experimental values and BEEF-vdW ensemble statistics are also shown. To test for normality in the BEEF-vdW ensemble, we use the skew and kurtosis to calculate the Jarque-Bera (JB) statistic 51, 52, given by JB $=\frac{n}{6}\left[S^{2}+\frac{K^{2}}{4}\right]$, where $n=2000$ is the sample size of the ensemble, $S$ is the sample skew, and $K$ is the sample kurtosis. Under the null hypothesis that the ensemble is Gaussian, JB is approximately described by a $\chi^{2}(2)$ distribution for large $n$. We choose to reject the null hypothesis at a $95\%$ confidence level ($\alpha=0.05$), corresponding to a $\chi^{2}$ value of approximately 6 (i.e., we label the ensemble as non-Gaussian if JB $>6$). We use the standard deviation ($\sigma$) and coefficient of variation, $\text{COV}=\sigma/\mu$, where $\mu$ is the mean of the distribution, of the ensembles as measures of uncertainty due to choice of XC functional. A low COV indicates that the functionals tend to agree in their frequency predictions, while a high COV indicates disagreement. In Section 3.4 we find that a COV of $0.25$ reasonably separates ensembles with low and high uncertainty.

As a representative example, consider the asymmetric stretching mode of NH3. The BEEF-vdW ensemble of frequencies, predictions using other functionals, and the experimental value ($443.5$ meV 54) are plotted in Figure 2. The internal coordinate decomposition of the mode in Table 1 indicates that it is a pure stretching mode. We report a BEEF-vdW value of $440.1$ meV and an ensemble standard deviation of $8.5$ meV. The ensemble bounds the experimental frequency to within one standard deviation, demonstrating that DFT can accurately predict this frequency. In addition, the ensemble bounds the frequencies predicted by the other XC functionals to within one standard deviation, which indicates that the ensemble can reproduce the predictions of other XC functionals. Based on the low COV of the BEEF-vdW ensemble ($0.02$), there is also little disagreement in ensemble functionals in predictin this frequency. The ensemble has low skew ($-0.08$) and kurtosis ($0.06$). These values yield a JB of $2.4$, so that the ensemble is Gaussian.

In the Supporting Information, we propagate the uncertainty estimates provided by the NH3 frequency ensembles for all NH3 modes to quantify uncertainty in NH3 ZPE. We also show how ensembles for N2, H2, and NH3 can be used to quantify uncertainty in predicting the vibrational entropy for the Haber process, 3H2+N${}_{2}\rightarrow$ NH3.

The data in Table 1 indicate that most of the 23 modes of the small molecules are pure stretching or pure bending. There are only a few modes with internal coordinate mixing, e.g., the third mode for H2CO, which is $0.90$ stretching and $0.10$ bending. In terms of mean absolute error (MAE) in comparison to the experimental measurements, BEEF-vdW is the best XC functional for predicting vibrational frequencies. The experimental frequencies for 16 of the 23 frequencies are bounded to within one standard deviation of the BEEF-vdW ensemble, with the largest deviation ($2.07\sigma$) coming from the sixth mode of H2CO. On average, the predictions via BEEF-vdW deviate $0.70\sigma$ from the experimental values.

The predictions of the other XC functionals are also generally bounded by the BEEF-vdW ensemble to within one standard deviation, with each XC functional predicting at most three frequencies outside of $\pm\sigma$. The average deviation is smallest for PW91 ($0.44\sigma$), followed by optPBE-vdW ($0.46\sigma$), PBE ($0.52\sigma$), RPBE ($0.62\sigma$), and PBEsol ($0.74\sigma$).

The largest COV for an ensemble is $0.10$, for the lowest frequency mode of NH3. The low COV values indicate agreement among the ensemble functionals in predicting the frequencies. Each ensemble has relatively low skew and kurtosis values, the highest values being a skew of $-0.70$ for the third H2CO mode and a kurtosis of $0.72$ for the first NH3 mode. Eight modes have JB $>6$ and are non-Gaussian. Each of these eight modes is a pure bend mode or a mix of bending and stretching.

The results presented in Table 1 indicate that for small molecules with simple vibrational motion, the BEEF-vdW ensemble can bound both the spread of predictions of other XC functionals and experimental frequencies. To test the robustness of the BEEF-vdW ensemble for molecules containing torsional vibrations, we next apply the method to HCOOH and C2H6 and the results are presented in Table 2. By the JB test, 12 of the 21 modes have non-Gaussian ensembles, none of which correspond to pure stretch modes.

Based on the MAE, BEEF-vdW yields the frequencies closest to the experimental measurements. For 15 of the 21 modes, the experimental frequencies are bounded to within one standard deviation of the BEEF-vdW prediction, with an average deviation of $0.74\sigma$. Amongst the XC functionals considered, PW91 again has the smallest deviation ($0.32\sigma$) and PBEsol again has the largest ($0.58\sigma$). These deviations are smaller in comparison to those for the eight smaller molecules and seem to indicate better agreement among XC functionals for predicting of these frequencies. The ensembles for several modes, however, indicate this is not always the case. The two lowest frequency modes for both HCOOH and C2H6 have ensemble COV higher than the highest value reported in Table 1, with a maximum of $0.69$ for the C2H6 torsional mode. This torsional mode also has the largest standard deviation, $23.8$ meV, of the modes presented in both tables. While the ensemble for this mode easily bounds the experimental and other DFT-predicted frequencies, its high COV indicates disagreement among GGA-level functionals in predicting this frequency. This mode is also noteworthy because some 5% of the ensemble functionals predict a imaginary frequencies. These these imaginary frequencies correspond to negative solutions to the eigenvalue problem in eq 2 and indicate that some ensemble XC functionals predict that this mode is not stable. These values are not included in calculation of any ensemble statistics. There is also a significant fraction of ensemble functionals, over 15%, which predict very small (nonzero but less than 2 meV) frequencies. These findings are in agreement with previous analysis showing the difficulty in describing the torsional energy landscape of C2H6 within GGA-level DFT56.

3.2 Rare-gas dimers

In Section 3.1, we determined that the proposed methodology can reasonably bound frequencies both from experiment and predictions using other XC functionals for a benchmark set of ten molecules. We next consider ten rare-gas dimers built from He, Ne, Ar, and Kr. These systems are diatomic molecules and therefore have only one stretching mode. Rare-gas dimers are bonded by dispersion interactions and have been studied to determine the applicability of DFT to simple van der Waals systems 57. Patton and Pederson 12 and Tao and Perdew 13 found that while GGA functionals corrected the overbinding tendency of the LDA functional, they overestimated the interaction strength when the outer electron shell consisted of $s$ electrons (as in He2) and underestimated interaction strength when the outer shell consisted of $p$ electrons (as in Ne2) 25. The rare-gas dimers offer an interesting case study for both the BEEF-vdW XC functional, which explicitly models non-local interactions, and the BEEF-vdW ensemble, as the vibrational frequencies are small enough that the uncertainty in the prediction approaches the magnitude of the frequency itself ($\sim 10$ meV).

We compare the BEEF-vdW ensemble of frequencies for each molecule to those predicted by the PBE, RPBE, optPBE-vdW 50, vdW-DF2 24, and BEEF-vdW XC functionals in Table 3. Experimental values are from Ogilvie and Wang 58, 59. We also include DFT values from Patton and Pedersen 12 and Tao and Perdew 13. The tendency of GGA functionals to correct the severe overestimation of frequencies by LDA is consistent with our PBE and RPBE frequencies. For example, Tao and Perdew reported a frequency of $122.10$ meV for He2 using LDA, while our PBE value of $8.4$ meV for the same quantity is much closer to the experimental value.

The three vdW functionals tested, optPBE-vdW, vdW-DF2, and BEEF-vdW, overestimate the frequency for each dimer. BEEF-vdW in particular performs poorly. For example, it predicts a frequency of $13.1$ meV for HeNe, which has an experimental value of $4.334$ meV, and overpredicts the experimental frequencies by an average of $4.7$ meV. The previously-observed tendency of GGA functionals to overestimate $s$-shell electron interaction energy and underestimate $p$-shell electron interaction energy 13, 60, 61 is also reproduced in our results. The exceptions are the vdW functionals, which consistently overestimate the vibrational frequency for every dimer. This tendency can also be observed in considering the binding energy for each dimer, which is calculated from the electronic energy difference between the dimer and its constituent monomers. A table of binding energies can be found in the Supporting Information. In considering the MAE in comparison to the experimental frequencies, the most accurate predictions come via Tao and Perdew using the meta-GGA level TPSSh functional ($1.012$ meV MAE). The best GGA-level functional is PBE ($1.339$ meV MAE for Tao and Perdew’s values, $1.9$ meV MAE in this work).

Statistics for the BEEF-vdW ensembles are presented on the right side of Table 3. By the JB test, most notably, none of the ensembles are Gaussian. With the exception of He2, the experimental value for each frequency is bounded to within two standard deviations of the BEEF-vdW result, with an average deviation of $1.12\sigma$. For each dimer, the majority of the frequencies in the ensemble overestimate the experimental value. The frequency predictions from the other tested XC functionals are also bounded by the ensemble to within two standard deviations, but agreement between BEEF-vdW and the other functionals is worse than for the benchmark molecules in Section 3.1. The smallest average deviation from the BEEF-vdW predictions came from optPBE-vdW ($0.44\sigma$), followed by vdW-DF2 ($0.60\sigma$), PW91 via Patton and Pederson ($0.73\sigma$), RPBE ($0.81\sigma$), PBE ($0.86\sigma$), PBE via Tao and Perdew ($1.08\sigma$), and TPSSh via Tao and Perdew ($1.25\sigma$). This higher spread is also present in the BEEF-vdW ensembles. The COV of every ensemble is higher than the COV of any frequency ensemble in Section 3.1, with the exception of the C2H6 torsional mode. The high COV values and consistent overprediction of the experimental frequencies reflect the difficulty that GGA-level XC functionals have in consistently and correctly describing the van der Waals interactions in rare-gas dimers.

3.3 S22

The rare-gas dimers from Section 3.2 may not be representative of non-covalent interactions in larger systems 30, 29. As such, we next apply our methodology to S22, a set of 22 non-covalently bonded molecular complexes 29. The complexes in S22 are bonded by hydrogen bonds, dispersion bonds, or a mix of the two. The size of each complex varies between 6 and 30 atoms. The complexes in S22 contain both intermolecular and intramolecular vibrational degrees of freedom. Note that there are always six intermolecular vibrational modes for each complex. By predicting their vibrational frequencies, it will be possible to examine a given XC functional’s ability to describe both covalent and non-covalent interactions in the same system.

The results for the vibrational frequencies of the water dimer are presented in Table 4 with internal coordinate decomposition for each mode. We include frequencies from Xu and Goddard 62 predicted with the GGA functional BLYP 63, 64 and the hybrid functional B3LYP 63, 64, 65, 66. From Dunn et al. 67, we include frequencies predicted with MP2 theory 68 using the aug-cc-pVDZ basis set. The results for the entire S22 dataset are available in the Supporting Information.

BEEF-vdW predicts the highest intramolecular frequencies (modes 7 to 12), with values comparable to the B3LYP frequencies reported by Xu and Goddard. With respect to the experimental intramolecular frequencies reported by Fredin et al.69, BEEF-vdW is the best functional out of the six tested in this work, with an MAE of 3.8 meV. The best results come via Xu and Goddard 62 using B3LYP. A comparison of intermolecular frequencies is not possible because, with the exception of mode 4, the experimental measurements have not been corrected for anharmonic effects, and therefore the reported frequencies are not harmonic. We do note, however, that BEEF-vdW strongly overestimates the intermolecular frequencies in comparison to the other XC functionals. All six experimental intramolecular frequencies are bounded to within one standard deviation of the BEEF-vdW value, with an average deviation of $0.49\sigma$. Ensembles for the intramolecular modes all have COV lower than $0.05$, which is consistent with ensembles for H2O reported in Table 1. Ensembles for intermolecular modes have relatively high COV, comparable to the COV values of the C2H6 torsional mode and the rare-gas dimer modes, with the highest (1.22) coming from the lowest frequency mode. Intermolecular modes of other S22 complexes also tend to have higher ensemble COV values. All six intermolecular ensembles for the water dimer are non-Gaussian, while this is true of only two of the six intramolecular modes. As with the C2H6 torsional mode, some BEEF-vdW ensemble functionals predict unstable imaginary frequencies for mode 1 (458 functionals) and mode 2 (23 functionals). Challenges in describing interaction of small molecule clusters with different XC approximations has been previously noted 70.

3.4 Comparison

Based on the decomposition of the normal modes in terms of stretching, bending, and torsional motion for every molecule and complex considered, we now examine the relationship between the physical motion of the mode and the spread of its BEEF-vdW ensemble, as measured by COV. In Figure 3(a), the bending component is plotted against the stretching component with marker color determined by the magnitude of the COV. Figure 3(b) is similar, with the bending component plotted versus the torsional component. The data plotted in Figures 3(a) and 3(b) indicate that COV tends to (i) increase as the component of bending or torsion increases and (ii) decrease as the stretching component increases. The rare-gas dimers, which have high COV stretch modes, are an exception.

The COV is plotted as a function of the BEEF-vdW ensemble mean frequency in Figure 4. The physical character of the mode is indicated by the marker color, with red, green, and blue representing pure stretch, bend, and torsion. The intermolecular modes for the S22 complexes are marked in black. The COV decreases as the mean frequency increases, with the S22 intermolecular modes having the largest spreads. Modes with frequencies higher than 200 meV are almost entirely stretch modes with low COV. Torsional and bending modes, rare-gas dimer stretching modes, and intermolecular modes tend to have frequencies lower than 200 meV and higher COV.

The COV for all modes, organized by dataset and molecule, are provided in Figure 5. Within each dataset, the molecules and complexes are listed in order of increasing mass. An analogous plot with a logarithmic vertical axis is provided in the Supporting Information. From Figures 4 and 5, a general ordering of mode types based on their COV can be observed. Stretching modes tend to have the lowest COV, followed by bending modes, torsional modes, and intermolecular modes. For every molecule in S22, nearly all intermolecular modes have higher COV than any other mode. The relatively high COV values for bending and torsional modes, in contrast with the stretching modes, indicates disagreement among XC functionals in describing these vibrations. As stretching modes are typically localized while bending and torsional modes are not74, our results suggest that XC functionals at the GGA level have difficulty consistently describing a molecule’s non-local vibrational behavior. The stretching modes of the rare-gas dimers have high COV and are an exception to this rule. These systems are bonded by van der Waals interactions, however, and the high COV values correctly reflect disagreement in XC functionals in describing non-covalent interactions. This inconsistency also causes the high COV values of the intermolecular S22 modes.

Based on Figure 5 there is a separation of modes at a COV value of $0.25$, which is marked by a horizontal dashed line. Modes above this boundary include the high COV stretch modes of the rare-gas dimers, most of the S22 intermolecular modes, and some bending and torsional modes from C2H6 and several S22 complexes. Based on this separation, we propose a criterion of a COV equal to $0.25$ for disagreement among XC functionals in predicting the frequency. For ensembles with COV larger than $0.25$, extra care should be taken in choosing an XC functional for the frequency prediction.

4 Summary

We presented a computationally efficient method to quantify the uncertainty due to the choice of the XC functional at the GGA level in DFT predictions of harmonic vibrational frequencies. To test the robustness of this method, we considered three sets of molecules of varying size and complexity. The first set consisted of small benchmark molecules (Tables 1 and 2), none of which were larger than 8 atoms. We found that the BEEF-vdW ensemble bounds most experimental frequencies and the frequency predictions of six other XC functionals to within one standard deviation (e.g., Figure 2). We then applied the method to ten rare-gas dimers (Table 3) and the S22 dataset (Table 4) of molecular complexes, which offered a variety of systems that differed in their bonding environments (i.e., covalent, hydrogen, van der Waals) and physical nature of the their vibrational modes (i.e., stretch, bend, and torsion, and intermolecular). Our results show that frequency predictions for modes with delocalized motion, such as torsion and bending, and modes involving non-covalent bonds are more sensitive to choice of XC functional in comparison to predictions for localized and covalent stretch modes (Figures 3, 4, and 5). Our proposed method can therefore be used to link DFT uncertainty to the physical behavior of the system.

{acknowledgement}

H. L. P. acknowledges support from a Presidential Fellowship at Carnegie Mellon University and helpful feedback from Dilip Krishnamurthy. V. V. acknowledges support from the National Science Foundation under award CBET-1554273.