Can Density Functional Theory Be Trusted for High-Order Electric

May 13, 2019 - ... value for the excess nuclear relaxation (NR) polarizability was found ...... E. Exchange-correlation potential with correct asympto...
0 downloads 0 Views 865KB Size
Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. 2019, 15, 3570−3579

Can Density Functional Theory Be Trusted for High-Order Electric Properties? The Case of Hydrogen-Bonded Complexes Robert Zaleśny,*,† Miroslav Medved’,‡ Sebastian P. Sitkiewicz,¶,§,∥ Eduard Matito,*,§,▽ and Josep M. Luis*,∥

Downloaded via UNIV OF SOUTHERN INDIANA on July 18, 2019 at 10:17:35 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



Department of Physical and Quantum Chemistry, Faculty of Chemistry, Wrocław University of Science and Technology, Wyb. Wyspiańskiego 27, PL−50370 Wrocław, Poland ‡ Department of Chemistry, Faculty of Natural Sciences, Matej Bel University, Tajovského 40, 974 01 Banská Bystrica, Slovak Republic ¶ Kimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/EHU), 20080 Donostia, Euskadi, Spain § Donostia International Physics Center (DIPC), Manuel Lardizabal Ibilbidea 4, 20018 Donostia, Euskadi, Spain ∥ Institute of Computational Chemistry and Catalysis and Department of Chemistry, University of Girona, Campus de Montilivi, 17003 Girona, Catalonia, Spain ▽ Ikerbasque Foundation for Science, 48011 Bilbao, Euskadi, Spain S Supporting Information *

ABSTRACT: This work reports on an extensive assessment of the performance of a wide palette of density functional approximations in predicting the (high-order) electric properties of hydrogen-bonded complexes. To this end, we compute the electronic and vibrational contributions to the electric polarizability and the first and second hyperpolarizabilities, using the CCSD(T)/aug-cc-pVTZ level of theory as reference. For all the studied properties, the average absolute errors below 20% can only be obtained using the CAMB3LYP functional, while LC-BLYP and MN15 are shown to be only slightly less accurate (average absolute errors not exceeding 30%). Among Minnesota density functionals, i.e., M06, M06-2X, and MN15, we only recommend the latter one, which quite accurately predicts the electronic and vibrational (hyper)polarizabilities. We also analyze the optimal tuning of the range-separation parameter μ for the LC-BLYP functional, finding that this approach does not bring any systematic improvement in the predictions of electronic and vibrational (hyper)polarizabilities and the accuracy of computed properties is largely system-dependent. Finally, we report huge errors in predicting the vibrational second hyperpolarizability by ωB97X, M06, and M06-2X functionals. Based on the explicit evaluation of anharmonic terms contributing to the second hyperpolarizability, this failure is traced down to a poor determination of third- and fourth-order energy derivatives with respect to normal modes. These results reveal serious flaws of some density functional approximations and suggest caution in selecting the appropriate functional to calculate not only electronic and vibrational (hyper)polarizabilities but also other molecular properties that contain vibrational anharmonic contributions.



nature,11,12 which gives even stricter requirements on the accuracy of applied computational methods.4,13 Density Functional Theory (DFT) has proven to be a powerful tool for solving various quantum mechanical problems in a costeffective way, but since the seminal works of Champagne and co-workers,14−16 DFT performance in the area of molecular nonlinear optics has been under active assessment. Density functional approximations (DFAs) with a large amount of Hartree−Fock exchange have been recommended to avoid the large self-interaction or delocalization error17 and the wrong exchange potential decay of generalized gradient approxima-

INTRODUCTION The bottom-up design of new molecular materials possessing large nonlinear optical (NLO) responses in condensed phases (liquids, crystals, thin films, etc.) remains a challenging task for computational chemistry due to the necessity of accurate descriptions of both electronic and vibrational (hyper)polarizabilities of individual molecules as well as an appropriate inclusion of the effects of intermolecular interactions in the bulk.1−5 Although the former usually requires high-level quantum-chemical methods and large basis sets,6−10 the molecular size often becomes a limiting factor for the application of such methods. On the other hand, the interactions with surrounding molecules are by definition long-range and have dynamic and cooperative (many-body) © 2019 American Chemical Society

Received: February 12, 2019 Published: May 13, 2019 3570

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation tions (GGAs).18 On the contrary, DFAs that do not include HF exchange often provide too delocalized structures, leading to the overestimation of conjugation and aromaticity.19−22 Long-range corrected DFAs (LC-DFAs),23−25 including a variable amount of HF exchange at different interelectronic ranges, were found to be successful in reproducing the results of high-level correlation methods and/or experimental values of both electronic and vibrational (hyper)polarizabilities of isolated molecules including charge-transfer, highly conjugated, and zwitterionic systems.17,26−31 However, as shown by some of the present authors, in some instances LC-DFAs do not bring improvement in the description of vibrational hyperpolarizabilities of π-conjugated systems over conventional functionals.32 LC-DFAs employ an attenuating function to determine the mixture of pure and HF exchange at each interelectronic separation. This attenuating function is often an error function that depends on three tunable parameters. The LC-DFAs that adjust the latter parameters to satisfy the ionization potential theorem (i.e., minimizing the sum of the highest-occupied Kohn−Sham orbital energy and the ionization potential) are known as optimal or tuned LC-DFAs (OLC-DFAs).33 OLC-DFAs provide an improved description of charge-transfer effects, UV−vis spectrum, and other properties involving frontier orbitals. However, the appropriateness of using OLC-DFAs to evaluate molecular NLO properties is still unclear.34−37 In the case of resonant NLO properties of molecules, like two-photon absorption, the treatment of vibrational contributions using DFT may lead to catastrophic predictions of vibrationally resolved spectra.38,39 Let us recall that the nonresonant vibrational contributions to hyperpolarizability were in many cases found to be of the same magnitude or even larger than their electronic counterparts.40−42 As a natural step toward more complex systems, the electronic NLO properties of various types of weakly bound systems including hydrogenbonded and van der Waals complexes were intensively studied by high-level post-Hartree−Fock ab initio as well as DFT methods.11,12,43−58 It was found that the properties of hydrogen bonded systems can be satisfactorily described by LC-DFAs and also by some global hybrids (e.g., M06-2X).49 In the case of van der Waals systems, these methods also performed well provided that the geometries used in the property calculations were optimized taking into account the dispersion interaction.56,58 The good performance of LC-DFAs could be traced back not only to an improved description of the NLO properties of the individual subsystems but also to a better description of the interaction-induced (excess) component of NLO properties involving the effects of intermolecular interactions in the complex.59 On the other hand, studies of the vibrational NLO properties of weakly bound complexes are much scarcer,51,59 although the presence of low-frequency vibrational modes in such systems suggests that the nuclear motions can have higher impact on the properties of weakly bound complexes than on the properties of the constituent isolated molecules. From a methodological point of view, it should be highlighted that the accurate description of vibrational NLO properties demands stricter requirements on the applied approach than the calculation of their electronic counterpart, since for the former the equilibrium geometry, vibrational frequencies and electrical and mechanical anharmonic corrections should be calculated at the same level of theory (see below). As already mentioned, the magnitude of vibrational contributions to hyperpolariz-

