Performance of density functional theory in computing nonresonant vibrational (hyper)polarizabilities

A set of exchange‐correlation functionals, including BLYP, PBE0, B3LYP, BHandHLYP, CAM‐B3LYP, LC‐BLYP, and HSE, has been used to determine static and dynamic nonresonant (nuclear relaxation) vibrational (hyper)polarizabilities for a series of all‐trans polymethineimine (PMI) oligomers containing up to eight monomer units. These functionals are assessed against reference values obtained using the Møller–Plesset second‐order perturbation theory (MP2) and CCSD methods. For the smallest oligomer, CCSD(T) calculations confirm the choice of MP2 and CCSD as appropriate for assessing the density functionals. By and large, CAM‐B3LYP is the most successful, because it is best for the nuclear relaxation contribution to the static linear polarizability, intensity‐dependent refractive index second hyperpolarizability, static second hyperpolarizability, and is close to the best for the electro‐optical Pockels effect first hyperpolarizability. However, none of the functionals perform satisfactorily for all the vibrational (hyper)polarizabilities studied. In fact, in the case of electric field‐induced second harmonic generation all of them, as well as the Hartree–Fock approximation, yield the wrong sign. We have also found that the Pople 6–31+G(d) basis set is unreliable for computing nuclear relaxation (hyper)polarizabilities of PMI oligomers due to the spurious prediction of a nonplanar equilibrium geometry. © 2013 Wiley Periodicals, Inc.


