Can Density Functional Theory Be Trusted for High-Order Electric

Subscriber access provided by UNIV AUTONOMA DE COAHUILA UADEC

Quantum Electronic Structure

Can Density Functional Theory Be Trusted for High-Order Electric Properties? The Case of Hydrogen-Bonded Complexes. Robert Zalesny, Miroslav Medved, Sebastian Sitkiewicz, Eduard Matito, and Josep María Luis J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00139 • Publication Date (Web): 13 May 2019 Downloaded from http://pubs.acs.org on May 13, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Can Density Functional Theory Be Trusted for High-Order Electric Properties? The Case of Hydrogen-Bonded Complexes. Robert Zaleśny,∗,† Miroslav Medveď,‡ Sebastian Sitkiewicz,¶ Eduard Matito,∗,§ and Josep M. Luis∗,k †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 kInstitute of Computational Chemistry and Catalysis and Department of Chemistry, University of Girona, Campus de Montilivi, 17003 Girona, Catalonia, Spain E-mail: [email protected]; [email protected]; [email protected] 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

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

studied properties, the average absolute errors below 20% can only be obtained using the CAM-B3LYP 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 range-separation parameter µ for 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 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.

2


Page 2 of 36

Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 becomes often 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) 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 cost-effective 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 error 17 and the wrong exchange potential decay of generalized gradient approximations (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 zwitter-ionic 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 of3



ten an error function that depends of 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 non-resonant 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 towards more complex systems, the electronic NLO properties of various types of weakly-bound systems including hydrogen-bonded and van der Waals complexes were intensively studied by high-level postHartree-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

4


Page 4 of 36

Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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). Let us recall that the non-resonant vibrational contributions to hyperpolarizability are in many cases of the same magnitude or even larger than their electronic counterparts. This, in general, might hold for static electric 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 study 59 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

5



Page 6 of 36

Theory and Computational Details Under the Born-Oppenheimer approximation, one may define the total property P as the sum of electronic (P el ), nuclear relaxation (P nr ) and curvature (P curv ) 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 P nr 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 (β) and second (γ) hyperpolarizability. As far as the nuclear relaxation contribution is concerned, following BHK we define: 66

µi (F, RF ) − µi (0, R0 ) =

X

a1ij Fj +

j

1X 1 1X 1 bijk Fj Fk + g Fj Fk Fl + . . . 2 jk 6 jkl ijkl

(2)

where µ(F, RF ) is an electronic property obtained at the field-relaxed geometry and µ(0, R0 ) is the same property at field-free conditions. The expansion coefficients yield the static properties: (0,0)

el nr el a1ij = αij + αij = αij + [µ2 ]ij

(0,0)

(3) (1,0)

(0,1)

nr el el b1ijk = βijk + βijk = βijk + [µα]ijk + [µ3 ]ijk + [µ3 ]ijk

(4)

el nr 1 + γijkl gijkl = γijkl (0,0)

(5) (0,0)

(1,0)

(0,1)

(2,0)

(0,2)

(1,1)

el = γijkl + [α2 ]ijkl + [µβ]ijkl + [µ2 α]ijkl + [µ2 α]ijkl + [µ4 ]ijkl + [µ4 ]ijkl + [µ4 ]ijkl

Each square bracket involves the products of normal coordinate derivatives of the electronic electrical properties indicated, as well as harmonic vibrational frequencies and anharmonic

6


Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 field-dependent 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 =

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

(6)

where k is related to the strength of the field used in the FF-NR calculations (see below). The value of P 0,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 a.u. In all calculations employing the Romberg– Rutishauser scheme we used the following sequence of fields ±Fi , ±2Fi , . . . , ±2k Fi , where k=7. All FF-NR results were re-checked using field-induced coordinates 61,72 (FIC) and Bishop-Kirtman perturbation theory 73,74 (BKPT) methods. In 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)