ability is often comparable to their electronic counterparts. This, in general, might hold for static electric dipole hyperpolarizabilities (studied in the present work).60 The role of the vibrational hyperpolarizabilities is also important for electro-optical phenomena like the Pockels and DC Kerr effects and for optical processes like the intensity dependent refractive index.61 On the contrary, for second and third harmonic generation hyperpolarizabilities, when the optical frequencies are much higher than the vibrational frequencies, the vibrational hyperpolarizabilities become far smaller than their electronic counterparts.62 In our pilot study59 we performed an interaction energy decomposition analysis of the physical origins of excess vibrational NLO contributions for the HCN dimer. In that work, we also employed three DFAs, namely BLYP and LC-BLYP (and its variant including density-dependent dispersion correction,63 LC-BLYP-dDsC) showing that LC-BLYP provided the most promising results compared to CCSD(T), whereas LC-BLYP-dDsC giving a very reasonable value for the excess nuclear relaxation (NR) polarizability was found numerically unstable for higher-order properties. The satisfactory performance of LC-DFAs observed for HCN···HCN encouraged us to undertake the present study, which aims at benchmarking a wide palette of DFAs in predicting anharmonic vibrational (and electronic) (hyper)polarizabilities of weakly bound molecular complexes. To this end, we have selected 9 hydrogen-bonded systems and used the results of highly accurate CCSD(T)/aug-cc-pVTZ calculations (employing frozen-core approximation),64 as reference for the assessment of diverse set of DFAs approaches. The present contribution focuses on molecular complexes, and it is thus a substantial extension to our previous paper devoted to electronic and vibrational (hyper)polarizabilities of organic π-conjugated molecules.32



THEORY AND COMPUTATIONAL DETAILS Under the Born−Oppenheimer approximation, one may define the total property P as the sum of electronic (Pel), nuclear relaxation (Pnr), and curvature (Pcurv) contributions:65 P = P el + P nr + P curv

(1)

Bishop, Hasan, and Kirtman (BHK) showed that one may employ a computationally efficient finite-field nuclear relaxation (FF-NR) approach to evaluate the nuclear-relaxation contribution to the vibrational (hyper)polarizabilities. The Pnr contains the lowest-order terms and is usually the largest contribution to the vibrational (hyper)polarizabitiles. In the present study, we focus solely on static properties, including polarizability (α), first hyperpolarizability (β), and second (γ) hyperpolarizability. As far as the nuclear relaxation contribution is concerned, following BHK we define66 μi (F , R F) − μi (0, R 0) =

∑ aij1Fj + j

1 + 6

1 ∑ gijkl FjFkFl + ...

1 2

∑ bijk1 FjFk jk

(2)

jkl

where μ(F,RF) is an electronic property obtained at the fieldrelaxed geometry, and μ(0,R0) is the same property at field-free conditions. The expansion coefficients yield the static properties: aij1 = αijel + αijnr = αijel + [μ2 ](0,0) ij 3571

(3)

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation

of HOMO and the first ionization potential of a system, the resulting functional being labeled as OLC-BLYP.91 We also implemented the constrain of satisfying the IP theorem on N+1 species, but we did not add these results because the EA was negative in all cases. The use of the IP theorem in these cases is, at the very least, questionable. However, in all cases the EA was so close to zero that it did not affect the parameter optimization. Therefore, even after considering the IP theorem on N+1 species for the optimization of the range-separation parameter, we obtained little improvement of the NLOPs with respect to the case where we only use the IP theorem on the N species. All calculations were performed with the aug-cc-pVTZ basis set.92 As demonstrated in our previous study, this basis set in conjunction with CCSD(T) allows for obtaining 1% accuracy in predictions of electronic and nuclear relaxation (hyper)polarizabilities of small molecules.93

1 bijk = βijkel + βijknr = βijkel + [μα ](0,0) + [μ3](1,0) + [μ3](0,1) ijk ijk ijk

(4) 1 el nr gijkl = γijkl + γijkl el 2 (1,0) 2 (0,1) (0,0) = γijkl + [α 2](0,0) ijkl + [μβ ]ijkl + [μ α ]ijkl + [μ α ]ijkl 4 (0,2) 4 (1,1) + [μ4 ](2,0) ijkl + [μ ]ijkl + [μ ]ijkl

(5)

Each square bracket involves the products of normal coordinate derivatives of the electronic electrical properties indicated, as well as harmonic vibrational frequencies and anharmonic force constants.67 Each term is identified by a pair of superscripts (n,m) denoting the order in electrical and mechanical anharmonicities, respectively. For the sake of brevity, in what follows the tensor indices will be omitted because we will focus on diagonal (longitudinal) properties. In applying the FF-NR formulas, one should take into account that the geometry relaxation must not result from rotations of the molecule through alignment of the permanent and(or) induced dipole moment in the field direction (indeed this may be the easiest way for the system to lower its energy). For that reason, the field-dependent optimization should be performed while strictly maintaining the Eckart conditions.68 Such fielddependent geometry optimizations can be carried out with the aid of the procedure developed by Luis et al.69 Moreover, in order to obtain control over the numerical differentiation error in the FF-NR approach, we employed the Romberg−Rutishauser scheme.70,71 The correction of order p to the derivative P is determined iteratively according to the formula P p,k =

4 p·P p − 1, k − P p − 1, k − 1 4p − 1



RESULTS AND DISCUSSION As highlighted in the Introduction, in the present study the focus is put on the electric properties of molecular complexes. To this end, we will present the results obtained for longitudinal electronic and vibrational (hyper)polarizabilities of the following set of hydrogen-bonded systems: HCN···HCl, HCN···HCN, HCN···HF, HCN···HNC, HNC···HCN, N2··· HF, OC···HF, FCCH···NCF, and HCCH···NCF. The reference values were determined at the CCSD(T)/aug-ccpVTZ level of theory. For the sake of computing various anharmonic terms, we also employed the MP2 method, and its performance will also be discussed in connection with the vibrational properties. Unless otherwise indicated, the properties (electronic and vibrational) and equilibrium geometries were determined at the very same level of theory. The summary of electric property calculations is presented in Figure 1. Note that, for the sake of clarity, errors which amount

(6)