Introduction
Reliable predictions of electric dipole (hyper)polarizablilities are critical for the rational design of materials possessing large nonlinear optical (NLO) response. [1,2] For a long time, it was common to investigate only electronic NLO properties leaving aside, as less important, the vibrational counterpart. Due to a large effort by several groups, however, the prominent role of the vibrational term in many instances is now well-established. Taking into account the coupling of electronic and nuclear degrees of freedom is necessary not only for resonant processes such as two-photon absorption, [3][4][5] simulation of vibronic profiles [6][7][8][9][10] and modeling of vibrational sum-frequency generation spectra, [11,12] but for non-resonant processes as well, which is the subject of this article.
In the case of the electronic component, it is possible to perform calculations using high-level correlated wavefunctions as the molecular geometry could be taken either from experiment or from more approximate methods. On the other hand, the determination of the vibrational component requires geometry optimization and the calculation of computationally demanding potential energy and electrical property surfaces. Thus, the more efficient Kohn-Sham density functional theory (KS-DFT) is an attractive alternative way to account for electron correlation effects. However, as shown by Champagne et al. [36,37] and others, [38][39][40] when conventional functionals are used KS-DFT tends to drastically overestimate electronic contributions to electric dipole (hyper)polarizabilities in spatially extended systems, especially (though not only) in P-conjugated compounds. The famous DFT overshoot problem, which has its roots primarily in the self-interaction error, [40,41] is believed to be alleviated by the recently introduced longrange corrected (LC), as well as screened, exchange-correlation (XC) functionals [42][43][44][45][46] (note that CAM-B3LYP and, e.g., LC-BLYP are both included under the LC rubric). Indeed, the recent work of Jacquemin et al., [47] who studied a series of polymethineimine (PMI) oligomers, tends to confirm this contention, at least for first hyperpolarizability. For that case, using LC-DFT they reported results that systematically approached those computed by Mïller-Plesset second-order perturbation theory (MP2) as the chain was elongated. On the other hand, in comparison with the Hartree-Fock (HF) approximation, the performance of these functionals for linear and second hyperpolarizabilities seems to be somewhat mixed. [48][49][50][51] Although the electronic contributions to resonant and nonresonant (hyper)polarizabilities have been studied extensively using LC-DFT, [49][50][51][52][53][54][55][56] there is not much known yet as far as predictions of the vibrational counterpart is concerned. [35,[57][58][59][60][61][62][63] In this study, a set of the most popular XC functionals, including recent LC and screened exchange schemes, are used to evaluate the electronic and vibrational contributions to electric dipole (hyper)polarizabilities of PMI oligomers (with number of monomer units ranging from 2 up to 8; cf. Fig. 1). These oligomers are good candidates for assessment purposes as they are known to be challenging for most XC functionals as far as electronic hyperpolarizabilities are concerned. The results of our DFT computations are compared with MP2 and CCSD reference data. Additionally, the HF functional is used to analyze the importance of electron correlation. Finally, as the conclusions presented by Jacquemin et al. [47] were drawn from calculations carried out using geometries optimized at the HF level of theory, the extent of correlation effects due to structural changes are investigated.

Methods and Computational Details
In the presence of an external electric field (F), the Cartesian component of the total molecular dipole moment l i may be expressed as a Taylor series which takes the form: [64] l where l 0 i is the i-th component of the permanent dipole moment; a ij 2x r ; x 1 ð Þ ; b ijk 2x r ; x 1 ; x 2 ð Þ and c ijkl 2x r ; x 1 ; ð x 2 ; x 3 Þ are components of the linear polarizability, first hyperpolarizability, and second hyperpolarizability, respectively; x r is the sum of the external field frequencies x i ; and K 2 ð Þ and K 3 ð Þ are factors required in order that all hyperpolarizabilities of the same order have the same static limit. For instance, in the presence of a time-dependent field, F5F 0 1F x cos xt ð Þ: [65] l Within the Born-Oppenheimer approximation, the molecular (non)linear optical properties in eq. (1) may be separated into pure electronic (P ele ) and pure vibrational (P v ) contributions, as well as the zero-point vibrational averaging correction: [2,14] where P5a; b; g. In the following subsections, a brief introduction to the computational schemes used herein to obtain these properties will be given together with a detailed description of the numerical protocols applied. The reader unfamiliar with this material can find more in the rich literature on the subject. [2,66] As far as the orientational averaging of static properties is concerned, we use: [14] a5 X i2 x;y;z f g

FULL PAPER
WWW.C-CHEM.ORG g5 X i;j f g2 x;y;z f g g iijj 5 (6) where jjljj is the norm of the dipole moment, b i 5 3 5 X j2 x;y;z f g b ijj (7) and the frequency designations have been suppressed for convenience.

Electronic contribution
In the case of DFT calculations, the electronic contribution to the (hyper)polarizability was evaluated either analytically (a, b) or fully numerically based on differentiation of the dipole moment with respect to an external electric field. In the latter event, we used the Romberg-Rutishauser algorithm (using the electric field amplitudes 62 n h, where h 5 0.0002 a.u. and n 5 0,1, …, 6). [67] At the MP2 level of theory, a semianalytical procedure was used, namely numerical differentiation of the analytical electronic polarizability was used to obtain b and c; at the CCSD level, we used the energy expansion to determine (hyper)polarizabilities in fully numerical fashion.

Vibrational contributions
Based on time-dependent perturbation theory, Bishop and Kirtman proposed a double (electrical and mechanical) perturbation theory (BKPT) treatment for the pure vibrational contributions to molecular (hyper)polarizabilities. [17] Within this approach, the pure vibrational contributions may be expressed in a short-hand "square bracket" notation as (frequency dependence is suppressed on the r.h.s. of the following equations, again for brevity of notation): Each square bracket involves products of normal coordinate derivatives of the electronic electrical properties indicated, as well as harmonic vibrational frequencies and anharmonic force constants. First derivatives of the electrical properties constitute the zeroth-order harmonic approximation; second derivatives are first-order in electrical anharmonicity, an so forth. Similarly, the quadratic terms in the potential energy expansion give rise to the mechanical zeroth-order harmonic model; cubic terms are first-order in mechanical anharmonicity, and so forth. Thus, each square bracket term may be identified by a pair of superscripts n; m ð Þ denoting the order in electrical and mechanical anharmonicity, respectively.
Within the double harmonic approximation, only four square bracket terms do not vanish. For static fields (x 1 5x 2 5x 3 50), they may be expressed as follows: where Q a ; x a denote the a-th normal coordinate and corresponding (circular) harmonic frequency, whereas P i;j;… stands for all permutations of i; j… . In this study, the diagonal components and average value of each of the above square brackets were evaluated. The derivatives of the electronic dipole moment and linear polarizability with respect to normal modes were obtained analytically together with the harmonic frequencies. At DFT and HF levels of theory, the derivatives of the linear polarizability with respect to normal modes were also computed analytically. For those methods, gradients of the electronic first hyperpolarizability were computed seminumerically, that is to say as derivatives of the electronic linear polarizability with respect to Cartesian coordinates were transformed to normal coordinates and differentiated with respect to the electric field. In doing so, the seven-point central difference algorithm was applied in the numerical differentiation and the accuracy of the procedure was compared with analytical derivatives at the HF level of theory. At MP2 level, the derivatives of the a ele and b ele with respect to the normal modes were computed numerically. At this level, a ele was obtained analytically, whereas b ele was computed by finite field differences of a ele . A check of the numerical stability of the lb ½ 0;0 ð Þ tensor computed at all levels of theory was obtained by separately computing symmetrically related components and verifying their equality, that is: As a variational alternative to the BKPT approach, it was shown by Bishop, Hasan and Kirtman (BHK) that one may use a finite-field nuclear relaxation (FF-NR) formalism to evaluate vibrational (hyper)polarizabilities approximately including some (low-order) anharmonicity. [18,20] Both static and dynamic properties are included, the latter within the infinite optical frequency approximation (squared vibrational frequencies are assumed to be negligible in comparison with squared applied field frequencies). According to this treatment: [38] P v 5P NR 1P curv 2P ZPVA (16) where P NR is the nuclear relaxation (NR) contribution. The combination P curv 2 P ZPVA , is of higher-order than the nuclear relaxation term and is neglected in this study.
As far as the NR contribution is concerned, following BHK we define: [20] where P E; R E ð Þ is an electronic property obtained at the fieldrelaxed geometry and P 0; R 0 ð Þ is the same property for fieldfree conditions. It is possible to expand the right hand side of this relation as a power series in the electric field. For example: where the expansion coefficients yield the static properties: As one may readily notice, having obtained the b 1 coefficient, together with the electronic and vibrational double harmonic (available through BKPT) terms, the first-order anharmonic vibrational contributions to the static first hyperpolarizability may be determined. Of course, one may also use b 1 directly as an estimate of the total static b. An analogous interpretation holds for g 1 , while the corresponding expansions of Da and Db yield dynamic vibrational hyperpolarizabilities in the infinite optical frequency approximation (see further below).
In applying these FF-NR formulae, one should not forget 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 optimizations can be carried out with the aid of the procedure developed by Luis et al. [19] The numerical differentiation needed to evaluate a 1 , b 1 , and g 1 was performed with the Romberg-Rutishauser algorithm (applying the electric field amplitudes 62 n h, where h 5 0.0002 a.u. and n 5 0,1, …, 6), thereby minimizing the contamination of higher power (in the field) terms. [67] Additionally, the accuracy of the FF-NR procedure was confirmed in each case by comparing l 2 ½ 0;0 ð Þ ij with the result of BKPT calculations. To reduce computational cost, the FF-NR treatment was applied only for the diagonal longitudinal (z) component of the static vibrational hyperpolarizabilities with the PMI oligomers oriented as shown in Figure 1 (terminal carbon atoms are aligned along the Cartesian z-direction).
On the other hand, for many of the vibrational dynamic properties, the double harmonic square bracket terms are sufficient to compute the nuclear relaxation contribution (there are no first-or second-order terms) within the infinite optical frequency approximation. Such properties include: the electrooptical Pockels effect (EOPE), intensity-dependent refractive index (IDRI), and electric field-induced second harmonic generation (ESHG). In square bracket notation, the average properties were determined as: [15] b NR EOPE 5 In addition, the nuclear relaxation static linear polarizability is completely determined by the double harmonic approximation [cf. eq. (11)], whereas the dynamic a NR vanishes in the infinite optical frequency approximation.

Software and density functionals
All calculations were performed with the Gaussian 09 suite of programs using default definitions of functionals (including those asymptotically corrected). [69] Geometries were optimized so that the root-mean-square (RMS) gradient was less than 10 26 Hartree/Bohr, and numerical integrations were done using a pruned (99,590) point grid. The SCF convergence was set at 10 212 on the density matrix RMS for smaller oligomers (PMI2, PMI3) and 10 210 for larger members of the set. This protocol is believed to be sufficient for description of even the low frequency modes. [70] From the plethora of available approximations to the XC energy, the set of functionals that were studied in this work includes popular representatives of the most commonly used DFT models. We used the BLYP generalized-gradient approximation (GGA) functional, together with a family of more sophisticated schemes based on it. This consists of the B3LYP [71] and BHandHLYP [69] global hybrids as well as the CAM-B3LYP[ [72] ] and LC-BLYP [43] LC hybrids. This combination allows us not only to assess the impact of exact exchange on the quality of the properties of interest, but also to gain an insight into the dependence of the vibrational (hyper)polarizabilites on the range in which the nonlocal exchange is included. The latter is particularly important, as previous studies indicated significant improvement in the quality of predictions of the electronic contributions to the molecular NLO properties with the asymptotically corrected XC potentials. Finally, we used the screened hybrid functional of Heyd-Scuseria-Ernzerhof (HSE) [73][74][75][76][77] where the nonlocal exchange energy is included only for small interelectronic distances, such that the asymptotic behaviour of the XC potential is dictated by the GGA approximation. For the sake of comparison, the parent functional of HSE, the PBE0 global hybrid [78,79] was also included.

Results and Discussion
Electronic and vibrational (hyper)polarizabilities, unless indicated otherwise, were calculated using the medium-size augcc-pVDZ basis set. This choice was dictated by the usual compromise between accuracy and computational feasibility; again, one should not overlook a more demanding protocol for calculations of the vibrational counterpart of NLO response. The recent work of Torrent-Sucarrat et al. [32] suggests that a smaller 6-311G(d) basis set might be satisfactory. However, it will be shown that the latter does not suffice for vibrational properties due to an inability to predict a sufficiently accurate potential energy hypersurface. A previous study of the basis set dependence for NLO properties of PMI was performed by Jacquemin et al. [80] However, their study dealt only with the static electronic first hyperpolarizability and, in addition, was limited to short chains (up to PMI3).
In this work, all structures optimized using the aug-cc-pVDZ basis set were found to be planar, all-trans conformers. This is in agreement with an earlier study of Medved et al. [81] On the other hand, when using the 6-311G(d) basis set, we have found in this work a spurious deviation from planarity for PMI4-PMI8.

Polarizability
Electronic and nuclear relaxation contributions to the static longitudinal polarizability, computed using the B3LYP, CAM-B3LYP, HF, MP2, and CCSD methods with the aug-cc-pVDZ basis set, are presented in Table 1 for PMI2-PMI4. These calculations on the smaller chains were done to provide a preliminary assessment of the various levels of treatment. Just two representative functionals, one hybrid (B3LYP) and one LChybrid (CAM-B3LYP), were chosen for that purpose. For PMI2, the polarizability values obtained at the MP2/aug-cc-pVTZ level (in parentheses) display very good agreement with the MP2/ aug-cc-pVDZ results, which suggests that the double-f basis is sufficient. Our data show that the relative error in the B3LYP electronic a value, in comparison with the CCSD reference, increases from PMI2 to PMI3. The same is true for the CAM-B3LYP values, although the relative errors are much smaller (roughly by a factor of two). On the contrary, although an increase in the relative error from PMI2 to PMI3 is also observed for the a NR , in that case, the B3LYP error is about half that of the CAM-B3LYP one. The performance of the HF method is surprisingly good for the electronic linear polarizability. However, HF gives the worst results of all methods for the NR contribution. In comparison with CCSD, MP2 has lower errors overall than either B3LYP or CAM-B3LYP. We have also determined the electronic and nuclear-relaxation polarizability of PMI2 at the CCSD(T)/aug-cc-pVDZ level of theory (cf. Table  1). The comparison with a el values determined using the MP2 and CCSD methods reveals that the former is slightly more successful (note that geometry differences are taken into account in our electronic values). In the case of a NR , both approaches predict equally small deviations from the reference value determined using the CCSD(T) method.
To examine the reliability of DFT calculations for the a NR of longer PMI chains, we present the average static vibrational (NR) linear polarizability of PMI2-PMI8 chains in Figure 2. In this case, we take the MP2 values as a reference because CCSD calculations are not feasible for the longer chains. Despite the arguments made above to justify this choice any conclusions that rely on MP2 as the reference should be regarded with caution. As one may readily notice in Figure 2 has the best agreement with MP2 (average error 5 11.6%). HF also performs very well for the a NR ; the average difference from MP2 is 13.5%. After CAM-B3LYP, the DFT approaches that agree Table 1. Electronic and nuclear relaxation contributions to the longitudinal components of static electric-dipole (hyper)polarizabilities (a zz , b zzz , and c zzzz ). First hyperpolarizability Table 1 contains a summary of static NR longitudinal (hyper)polarizabilities for the smaller oligomers (n 5 2-4) obtained using the aug-cc-pVDZ basis set. For comparison purposes, the corresponding electronic (hyper)polarizabilities are included as well. As in the case of the linear polarizability, the small differences between the aug-cc-pVDZ and aug-cc-pVTZ (in parentheses) first hyperpolarizabilities for PMI2 support the choice of the former basis set for calculations for larger oligomers. We begin here with a brief discussion of the electronic first hyperpolarizability. Our MP2 results in Table 1 and Figure 3 are in agreement with those of Jacquemin, et al. [47] For short chains (PMI2 and PMI3), the longitudinal component is overestimated (in magnitude) at the MP2 level of theory as compared to the CCSD (and CCSD(T) for PMI2) reference. The percent difference decreases significantly; however, on going from PMI2 to PMI3. As opposed to HF, B3LYP, and CAM-B3LYP, it is noteworthy that MP2 correctly predicts the negative sign of the electronic longitudinal first hyperpolarizability for PMI2 and PMI3.
For longer chains, average values of the static electronic first hyperpolarizability calculated by means of the seven DFT schemes studied here, as well as the HF and MP2 methods, are presented in Figure 3. In agreement with the results of Table 1, all DFT methods and HF predict that b ele has the opposite sign to that of the reference MP2 value. These methods also predict that the average values change in the opposite direction from MP2 as a function of increasing chain length (the sign of b ele is opposite to that of the longitudinal b ele because l z is negative [see eq. (5)]. The only exception is LC-BLYP; in that case, b ele initially becomes more negative for small chains but begins to increase for PMI7 and PMI8. The general behavior just described is quite similar to that found by Jacquemin, et al. over the range of chain lengths considered here, although geometry and basis set considered by these authors are different. It is only for longer chains that they find CAM-B3LYP and B3LYP begin to approach the MP2 value (in their case for the longitudinal electronic first hyperpolarizability). [47] We have examined the effect of geometry, for the PMI2-PMI6 oligomers using the CAM-B3LYP functional at the MP2 optimized geometries. Note in Figure 3, the rather insignificant change in the property value. Thus, the sign disagreement between the DFT and MP2 descriptions of the electronic b for  PMI n oligomers (n 5 2-6 and, possibly, n 5 7-8) is not due to the choice of geometry. Now we turn to the NR vibrational first hyperpolarizabilities. The comparison between MP2 and CCSD in Table 1 shows good agreement between both methods, especially for PMI3. This suggests that MP2 is a reasonable choice as the reference for NR first hyperpolarizability values. On the contrary, HF, B3LYP, and CAM-B3LYP strongly underestimate the static longitudinal b NR (as compared to CCSD) with relative errors between 53 and 134%. The choice of MP2 is further supported by the results of b NR calculations for PMI2 at the CCSD(T)/ aug-cc-pVDZ level of theory. The jb NR MP2 j are about one order of magnitude smaller than jb NR As far as b NR for longer chains is concerned, the EOPE is simplest to deal with because the infinite optical frequency approximation (IOFA) causes the anharmonic contribution to vanish, thereby leaving only the double harmonic term [cf. eq. (22)]. It happens that the behavior of the DFT functionals bears a similarity to what was observed for a NR (cf. Fig. 4 For the static NR first hyperpolarizability, there is also a firstorder anharmonic contribution. In carrying out the required field-dependent geometry optimizations, these calculations proved to be very expensive using the MP2 method for the longer PMI molecules. Thus, our conclusions regarding the DFT performance are limited to the PMI2-PMI4 values presented in Table 1. Contrary to the electronic static longitudinal b, the MP2 and CCSD values for the corresponding NR property are quite similar to one another for both PMI2 and PMI3. All three oligomers in the table give B3LYP and CAM-B3LYP results that are roughly 1.5-2.0 times larger in magnitude than the corresponding MP2 value. The HF value is even somewhat larger in magnitude. It is clear from the table that the double harmonic term, la ½ 0;0 ð Þ zzz , makes the dominant contribution for all methods and is the major source of the difference between MP2 and the DFT values. The DFT electrical anharmonicity contribution amplifies this discrepancy, whereas the larger mechanical anharmonicity term improves the agreement. Table 1 also contains the property values determined for PMI4 at the MP2/6-311G(d) level of theory. At that level, the values of b NR (and c NR as well) are strikingly different from the MP2/aug-cc-pVDZ results. As noted above, using the 6-311G(d) basis set, nonplanar field-free equilibrium geometries are obtained for PMI4-PMI8. Furthermore, the field-dependent geometry optimizations, even in a weak-field regime, lead to highly distorted structures not observed when using the augcc-pVDZ basis set. Such spurious nonplanar geometries predicted by the 6-311G(d) basis set lead to longitudinal b NR values with the wrong sign and a strongly overestimated longitudinal c NR . These strange finite-field NR results were checked using analytical BKPT formulas and the field-induced coordinate method, [82,83] which gave almost identical values of b NR (within the range of numerical error estimated as 65 a.u.). The inability of the 6-311G(d) basis set to predict, even qualitatively, the anharmonic contribution to b NR is probably due to the fact that members of the family of basis sets developed by Pople and coworkers, if not well-balanced, will tend to give a poor description of potential energy hypersurfaces. [84,85] Second hyperpolarizability Finally, we come to the second hyperpolarizability. Two main conclusions can be drawn from our results for the static electronic longitudinal second hyperpolarizability presented in Table 1). First, electron correlation effects are important. For both PMI2 and PMI3, the relative difference between the HF value and the CCSD reference is 36%. A similar conclusion was drawn by Medved' et al. for PMI chains up to PMI6. [81] For PMI2, there is an excellent agreement between c NR value computed at CCSD(T)/aug-cc-pVDZ level of theory and the data obtained at the MP2 and CCSD levels. Second, as the chain is elongated, the value predicted by the B3LYP functional becomes worse as compared to the CCSD or MP2 method (for PMI2, PMI3, and PMI4, the relative differences with respect to the MP2 (CCSD) values are 6% (21%), 25% (49%), and 63%, respectively). As we have come to expect, the Coulomb-attenuated correction (CAM-B3LYP) leads to a substantial improvement of the B3LYP predictions. In this case, the PMI2, PMI3, and PMI4 relative differences with respect to MP2 (CCSD) are reduced to 2% (11%), 14% (5%) and 6%, respectively. Although there is a significant difference between MP2 and CCSD, the same trends are revealed in either comparison. Once more, we add the proviso that the above statements may need to be revised when longer chains are examined.
As far as the NR vibrational contribution to the longitudinal second hyperpolarizabilities are concerned, it is again convenient to begin with the dynamic properties for which there is no anharmonicity correction in the IOFA, namely c IDRI and c ESHG . In the case of the average c IDRI , the BLYP, B3LYP, HSE, and PBE0 functionals (see Fig. 5) strongly overestimate this property in comparison with MP2 and, what is worse, the relative deviation rapidly increases with increasing chain length. Again the BLYP values show the poorest agreement; the relative difference with regard to the MP2 values is 71.9 and 202.3% for PMI2 and PMI7, respectively. As in the case of the other double harmonic properties treated thus far, the CAM-B3LYP functional gives the best agreement. In fact, discounting the two smallest oligomers, the relative difference from MP2 is always smaller than 10.1%, with an average difference of only 5.2%. The other two DFT approaches that differ least from MP2 are BHandHLYP and LC-BLYP, with average relative differences of 11.5 and 16.0%, respectively.
The scenario is clearly worse for the ESHG process. As one can see from Figure 6, all DFT functionals, and the HF approximation, strongly underestimate the MP2 value. For chains longer than PMI2, they even fail to predict the correct sign. To further validate the MP2 method as the reference, the c ESHG value was computed for PMI3 at the CCSD/aug-cc-pVDZ level of theory. The latter was found to be 537 a.u., which may be compared with the MP2 value of 773 a.u. This validates the sign of the MP2 result. Because the NR IOFA for c ESHG involves the first derivative of the electronic b with respect to normal coordinates, [see eq. (14)], the erroneous sign is related to the result, discussed above, that neither HF nor any of the DFT methods used here reproduces the proper sign of b ele (see  Finally, the NR static second hyperpolarizability is the most difficult property to determine. It contains both first-and second-order anharmonicity contributions. From the results reported in Table 1, the longitudinal c NR (0;0,0,0) is found to be about the same for both ab initio methods as well as the DFT treatment with the B3LYP and CAM-B3LYP functionals (for PMI2 and PMI3 the relative differences of all methods considered with respect to the CCSD values are smaller than 9.5%). There are two square bracket terms that combine to constitute the double harmonic contribution. The a 2 ½ 0;0 ð Þ term obtained for the DFT functionals overshoots the MP2 value, although the ratio is less than about 1.5 with CAM-B3LYP and decreases to 1.1 in that case for PMI4. On the other hand, the lb ½ 0;0 ð Þ term obtained with the DFT functionals has the wrong sign. This is undoubtedly a reflection of the fact that the static electronic b has the wrong sign as discussed above. It is clear that the satisfactory performance of HF, B3LYP, and CAM-B3LYP for the total c NR 0; 0; 0; 0 ð Þis due to a fortuitous error cancellation where the anharmonic corrections are substantially overestimated and alleviate the poor performance of the lb ½ 0;0 ð Þ term. As judged by the values determined using the MP2 method, the first-and second-order anharmonic contributions to the static nuclear relaxation second hyperpolarizability are substantially more important than the double harmonic contribution. By comparing HF with MP2 and CCSD, it might seem that correlation effects are relatively unimportant for the static c NR . However, as pointed out above, the good performance of HF for the total value is just the fruit of a fortunate cancellation of errors. Indeed, correlation is essential to accurately reproduce both the harmonic and anharmonic contributions to the static c NR . Finally, we note that the static c NR computed at the MP2/6-311G(d) level for PMI4 is 20 times larger than the MP2/aug-cc-pVDZ value. As discussed above, the spurious nonplanar equilibrium geometry obtained with the 6-311G(d) basis set is the origin of this discrepancy.

Summary and Conclusions
We have used seven of the most popular XC functionals, within the Kohn-Sham formulation of density functional theory, to assess their performance in computing nonresonant NR vibrational (hyper)polarizabilities for a series of (PMI) n oligomers (n 5 2-8). The properties studied included not only static (hyper)polarizabilities but also several dynamic (hyper)polarizabilities within the infinite optical frequency approximation. For reference purposes, the CCSD and MP2 (for n > 3) methods were used together with the aug-cc-pVDZ basis set. An attempt to use the smaller 6-311G(d) basis to compute NR hyperpolarizabilities proved to be unreliable because, for n > 3, the optimum geometry is erroneously calculated to be nonplanar.
The performance of the various functionals was found to be somewhat mixed. CAM-B3LYP yields the best overall results for the DFT functionals considered in this paper. It is the best as far as a NR 0; 0 ð Þ, c NR IDRI and c NR 0; 0; 0; 0 ð Þ are concerned and gives good agreement with the reference values. For c NR 0; 0; 0; 0 ð Þ , however, this success is due to a cancellation of errors between (smaller) harmonic and (larger) anharmonic contributions to the total result. As far as ESHG is concerned, in contrast to the other second hyperpolarizabilities, none of the functionals, or the HF method, predict the correct sign of the static NR vibrational c. This behavior is related to the inability of all functionals to predict the correct sign of the static electronic b which, in turn, agrees with previous reports. With regard to first hyperpolarizabilities, although the CAM-B3LYP functional does not give the best agreement with the reference value for the NR Pockels effect, it, nonetheless, does not differ substantially from that value. On the other hand, the CAM-B3LYP static b NR significantly overshoots the reference value, at least for the small (n 5 2-4) chains that were studied. In this case, an overshoot occurs, but for other properties the error corresponds to an undershoot and, thus, there is no systematic pattern. The effect of anharmonicity is quite small for the static b NR but, on the contrary, it is dominant for the static NR second hyperpolarizability.