7



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 LC-BLYP (µ=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 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 to obtain 1% accuracy in predictions of electronic and nuclear relaxation (hyper)polarizabilities of small molecules. 93

Results and Discussion As highlighted in 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-cc-pVTZ level of theory. For the sake of computing various anharmonic terms,

8


Page 8 of 36

Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 Fig. 1. Note that, for the sake of clarity, errors which amount to 100% and larger are marked with the same colour (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 in-depth analysis of the vibrational counterparts. All the underlying numerical data are included in the Supporting Information. The average and maximum absolute relative errors in predicting αel , β el and γ el are shown in Figures 2a, 2b and 2c, 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%. 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 meta-GGA 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 Supporting Information contains the comparison of

9



electronic properties computed with the set of DFAs at the equilibrium geometries obtained using CCSD(T) method. The comparison of the results presented in Table S9 with their counterparts presented in Figure 2 allows to draw the conclusion that the choice of reference geometry has 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, 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 worse 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. 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 and anharmonic terms to γ nr at the CCSD(T) level are unfeasible, and then this approach cannot be used to decompose the nuclear relaxation contributions. Opportunely, the accuracy of 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

10


Page 10 of 36

Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of some of the DFAs. Figure 4 contains the net harmonic and 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 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 can not be predicted only by analyzing the effect of the anharmonicity in the vibrational frequencies and intensities (see Table S14 in 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: n nr nr (DFA) − γi,har (MP2) 1 X γi,har ∆har (DFA) = × 100% n i=1 γinr (MP2)

(7)

n nr nr X γ (DFA) − γ (MP2) 1 i,anh i,anh ∆anh (DFA) = × 100% nr n γi (MP2)

(8)

i=1

∆har and ∆anh are shown in Figure 5. The analysis of these figures suggests that the huge errors in predicting γ nr by M06, M06-2X and ωB97X (see Fig. 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

11



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 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 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 perturbation theory formalism, this term encompasses a40 ∂Q ∂Q ∂Q ∂Q ∂µ ∂µ ∂µ ∂µ and a30 a30 ∂Q contributions, where a30 and a40 denote third- and fourth-order energy ∂Q ∂Q ∂Q

derivatives with respect to normal modes. In an attempt to pinpoint the origins of large error in [µ4 ](0,2) term for 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 [µ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 these findings are in line with the results presented herein and work by Biczysko et al. 97

The satisfactory performance of the LC-BLYP functional (for each of studied properties the relative errors do not exceed 30%) encouraged us to analyze the effect of tuning the rangeseparation 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

12


Page 12 of 36

Page 13 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 µ 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 tuning of µ 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 CCSD(T) reference increasing by 3-15% with respect to LC-BLYP results. The tuning of µ parameter, by and large, does not improve significantly the predictions of γ el , i.e. for HCN· · · HCl and HCN· · · HCN there is only insignificant decrease in relative error while for HCN· · · HF we find 10% increase in relative error with respect to the CCSD(T) reference. In the case of nuclear-relaxation contributions, the pattern of changes upon using system-specific µ 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 µ 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 dispersion-corrected functional, LC-BLYP-dDsC. These preliminary results indicate that electronic and vibrational hyperpolarizabilities are not much improved using dispersion-

13



corrections, at least for the systems considered in the present work.

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 is usually 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 CAM-B3LYP 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, M062X), 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 range-separation parameter µ for 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 non-resonant (hyper)polarizabilities of hydrogen-bonded complexes, we do not

14


Page 14 of 36

Page 15 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 towards thorough assessment of density functional theory in predicting high-order electric properties including anharmonic vibrational corrections.

Supporting Information Available The numerical data used to perform statistical analysis discussed in the manuscript. The dependence of longitudinal electronic (hyper)polarizabilities on the geometry of selected molecular complexes. Selected anharmonic contributions to nuclear-relaxation second hyperpolarizability computed using the M06, M06-2X and CAM-B3LYP functionals.

Acknowledgement 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) and CMST COST Action CM1405 MOLIM: MOLecules In Motion. J.M.L. and E.M. are grateful for financial support from the Spanish MINECO CTQ2014-52525-P (E.M. and J.M.L.) 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). Authors gratefully acknowledge Wroclaw Centre for Networking and Supercomputing for the generous allotment of computer time.

15



References (1) Bella, S. D.; 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 H-bonded 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. 2014, 10, 2114–2124. (5) Medveď, M.; Budzák, Š.; Bartkowiak, W.; Reis, H. Solvent effects on molecular electric properties; Handbook of Computational Chemistry; Springer International Publishing, 2017; Vol. 17; pp 741–795. (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) Cacheiro, J. L.; 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.

16


Page 16 of 36

Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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 coupled-cluster hyperpolarizabilities. J. Chem. Phys. 2009, 130, 194108. (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

17



Page 18 of 36

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 1994, 49, 2421. (19) Wannere, C. S.; Sattelmeyer, K. W.; Schaefer III, H. F.; 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 longrange 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.

18


Page 19 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(23) Savin, A. A combined density functional and configuration interaction method. Int. J. Quantum Chem. 1988, 34, 59–69. (24) Savin, A. Density functional methods in chemistry; Springer, 1991; pp 213–230. (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.; Medveď, 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 uthi, 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 donor-acceptor-substituted azobenzene. J. Chem. Phys. 2010, 133, 244308. (29) Medveď, 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

19



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. Ann. 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.

20


Page 20 of 36

Page 21 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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. Comp. 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, and second virial coefficients. Phys. Chem. Chem. Phys. 2009, 11, 9871–9883. (46) Baranowska, A.; Zawada, A.; Fernandez, B.; Bartkowiak, W.; Kedziera, D.; Kaczmarek-

21



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 interactioninduced 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 interaction-induced electric properties of molecular complexes. Molecular H-bonded 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. Quant. Chem. 2012, 112, 2231–2241. (51) Avramopoulos, A.; Papadopoulos, M. G.; Sadlej, A. J. Relativistic effects on interactioninduced 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

22


Page 22 of 36

Page 23 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 interaction-induced 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. Comp. Chem. 2013, 34, 275–283. (56) Medveď, 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. 2015, 116, 103–112. (59) Zalesny, R.; Garcia-Borr‘as, M.; G’ora, 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. 23



(62) Bishop, D. M.; Kirtman, B.; Kurtz, H. A.; Rice, J. E. Calculation of vibrational dynamic hyperpolarizabilities for H2 O, CO2 and NH3 . J. Chem. Phys. 1993, 98, 8024–8030. (63) Steinmann, S. N.; Corminboeuf, C. Comprehensive benchmarking of a densitydependent 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 interactioninduced 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. 1998, 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.

24


Page 24 of 36

Page 25 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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 1988, 38, 3098–3100. (77) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 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 exchange-correlation functional using the coulomb-attenuating method (CAM-B3LYP). 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.

25



(81) Chai, J.-D.; Head-Gordon, M. Systematic optimization of long-range 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 Kohn-Sham 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. (85) Ernzerhof, M.; Scuseria, G. E. Assessment of the Perdew-Burke-Ernzerhof exchangecorrelation functional. J. Chem. Phys. 1999, 110, 5029–5036. (86) J. Heyd, G. E. Scuseria and M. Ernzerhof, 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.

26


Page 26 of 36

Page 27 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(90) T. M. Henderson, A. F. Izmaylov, G. Scalmani and G. E. Scuseria, 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 Jr., T. H. 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.; Medveď, 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 longrange/dispersion corrected density functionals. Chem. Phys. Lett. 2009, 475, 105–110.

27



Page 28 of 36

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

γ el

αnr

β nr

γ nr

26.10 23.21

43.98 42.57

14.51 12.86

15.93 12.39

20.96 25.62

12.00 20.00

HCN· · · HCl OLC-BLYP (µ=0.4502) LC-BLYP (µ=0.4700)

-0.17 -0.83

20.25 16.96

-14.47 -17.11

21.93 19.82

HCN· · · HCN OLC-BLYP (µ=0.4553) LC-BLYP (µ=0.4700)

-1.82 -2.22

35.71 32.50

-18.46 -20.00

14.01 12.95

HCN· · · HF OLC-BLYP (µ=0.5449) LC-BLYP (µ=0.4700)

-3.24 -1.04

-64.86 -51.35

28

-26.67 -16.67


12.03 15.41

Page 29 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 colour.

29



Figure 2: Average and maximum absolute relative error in predicting electronic polarizability (a), first (b) and second (c) hyperpolarizability for a set of 9 molecular complexes with respect to CCSD(T)/aug-cc-pVTZ reference values. 30 ACS Paragon Plus Environment

Page 30 of 36

Page 31 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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


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

32


Page 32 of 36

Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: Average and maximum absolute relative error 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 definition of the relative error.

33



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.

34


Page 34 of 36

Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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

35



TOC Graphic

36


Page 36 of 36

Can Density Functional Theory Be Trusted for High-Order Electric

Recommend Documents