where k is related to the strength of the field used in the FFNR calculations (see below). The value of P0,k is determined from the finite−difference expression for the derivative of the energy (E) or dipole moment with respect to the electric field (F). The smallest electric field amplitude Fi used was equal to 0.0002 au. In all calculations employing the Romberg− Rutishauser scheme, we used the following sequence of fields ±Fi, ±2Fi, ..., ±2kFi, where k = 7. All FF-NR results were rechecked using field-induced coordinates61,72 (FIC) and Bishop−Kirtman perturbation theory73,74 (BKPT) methods. In a few problematic cases, for which FF-NR calculations were burdened by large numerical errors, we report BKPT values for total nuclear-relaxation properties. Moreover, the latter method was used to break down nuclear-relaxation hyperpolarizabilities into anharmonic contributions. FF-NR, FIC, and BKPT calculations were performed using custom computer routines based on the energies (FF-NR method) and property derivatives (FIC and BKPT methods) computed using the GAUSSIAN package (ver. 09 and 16).75 Tight convergence criteria for energy and density were imposed (energy convergence criterion was set to 10−12 hartree), and we used a pruned superfine grid ((175,974) for first-row atoms and (250,974) for atoms in the second and third rows). The DFT calculations were performed employing 10 functionals, i.e., BLYP,76,77 B3LYP,78 CAM-B3LYP,79 LCBLYP (μ = 0.47),80 ωB97X,81 M06,82 M06-2X,82 MN15,83 PBE0,84,85 and HSE06.86−90 The range-separation parameter of LC-BLYP, μ, was tuned to minimize the sum of the energy

Figure 1. Average absolute error in predicting electric properties for a set of 9 molecular complexes. Errors which amount to 100% (and larger) are marked with the same color.

to 100% and larger are marked with the same color (red). It turns out that, except for the electronic polarizability αel, the average absolute errors for all properties show large variations. By and large, the errors of many DFAs for the electronic as well as vibrational counterparts exceed 25%. To shed light on these striking unsystematic trends we will start discussing the electronic contributions, which will be followed by the indepth analysis of the vibrational counterparts. All the underlying numerical data are included in the Supporting Information. 3572

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation

geometries obtained using the CCSD(T) method. The comparison of the results presented in Table S9 with their counterparts presented in Figure 2 allows for drawing the conclusion that the choice of reference geometry has an insignificant impact on the computed electronic properties. The average and maximum absolute relative errors in predicting αnr, βnr, and γnr are shown in Figures 3a, 3b, and 3c,

The average and maximum absolute relative errors in predicting αel, βel, and γel are shown in Figures 2a, 2b, and 2c,

Figure 2. Average and maximum absolute relative errors in predicting electronic polarizability (a), first hyperpolarizability (b), and second hyperpolarizability (c) for a set of 9 molecular complexes with respect to CCSD(T)/aug-cc-pVTZ reference values. Figure 3. Average and maximum absolute relative errors in predicting nuclear-relaxation polarizability (a), first hyperpolarizability (b), and second hyperpolarizability (c) for a set of 9 molecular complexes with respect to CCSD(T)/aug-cc-pVTZ reference values.

respectively. The average absolute errors in αel span from 0.80% (CAM-B3LYP) up to 7.99% (BLYP), and the maximum errors for all DFAs but BLYP are below 5%. In the case of βel, the average errors are much larger, and they are in the range of 13.40% (M06-2X)−118.06% (BLYP), the maximum error for DFAs being at least 45%. The BLYP functional is also the worst in predicting γel (128.19%). The property in question is the most accurately predicted by the MN15 functional, with an average absolute error of only 3.96%. The results presented in Figure 2, taken together, indicate that the pure generalized gradient approximation (GGA) BLYP functional is systematically the worst in predicting electronic (hyper)polarizabilities, and this is in line with other works.14,15,94 On the other hand, absolute average errors smaller than 20% can only be achieved by using the long-range corrected CAM-B3LYP and metaGGA M06-2X functionals. Note that the MP2 method on average performs slightly worse than these two functionals in predicting electronic (hyper)polarizabilities. As explained in the preceding paragraph, all the properties (electronic and vibrational) were obtained at the equilibrium geometries corresponding to the indicated level of theory. Table S9 in the Supporting Information contains the comparison of electronic properties computed with the set of DFAs at the equilibrium

respectively. In the case of αnr, the average errors are in the range 3.60% (ωB97X)−31.36% (BLYP). It is somewhat unexpected that the best functional for computing αnr shows the worst performance in predicting βnr, i.e., the corresponding average absolute error being 85.57%. The smallest average error for this property is found for the CAM-B3LYP functional (only 13.73%). Likewise, CAM-B3LYP and ωB97X deliver the smallest and the largest average errors of γnr, 13.60% and 1767.63%, respectively. The Supporting Information contains the plots showing system-specific properties (see Figure S1). It follows from these plots that only for electronic polarizability all studied DFAs predict the same ordering of molecular complexes as CCSD(T). Nevertheless, at this point it is interesting to remark that CAM-B3LYP also reproduces CCSD(T) ordering of molecular complexes for γel, αnr, and βnr and only presents one error in the ordering for βel and γnr. The huge errors in predicting γnr are striking and deserve a deeper analysis. However, computations of explicit harmonic 3573

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation and anharmonic terms to γnr at the CCSD(T) level are unfeasible, and thus this approach cannot be used to decompose the nuclear relaxation contributions. Opportunely, the accuracy of the MP2 method in predicting the properties in question is much better than the accuracy obtained by using any of the previously discussed DFAs (see Figure 3). The MP2 average relative errors for αnr, βnr, and γnr are 2.37%, 8.72%, and 9.95%, respectively. Note that in the case of γnr the MP2 method yields the smallest maximum and average errors. It can thus be concluded that, given the huge errors in γnr of some DFAs, MP2 can be used as a reference method for the analysis of the sources of the poor performance of some of the DFAs. Figure 4 contains the net harmonic and

Δ̅ anh (DFA) =

n

∑ i=1

(DFA) − γi nr (MP2)] [γi nr ,har ,har γi nr(MP2)

i=1

(DFA) − γi nr (MP2)] [γi nr ,anh ,anh γi nr(MP2)

× 100%

Δ̅ har and Δ̅ anh are shown in Figure 5. The analysis of these figures suggests that the huge errors in predicting γnr by M06,

Figure 5. Average and maximum absolute relative errors in predicting harmonic (top) and anharmonic (bottom) contribution to γnr for a set of 9 molecular complexes with respect to MP2/aug-cc-pVTZ. See eq 7 for the definition of the relative error.

M06-2X, and ωB97X (see Figure 3c) can indeed be associated with anharmonic terms contributing to γnr. In order to pinpoint the sources of these errors to specific anharmonic terms, shown in eq 5, the absolute relative errors for three functionals (M06, M06-2X, and ωB97X) were computed. Figure 6 also includes the absolute relative error in γnr with respect to MP2. Two interesting conclusions can be drawn: i) the largest errors are found for [μ4] terms and ii) the ωB97X functional predicts these terms with largest errors. To shed more light on the poor performance of the ωB97X functional, in Figure 7 the values of three [μ4] terms are shown for each studied molecular

anharmonic contributions of βnr and γnr for all systems computed with the MP2 method. The anharmonicity of the first hyperpolarizability, at best, reaches up to 50% (HCN··· HF) of the total nuclear relaxation contribution to β. In the case of second hyperpolarizability, however, the anharmonic corrections become more important and prevail over the harmonic terms. Except for HCCH···NCF (34.48%) and FCCH···NCF (49.63%) complexes, the anharmonicity contributions exceed 50% of the total γnr value and reach up to 86.99% (HCN···HCl). The relative weight of the anharmonicity in the vibrational hyperpolarizabilities cannot be predicted only by analyzing the effect of the anharmonicity in the vibrational frequencies and intensities (see Table S14 in the Supporting Information).95 Having established that anharmonic contributions to γnr prevail over the harmonic ones, we will now make an attempt to determine which terms of γnr (see eq 5) cause the bad performance of some DFAs in predicting this property. The average relative error in predicting γnr (for the set of n = 9 complexes) can be split into harmonic and anharmonic parts 1 n

n



(8)

Figure 4. Percentage weight of net harmonic and anharmonic contributions to nuclear-relaxation first and second hyperpolarizability computed at the MP2/aug-cc-pVTZ level of theory.

Δ̅ har (DFA) =

1 n

× 100%

Figure 6. Average absolute relative error for three functionals in predicting anharmonic contribution to γnr for a set of 9 molecular complexes with respect to MP2/aug-cc-pVTZ.

(7) 3574

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation

these findings are in line with the results presented herein and the work by Biczysko et al.97 The satisfactory performance of the LC-BLYP functional (for each of the studied properties the relative errors do not exceed 30%) encouraged us to analyze the effect of tuning the range-separation parameter μ. Table 1 contains the comparison of percentage relative errors for LC-BLYP and OLC-BLYP functionals for three selected molecular complexes. Whereas for the HCN···HF complex the tuned μ value is rather large and amounts to 0.545, it remains very close to a default value (0.47) for HCN···HCl (0.450) and HCN···HCN (0.455). As demonstrated in Table 1, the optimal tuning of the μ parameter brings either little improvement (HCN···HCl) or worsens (HCN···HCN, HCN···HF) the predictions of the electronic polarizability, in line with previous findings.35 For this property, the overall changes in relative errors with respect to CCSD(T) values due to the tuning of the μ parameter do not exceed 2.5%. In the case of βel, OLC-DFT leads to the deterioration of predictions of this property for all three molecular complexes, the corresponding relative errors with respect to the CCSD(T) reference increasing by 3−15% with respect to LC-BLYP results. The tuning of the μ parameter, by and large, does not improve significantly the predictions of γel, i.e., for HCN···HCl and HCN···HCN there is only an insignificant decrease in the relative error while for HCN··· HF we find a 10% increase in the relative error with respect to the CCSD(T) reference. In the case of nuclear-relaxation contributions, the pattern of changes upon using systemspecific μ values does not follow the trends observed for the electronic counterparts. The changes in relative errors for αnr, on passing from μ = 0.47 to the tuned value, are not large (2− 3%) and go in both directions, either to worsening (HCN··· HCl, HCN···HCN) or to improving (HCN···HF) the predictions of this property. Note that this trend is closely paralleled in the case of βnr and γnr, although the relative errors are slightly larger for HCN···HF. To sum up, the optimal tuning of the μ parameter does not bring any systematic improvement in the predictions of electronic and vibrational (hyper)polarizabilities, and, as demonstrated by the results assembled in Table 1, the success of this approach depends heavily on the property and the system. Table S10 in the Supporting Information contains also the results of calculations using the nonlocal van der Waals density LC-VV10 functional, which is recommended for weakly bound complexes. However, this functional does not bring any systematic improvement over LC-BLYP in predicting electronic and nuclear relaxation (hyper)polarizabilities. These results are in line with previously published data for the HCN···HCN complex using the

Figure 7. Selected anharmonic contributions to nuclear-relaxation second hyperpolarizability (given in au). The meaning of the label (*) is explained in the text.

complex. For each complex the sign of [μ4](1,1) is correctly predicted, whereas [μ4](0,2) and [μ4](2,0) are predicted for several complexes with a wrong sign. The overall performance of the ωB97X functional is thus largely dependent on the error compensation. This also holds for the M06 and M06-2X functionals (see Figures S2 and S3 in the Supporting Information). Figure S4 in the Supporting Information shows that the CAM-B3LYP functional does not suffer from similar error cancellations, i.e., for all studied complexes the signs of all the computed [μ4] terms are correctly predicted. Figure 7 shows that the most problematic complex is HCN··· HCl. In particular, the [μ4](0,2) term is largely overestimated by the ωB97X functional. Within the Bishop−Kirtman perturba∂μ ∂μ ∂μ ∂μ tion theory formalism, this term encompasses a40 ∂Q ∂Q ∂Q ∂Q ∂μ ∂μ ∂μ ∂μ

and a30a30 ∂Q ∂Q ∂Q ∂Q contributions, where a30 and a40 denote third- and fourth-order energy derivatives with respect to normal modes. In an attempt to pinpoint the origins of large error in the [μ4](0,2) term for the HCN···HCl complex, we have performed calculations combining a30 and a40 using MP2 with dipole moment derivatives calculated using the ωB97X functional. As it is shown in Figure 7, the value of the [μ4](0,2) term is now only slightly underestimated, linking the huge ωB97X errors with the calculation of the third- and fourth-order energy derivatives. It has been recently shown by Ravichandran and Banik that some of the functionals studied herein (M06 and M06-2X) exhibit significant errors in predicting anharmonic vibrational frequencies.96 This poor performance was linked to anharmonic force constants, and

Table 1. Percentage Relative Errors in Longitudinal Electronic and Nuclear-Relaxation (Hyper)polarizabilities of Molecular Complexes with Respect to CCSD(T)/aug-cc-pVTZ Reference Values αel

βel

OLC-BLYP (μ = 0.4502) LC-BLYP (μ = 0.4700)

−0.17 −0.83

20.25 16.96

OLC-BLYP (μ = 0.4553) LC-BLYP (μ = 0.4700)

−1.82 −2.22

35.71 32.50

OLC-BLYP (μ = 0.5449) LC-BLYP (μ = 0.4700)

−3.24 −1.04

−64.86 −51.35

γel HCN···HCl −14.47 −17.11 HCN···HCN −18.46 −20.00 HCN···HF −26.67 −16.67 3575

αnr

βnr

γnr

21.93 19.82

26.10 23.21

43.98 42.57

14.01 12.95

14.51 12.86

15.93 12.39

12.03 15.41

20.96 25.62

12.00 20.00

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation dispersion-corrected functional, LC-BLYP-dDsC. These preliminary results indicate that electronic and vibrational hyperpolarizabilities are not much improved using dispersion corrections, at least for the systems considered in the present work.





AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. *E-mail: [email protected].

SUMMARY AND CONCLUSIONS In the present study we performed an extensive assessment of the performance of DFAs in predicting molecular (hyper)polarizabilities of 9 hydrogen-bonded complexes using the CCSD(T)/aug-cc-pVTZ level of theory as reference. The studied properties included electronic and vibrational contributions to the polarizability and the first and second hyperpolarizabilities. The vibrational counterpart was approximated by nuclear-relaxation terms which contain the low-order anharmonic corrections, and it usually represents the dominant contribution to the vibrational (hyper)polarizabilities. The results presented in this work can be summarized as follows. Considering all the studied properties, i.e., the electronic and nuclear-relaxation (hyper)polarizabilities, average absolute errors below 20% can only be obtained using the CAMB3LYP functional. LC-BLYP and MN15 (errors do not exceed 30%) and PBE0 and HSE06 (errors reach up to 35%) are slightly less accurate functionals. Among Minnesota density functionals, we recommend MN15 which, in contrast to its forerunners (M06, M06-2X), quite accurately predicts the electronic and vibrational (hyper)polarizabilities. The dramatic failure of some DFAs (e.g., ωB97X, M06, and M06-2X) in predicting the vibrational second hyperpolarizability was traced down to a poor determination of third- and fourth-order energy derivatives with respect to normal modes. This observation might have serious implications for the accuracy of simulations of, e.g., IR spectra including anharmonicity. In this work we also attempted the optimal tuning of the rangeseparation parameter μ for the LC-BLYP functional, but it does not bring any systematic improvement in predictions of electronic and vibrational (hyper)polarizabilities. On the contrary, the accuracy of computed properties heavily depends on the system and the property, and, in many cases, we have observed a deterioration of the results with respect to default value (μ = 0.47). Hence, at least for nonresonant (hyper)polarizabilities of hydrogen-bonded complexes, we do not recommend this approach. It should be highlighted that the conclusions presented herein were drawn based on the results obtained for a set of small-sized hydrogen-bonded molecular complexes. The validity of these observations for larger molecular systems which could be problematic for commonly used DFAs, e.g., organic π-conjugated charge-transfer complexes, is to be determined. This work is a first step toward thorough assessment of density functional theory in predicting high-order electric properties including anharmonic vibrational corrections.



ability computed using the M06, M06-2X, and CAMB3LYP functionals (PDF)

ORCID

Robert Zaleśny: 0000-0001-8998-3725 Eduard Matito: 0000-0001-6895-4562 Josep M. Luis: 0000-0002-2880-8680 Funding

R.Z. acknowledges financial support by the Polish National Science Centre (Grant No. 2015/19/B/ST4/01881). M.M. acknowledges the Slovak Research and Development Agency (project no. APVV-15-0105). J.M.L. and E.M. are grateful for financial support from the Spanish MINECO CTQ201452525-P (E.M. and J.M.L.), PGC2018-098212-BC21 (E.M.), PGC2018-098212-B-C22 (JML) and EUIN2017-88605 (E.M.), and to the Catalan DIUE 2014SGR931 (J.M.L.). S.P.S. acknowledges the Basque Government for funding through a predoctoral fellowship (PRE 2018 2 0200). The authors gratefully acknowledge Wroclaw Centre for Networking and Supercomputing for the generous allotment of computer time. Notes

The authors declare no competing financial interest.



REFERENCES

(1) Di Bella, S.; Ratner, M. A.; Marks, T. J. Design of chromophoric molecular assemblies with large second-order optical nonlinearities. A theoretical analysis of the role of intermolecular interactions. J. Am. Chem. Soc. 1992, 114, 5842−5849. (2) Facchetti, A.; Annoni, E.; Beverina, L.; Morone, M.; Zhu, P.; Marks, T. J.; Pagani, G. A. Very large electro-optic responses in Hbonded heteroaromatic films grown by physical vapour deposition. Nat. Mater. 2004, 3, 910−917. (3) Datta, A.; Pati, S. K. Dipolar interactions and hydrogen bonding in supramolecular aggregates: understanding cooperative phenomena for 1st hyperpolarizability. Chem. Soc. Rev. 2006, 35, 1305−1323. (4) Seidler, T.; Stadnicka, K.; Champagne, B. Evaluation of the linear and second-order nlo properties of molecular crystals within the local field theory: Electron correlation effects, choice of XC functional, ZPVA contributions, and impact of the geometry in the case of 2-methyl-4-nitroaniline. J. Chem. Theory Comput. 2014, 10, 2114−2124. (5) Medved’, M.; Budzák, Š .; Bartkowiak, W.; Reis, H. Solvent effects on molecular electric properties. In Handbook of Computational Chemistry; Springer International Publishing: 2017; Vol. 17, pp 741− 795, DOI: 10.1007/978-3-319-27282-5_44. (6) Shelton, D. P.; Rice, J. E. Measurements and calculations of the hyperpolarizabilities of atoms and small molecules in the gas phase. Chem. Rev. 1994, 94, 3−29. (7) López Cacheiro, J.; Fernandez, B.; Marchesan, D.; Coriani, S.; Hättig, C.; Rizzo, A. Coupled cluster calculations of the ground state potential and interaction induced electric properties of the mixed dimers of helium, neon and argon. Mol. Phys. 2004, 102, 101−110. (8) Pecul, M.; Pawłowski, F.; Jørgensen, P.; K ohn, A.; Hättig, C. High-order correlation effects on dynamic hyperpolarizabilities and their geometric derivatives: A comparison with density functional results. J. Chem. Phys. 2006, 124, 114101. (9) Hammond, J. R.; Kowalski, K. Parallel computation of coupledcluster hyperpolarizabilities. J. Chem. Phys. 2009, 130, 194108.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.9b00139. Numerical data used to perform statistical analysis discussed in manuscript, dependence of longitudinal electronic (hyper)polarizabilities on geometry of selected molecular complexes, selected anharmonic contributions to nuclear-relaxation second hyperpolariz3576

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation (10) de Wergifosse, M.; Champagne, B. Electron correlation effects on the first hyperpolarizability of push−pull pi-conjugated systems. J. Chem. Phys. 2011, 134, 074113. (11) Skwara, B.; Bartkowiak, W.; Zawada, A.; Góra, R. W.; Leszczynski, J. On the cooperativity of the interaction-induced (hyper)polarizabilities of the selected hydrogen-bonded trimers. Chem. Phys. Lett. 2007, 436, 116−123. (12) Góra, R. W.; Zaleśny, R.; Zawada, A.; Bartkowiak, W.; Skwara, B.; Papadopoulos, M. G.; Silva, D. L. Large changes of static electric properties induced by hydrogen bonding: An ab initio study of linear HCN oligomers. J. Phys. Chem. A 2011, 115, 4691−4700. (13) Kanis, D. R.; Ratner, M. A.; Marks, T. J. Design and construction of molecular assemblies with large second-order optical nonlinearities. Quantum chemical aspects. Chem. Rev. 1994, 94, 195− 242. (14) Champagne, B.; Perpète, E. A.; van Gisbergen, S. J. A.; Baerends, E. J.; Snijders, J. G.; Soubra-Ghaoui, C.; Robins, K. A.; Kirtman, B. Assessment of conventional density functional schemes for computing the polarizabilities and hyperpolarizabilities of conjugated oligomers: An ab initio investigation of polyacetylene chains. J. Chem. Phys. 1998, 109, 10489−10498. (15) Champagne, B.; Perpète, E. A.; Jacquemin, D.; van Gisbergen, S. J. A.; Baerends, E. J.; Soubra-Ghaoui, C.; Robins, K. A.; Kirtman, B. Assessment of conventional density functional schemes for computing the dipole moment and (hyper)polarizabilities of push−pull π− conjugated systems. J. Phys. Chem. A 2000, 104, 4755−4763. (16) Kirtman, B.; Bonness, S.; Ramirez-Solis, A.; Champagne, B.; Matsumoto, H.; Sekino, H. Calculation of electric dipole (hyper)polarizabilities by long−range−correction scheme in density functional theory: A systematic assessment for polydiacetylene and polybutatriene oligomers. J. Chem. Phys. 2008, 128, 114108. (17) Garza, A. J.; Wazzan, N. A.; Asiri, A. M.; Scuseria, G. E. Can short- and middle-range hybrids describe the hyperpolarizabilities of long-range charge-transfer compounds? J. Phys. Chem. A 2014, 118, 11787−11796. (18) Van Leeuwen, R.; Baerends, E. Exchange-correlation potential with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1994, 49, 2421. (19) Wannere, C. S.; Sattelmeyer, K. W.; Schaefer, H. F., III; Schleyer, P. v. R. Aromaticity: The alternating C−C bond length structures of [14]-,[18]-, and [22] annulene. Angew. Chem., Int. Ed. 2004, 43, 4200−4206. (20) Torrent-Sucarrat, M.; Navarro, S.; Cossío, F. P.; Anglada, J. M.; Luis, J. M. Relevance of the DFT method to study expanded porphyrins with different topologies. J. Comput. Chem. 2017, 38, 2819−2828. (21) Szczepanik, D. W.; Solà, M.; Andrzejak, M.; Pawełek, B.; Dominikowska, J.; Kukułka, M.; Dyduch, K.; Krygowski, T. M.; Szatylowicz, H. The role of the long-range exchange corrections in the description of electron delocalization in aromatic species. J. Comput. Chem. 2017, 38, 1640−1654. (22) Casademont-Reig, I.; Woller, T.; Contreras-García, J.; Alonso, M.; Torrent-Sucarrat, M.; Matito, E. New electron delocalization tools to describe the aromaticity in porphyrinoids. Phys. Chem. Chem. Phys. 2018, 20, 2787−2796. (23) Savin, A. A combined density functional and configuration interaction method. Int. J. Quantum Chem. 1988, 34, 59−69. (24) Savin, A. Correlation Contributions from Density Functionals. In Density functional methods in chemistry; Springer: 1991; pp 213− 230, DOI: 10.1007/978-1-4612-3136-3_14. (25) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A long-range correction scheme for generalized-gradient-approximation exchange functionals. J. Chem. Phys. 2001, 115, 3540−3544. (26) Jacquemin, D.; Perpète, E. A.; Medved’, M.; Scalmani, G.; Frisch, M. J.; Kobayashi, R.; Adamo, C. First hyperpolarizability of polymethineimine with long-range corrected functionals. J. Chem. Phys. 2007, 126, 191108. (27) Limacher, P. A.; Mikkelsen, K. V.; Lüthi, H. P. On the accurate calculation of polarizabilities and second hyperpolarizabilities of

polyacetylene oligomer chains using the CAM-B3LYP density functional. J. Chem. Phys. 2009, 130, 194114. (28) Zaleśny, R.; Bulik, I. W.; Bartkowiak, W.; Luis, J. M.; Avramopoulos, A.; Papadopoulos, M. G.; Krawczyk, P. Electronic and vibrational contributions to first hyperpolarizability of donoracceptor-substituted azobenzene. J. Chem. Phys. 2010, 133, 244308. (29) Medved’, M.; Budzák, Š .; Pluta, T. Static NLO responses of fluorinated polyacetylene chains evaluated with long-range corrected density functionals. Chem. Phys. Lett. 2011, 515, 78−84. (30) Sun, H.; Autschbach, J. Influence of the delocalization error and applicability of optimal functional tuning in density functional calculations of nonlinear optical properties of organic donor−acceptor chromophores. ChemPhysChem 2013, 14, 2450−2461. (31) Garza, A. J.; Osman, O. I.; Wazzan, N. A.; Khan, S. B.; Asiri, A. M.; Scuseria, G. E. Prediction of the linear and nonlinear optical properties of tetrahydronaphthalone derivatives via long-range corrected hybrid functionals. Mol. Phys. 2014, 112, 3165−3172. (32) Bulik, I. W.; Zaleśny, R.; Bartkowiak, W.; Luis, J. M.; Kirtman, B.; Scuseria, G. E.; Avramopoulos, A.; Reis, H.; Papadopoulos, M. G. Performance of density functional theory in computing nonresonant vibrational (hyper)polarizabilities. J. Comput. Chem. 2013, 34, 1775− 1784. (33) Baer, R.; Livshits, E.; Salzner, U. Tuned range-separated hybrids in density functional theory. Annu. Rev. Phys. Chem. 2010, 61, 85−109. (34) Autschbach, J.; Srebro, M. Delocalization error and ‘‘functional tuning’’ in Kohn−Sham calculations of molecular properties. Acc. Chem. Res. 2014, 47, 2592−2602. (35) Garrett, K.; Sosa Vazquez, X.; Egri, S. B.; Wilmer, J.; Johnson, L. E.; Robinson, B. H.; Isborn, C. M. Optimum exchange for calculation of excitation energies and hyperpolarizabilities of organic electro-optic chromophores. J. Chem. Theory Comput. 2014, 10, 3821−3831. (36) Oviedo, M. B.; Ilawe, N. V.; Wong, B. M. Polarizabilities of πconjugated chains revisited: Improved results from broken-symmetry range-separated DFT and new CCSD(T) benchmarks. J. Chem. Theory Comput. 2016, 12, 3593−3602. (37) Karolewski, A.; Kronik, L.; Kümmel, S. Using optimally tuned range separated hybrid functionals in ground-state calculations: consequences and caveats. J. Chem. Phys. 2013, 138, 204115. ́ (38) Zaleśny, R.; Tian, G.; Hättig, C.; Bartkowiak, W.; Ågren, H. Toward assessment of density functionals for vibronic coupling in two-photon absorption: A case study of 4-nitroaniline. J. Comput. Chem. 2015, 36, 1124−1131. (39) Bednarska, J.; Zaleśny, R.; Tian, G.; Murugan, N. A.; Ågren, H.; Bartkowiak, W. Nonempirical simulations of inhomogeneous broadening of electronic transitions in solution: Predicting band shapes in one- and two-photon absorption spectra of chalcones. Molecules 2017, 22, 1643. (40) Kirtman, B.; Champagne, B.; Luis, J. M. Efficient treatment of the effect of vibrations on electrical, magnetic, and spectroscopic properties. J. Comput. Chem. 2000, 21, 1572−1588. (41) Torrent-Sucarrat, M.; Solà, M.; Duran, M.; Luis, J. M.; Kirtman, B. Basis set and electron correlation effects on initial convergence for vibrational nonlinear optical properties of conjugated organic molecules. J. Chem. Phys. 2004, 120, 6346−6355. (42) Luis, J. M.; Torrent-Sucarrat, M.; Christiansen, O.; Kirtman, B. Variational calculation of static and dynamic vibrational nonlinear optical properties. J. Chem. Phys. 2007, 127, 084118. (43) Wang, B.-Q.; Li, Z.-R.; Wu, D.; Hao, X.-Y.; Li, R.-J.; Sun, C.-C. Ab initio study of the interaction hyperpolarizabilities of H-bond dimers between two π-systems. J. Phys. Chem. A 2004, 108, 2464− 2468. (44) Skwara, B.; Kaczmarek, A.; Góra, R. W.; Bartkowiak, W. On decomposition of interaction-induced electric properties of HF dimer. Chem. Phys. Lett. 2008, 461, 203−206. (45) Baranowska, A.; Fernandez, B.; Rizzo, A.; Jansik, B. The CO-Ne van der Waals complex: ab initio intermolecular potential energy, interaction induced electric dipole moment and polarizability surfaces, 3577

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation and second virial coefficients. Phys. Chem. Chem. Phys. 2009, 11, 9871−9883. (46) Baranowska, A.; Zawada, A.; Fernandez, B.; Bartkowiak, W.; Kedziera, D.; Kaczmarek-Kedziera, A. Interaction-induced electric properties and cooperative effects in model systems. Phys. Chem. Chem. Phys. 2010, 12, 852−862. (47) Baranowska, A.; Fernandez, B.; Sadlej, A. J. Importance of electron correlation effects and basis set superposition error in calculations of interaction energies and interaction-induced electric properties in hydrogen-bonded complexes: a model study. Theor. Chem. Acc. 2011, 128, 555−561. (48) Zawada, A.; Góra, R. W.; Mikołajczyk, M. M.; Bartkowiak, W. On the calculations of interaction energies and induced electric properties within the polarizable continuum model. J. Phys. Chem. A 2012, 116, 4409−4416. (49) Zawada, A.; Kaczmarek-Kędziera, A.; Bartkowiak, W. On the potential application of DFT methods in predicting the interactioninduced electric properties of molecular complexes. Molecular Hbonded chains as a case of study. J. Mol. Model. 2012, 18, 3073−3086. (50) Maroulis, G. Quantifying the performance of conventional DFT methods on a class of difficult problems: The interaction (hyper)polarizability of two water molecules as a test case. Int. J. Quantum Chem. 2012, 112, 2231−2241. (51) Avramopoulos, A.; Papadopoulos, M. G.; Sadlej, A. J. Relativistic effects on interaction-induced electric properties of weakly interacting systems: The HF-AuH dimer. J. Chem. Phys. 2002, 117, 10026−10038. (52) Baranowska-Łączkowska, A.; Fernandez, B.; Rizzo, A.; Jansik, B. New basis sets for the evaluation of the CO-Ne van der Waals complex interaction induced electric dipole moment and polarizability surfaces. Mol. Phys. 2012, 110, 2503−2512. (53) Czyżnikowska, Ż .; Góra, R. W.; Zaleśny, R.; Bartkowiak, W.; Baranowska-Łączkowska, A.; Leszczynski, J. The effect of intermolecular interactions on the electric dipole polarizabilities of nucleic acid base complexes. Chem. Phys. Lett. 2013, 555, 230−234. (54) Góra, R. W.; Błasiak, B. On the origins of large interactioninduced first hyperpolarizabilities in hydrogen-bonded π-electronic complexes. J. Phys. Chem. A 2013, 117, 6859−6866. (55) Baranowska-Łączkowska, A.; Fernandez, B.; Zaleśny, R. New basis sets for the evaluation of interaction-induced electric properties in hydrogen-bonded complexes. J. Comput. Chem. 2013, 34, 275−283. (56) Medved’, M.; Budzak, Š .; Laurent, A. D.; Jacquemin, D. Direct and indirect effects of dispersion interactions on the electric properties of weakly bound complexes. J. Phys. Chem. A 2015, 119, 3112−3124. (57) Suponitsky, K. Y.; Masunov, A. E. Supramolecular step in design of nonlinear optical materials: Effect of π−π stacking aggregation on hyperpolarizability. J. Chem. Phys. 2013, 139, 094310. (58) Fominykh, O. D.; Sharipova, A. V.; Balakina, M. Y. The choice of appropriate density functional for the calculation of static first hyperpolarizability of azochromophores and stacking dimers. Int. J. Quantum Chem. 2016, 116, 103−112. (59) Zalesny, R.; Garcia-Borràs, M.; Góra, R. W.; Medved’, M.; Luis, J. M. On the physical origins of interaction-induced vibrational (hyper)polarizabilities. Phys. Chem. Chem. Phys. 2016, 18, 22467− 22477. (60) Zaleśny, R. Anharmonicity contributions to the vibrational first and second hyperpolarizability of para-disubstituted benzenes. Chem. Phys. Lett. 2014, 595−596, 109−112. (61) Luis, J. M.; Duran, M.; Kirtman, B. Field-induced coordinates for the determination of dynamic vibrational nonlinear optical properties. J. Chem. Phys. 2001, 115, 4473−4483. (62) Bishop, D. M.; Kirtman, B.; Kurtz, H. A.; Rice, J. E. Calculation of vibrational dynamic hyperpolarizabilities for H2O, CO2 and NH3. J. Chem. Phys. 1993, 98, 8024−8030. (63) Steinmann, S. N.; Corminboeuf, C. Comprehensive benchmarking of a density-dependent dispersion correction. J. Chem. Theory Comput. 2011, 7, 3567−3577.

(64) Zaleśny, R.; Medved’, M.; Góra, R. W.; Reis, H.; Luis, J. M. Partitioning of interaction-induced nonlinear optical properties of molecular complexes. I. Hydrogen-bonded systems. Phys. Chem. Chem. Phys. 2018, 20, 19841−19849. (65) Bishop, D. M. Molecular vibration and nonlinear optics. Adv. Chem. Phys. 2007, 104, 1−40. (66) Bishop, D. M.; Hasan, M.; Kirtman, B. A simple method for determining approximate static and dynamic vibrational hyperpolarizabilities. J. Chem. Phys. 1995, 103, 4157. (67) Luis, J. M.; Martí, J.; Duran, M.; Andrés, J. L.; Kirtman, B. Nuclear relaxation contribution to static and dynamic (infinite frequency approximation) nonlinear optical properties by means of electrical property expansions: Application to HF, CH4, CF4 and SF6. J. Chem. Phys. 1998, 108, 4123−4130. (68) Eckart, C. Operator calculus and the solution of the equations of quantum dynamics. Phys. Rev. 1926, 28, 711−726. (69) Luis, J. M.; Duran, M.; Andrés, J. L.; Champagne, B.; Kirtman, B. Finite field treatment of vibrational polarizabilities and hyperpolarizabilities: On the role of the Eckart conditions, their implementation, and their use in characterizing key vibrations. J. Chem. Phys. 1999, 111, 875−884. (70) Rutishauser, H. Ausdehnung des Rombergschen Prinzips. Numerische Mathematik 1963, 5, 48−54. (71) Medved’, M.; Stachová, M.; Jacquemin, D.; André, J.-M.; Perpète, E. A. A generalized Romberg differentiation procedure for calcultion of hyperpolarizabilities. J. Mol. Struct.: THEOCHEM 2007, 847, 39−46. (72) Luis, J. M.; Duran, M.; Champagne, B.; Kirtman, B. Determination of vibrational polarizabilities and hyperpolarizabilities using field-induced coordinates. J. Chem. Phys. 2000, 113, 5203− 5213. (73) Bishop, D.; Kirtman, B. A peturbation method for calculating vibrational dynamic dipole polarizabilities and hyperpolarizabilities. J. Chem. Phys. 1991, 95, 2646−2658. (74) Bishop, D. M.; Luis, J. M.; Kirtman, B. Additional compact formulas for vibrational dynamic dipole polarizabilities and hyperpolarizabilities. J. Chem. Phys. 1998, 108, 10008−10012. (75) Frisch, M. J. et al. Gaussian 2009; Gaussian Inc.: Wallingford, CT, 2009. (76) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (77) Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (78) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (79) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchangecorrelation functional using the coulomb-attenuating method (CAMB3LYP). Chem. Phys. Lett. 2004, 393, 51−57. (80) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A long-range-corrected time-dependent density functional theory. J. Chem. Phys. 2004, 120, 8425−8433. (81) Chai, J.-D.; Head-Gordon, M. Systematic optimization of longrange corrected hybrid density functionals. J. Chem. Phys. 2008, 128, 084106. (82) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: Two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215−241. (83) Yu, H. S.; He, X.; Li, S. L.; Truhlar, D. G. MN15: A KohnSham global-hybrid exchange-correlation density functional with broad accuracy for multi-reference and single-reference systems and noncovalent interactions. Chem. Sci. 2016, 7, 5032−5051. (84) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. 3578

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579

Article

Journal of Chemical Theory and Computation (85) Ernzerhof, M.; Scuseria, G. E. Assessment of the Perdew-BurkeErnzerhof exchange-correlation functional. J. Chem. Phys. 1999, 110, 5029−5036. (86) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207−8215. (87) Heyd, J.; Scuseria, G. E. Assessment and validation of a screened Coulomb hybrid density functional. J. Chem. Phys. 2004, 120, 7274−7280. (88) Heyd, J.; Scuseria, G. E. Efficient hybrid density functional calculations in solids: Assessment of the Heyd−Scuseria−Ernzerhof screened Coulomb hybrid functional. J. Chem. Phys. 2004, 121, 1187− 1192. (89) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)]. J. Chem. Phys. 2006, 124, 219906. (90) Henderson, T. M.; Izmaylov, A. F.; Scalmani, G.; Scuseria, G. E. Can short-range hybrids describe long-range-dependent properties? J. Chem. Phys. 2009, 131, 044108. (91) Stein, T.; Eisenberg, H.; Kronik, L.; Baer, R. Fundamental gaps in finite systems from eigenvalues of a generalized Kohn-Sham method. Phys. Rev. Lett. 2010, 105, 266802. (92) Dunning, T. H., Jr. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (93) Zaleśny, R.; Baranowska-Łączkowska, A.; Medved’, M.; Luis, J. M. Comparison of property−oriented basis sets for the computation of electronic and nuclear relaxation hyperpolarizabilities. J. Chem. Theory Comput. 2015, 11, 4119−4128. (94) Jacquemin, D.; Perpéte, E. A.; Scalmani, G.; Frisch, M. J.; Kobayashi, R.; Adamo, C. Assessment of the efficiency of long−range corrected functionals for some properties of large compounds. J. Chem. Phys. 2007, 126, 144105. (95) Barone, V. Anharmonic vibrational properties by a fully automated second-order perturbative approach. J. Chem. Phys. 2005, 122, 014108. (96) Ravichandran, L.; Banik, S. Performance of different density functionals for the calculation of vibrational frequencies with vibrational coupled cluster method in bosonic representation. Theor. Chem. Acc. 2018, 137, 1−14. (97) Biczysko, M.; Panek, P.; Barone, V. Toward spectroscopic studies of biologically relevant systems: Vibrational spectrum of adenine as a test case for performances of long-range/dispersion corrected density functionals. Chem. Phys. Lett. 2009, 475, 105−110.

3579

DOI: 10.1021/acs.jctc.9b00139 J. Chem. Theory Comput. 2019, 15, 3570−3579