On the physical origins of interaction-induced vibrational (hyper)polarizabilities

This paper presents the results of a pioneering exploration of the physical origins of vibrational contributions to the interaction-induced electric properties of molecular complexes. In order to analyze the excess nuclear relaxation (hyper)polarizabilities, a new scheme was proposed which relies on the computationally efficient Bishop–Hasan–Kirtman method for determining the nuclear relaxation contributions to electric properties. The extension presented herein is general and can be used with any interactionenergy partitioning method. As an example, in this study we employed the variational–perturbational interaction-energy decomposition scheme (at the MP2/aug-cc-pVQZ level) and the extended transition state method by employing three exchange–correlation functionals (BLYP, LC-BLYP, and LC-BLYP-dDsC) to study the excess properties of the HCN dimer. It was observed that the first-order electrostatic contribution to the excess nuclear relaxation polarizability cancels with the negative exchange repulsion term out to a large extent, resulting in a positive value of Da due to the contributions from the delocalization and the dispersion terms. In the case of the excess nuclear relaxation first hyperpolarizability, the pattern of interaction contributions is very similar to that for Da, both in terms of their sign as well as relative magnitude. Finally, our results show that the LC-BLYP and LC-BLYP-dDsC functionals, which yield smaller values of the orbital relaxation term than BLYP, are more successful in predicting excess properties.


Introduction
Great attention has been paid to the materials characterized by a large nonlinear optical response, as they are commonly used in optical rectification, modern waveguides, and optical communication technology and for optical switching and signal processing. 1 At the molecular level, the linear and nonlinear optical properties are given by the linear polarizability (a) and first (b) and second (g) hyperpolarizabilities, respectively.3][4][5] The latter effect can be attributed to the interaction-induced (excess) properties, which are defined as the difference between a property of the complex and the sum of the properties of the noninteracting subsystems.The interest in excess properties and collision-induced spectra has evolved over the last two decades from pairs of interacting atoms [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] to complexes involving molecules.  An eensive body of literature exists on the purely electronic interactioninduced electric properties, including various aspects of their theoretical treatment.Also, several methods were proposed for partitioning the excess electric properties on the grounds of intermolecular perturbation theory, as it leads to a better understanding of their origins. 7,10The most comprehensive approach was proposed by Moszyn ´ski et al., 10 who formulated a symmetryadapted perturbation theory of interaction-induced electric properties.A similar approach was followed by other authors, who carried out an exploration of the contributions to purely electronic excess properties. 23,27,39,40Below, we shall briefly summarize the most relevant findings.
Skwara et al. 23 performed the decomposition of the excess electric properties (dipole moment, polarizability and first hyperpolarizability) of the linear HF dimer into interaction energy components.They highlighted a subtle balance of the electrostatic and exchange repulsion contributions.Go ´ra et al. studied the excess electric properties of a set of hydrogen-bonded systems, encompassing quasi-linear dimers of hydrogen cyanide, urea, diformamide, 4-pyridone, 4-nitroaniline, and the complex of hydrogen fluoride with nitroacetylene. 27,40These authors found that the origins of the interaction-induced dipole moment are purely electrostatic in nature with only a minor contribution from the charge-delocalization term and negligible remaining components.On the other hand, the origins of interactioninduced polarizabilities were slightly more complicated even though they also seem to be dominated by the first-order electrostatic term. 40In particular, in some of the studied complexes with bifurcated hydrogen bonds the exchange-repulsion effects were found to quench completely the associated first-order electrostatic term.Thus, the excess polarizabilities of smaller complexes were shown to be dominated by induction (polarization) and exchange-induction effects. 40In the case of excess first hyperpolarizability Go ´ra and Błasiak reported that its origin was much more system-dependent than other electric properties.For example, in the case of the urea dimer, the excess first hyperpolarizability was found to be diminished due to induction and exchangeinduction terms, whereas in diformamide and 4-pyridone dimers the observed reduction was primarily due to exchange-repulsion effects. 40Even the dispersion contribution, which is in most cases negligible, may sometimes become unusually large, e.g. it contributes more than 15% to the excess electronic first hyperpolarizability of the p-aminobenzoic acidÁ Á Ábenzene complex. 43xcess electronic polarizabilities for a large set of nucleic acid base pair conformations were studied by Czyz ˙nikowska et al. 39 It turned out that the physical origins of this property are quite different for stacked and hydrogen-bonded base pairs. 39In the case of the stacked cytosine dimers the leading contribution to excess polarizability was due to the electrostatic contribution.The interaction pattern for Watson-Crick guanine-cytosine and adenine-thymine pairs was found to be different, i.e. the excess polarizability was dominated by exchange and induction interactions and electrostatics played only a minor role. 39The effects of many-body interactions on excess properties were also thoroughly examined. 22,27inally, it should not be overlooked that several approaches to estimate partial contributions of atoms or molecular fragments to (hyper)polarizabilities have also been developed.5][46] Their approach allows the analysis of local induced polarizabilities and charge flow between different regions of a molecule.Ye and Autschbach have investigated local contributions to the first hyperpolarizability of donor-acceptor substituted molecules by means of the Mulliken type and natural bond orbital (NBO) analyses. 47Their results highlight the contributions from atoms and bonds on different functional groups to the total value of the first hyperpolarizability.9][50] In their analysis the hyperpolarizabilities are decomposed into local and charge-transfer parts.A quite different approach is the concept of the (hyper)polarizability densities developed by Nakano and coworkers, [51][52][53] which are defined as the derivatives of the Mulliken charge densities with respect to the electric field.Their (hyper)polarizability density plots allow the visualization of the local contributions to the total properties.
At this point it is important to highlight that rather insufficient attention was paid to the vibrational contributions to interaction-induced electric properties.A preliminary treatment of this subject can be found in the paper of Avramopoulos et al., 54 who examined the influence of relativistic effects on the interaction-induced dipole moment and polarizability of the HF-AuH complex including the effect of molecular vibrations.These authors demonstrated, based on the low-order perturbation method of Bishop and Kirtman (BKPT), 55 that the vibrational contribution to the induced polarizability is approximately an order of magnitude larger than its electronic counterpart. 54It should be noted, though, that they argued that this finding may significantly depend on approximations used in the low-order perturbation treatment.
Our aim in this work is to perform a pioneering exploration of the physical origins of nuclear relaxation contributions to the excess electric properties of molecular complexes.The HCN dimer is chosen as a model system for this purpose, mainly because its electronic excess properties were already thoroughly studied in our previous paper. 27Theory and computational details This section is organized as follows.First, we outline the computational treatment of vibrational (hyper)polarizabilities. Next, we describe a framework for the partitioning of the excess vibrational (hyper)polarizabilities.Finally, we present the interaction energy decomposition schemes employed in this work.

Vibrational (hyper)polarizabilities
Within the Born-Oppenheimer approximation, one may define the total property P as the sum of electronic (P e ), nuclear relaxation (P nr ) and curvature (P curv ) contributions: 56 P = P e + P nr + P curv (1) P nr and P curv arise from the change of the electronic and vibrational energies caused by the field-induced relaxation of the equilibrium geometry, respectively.The nuclear relaxation contributions contain all the leading terms of the BKPT formulas, including the anharmonic corrections of the lowest order.The curvature contributions, which include the rest of the high-order anharmonicity, are far more computationally expensive than P nr and they are not computed here.Bishop, Hasan and Kirtman (BHK) demonstrated that one may employ a computationally efficient finite-field nuclear relaxation (FF-NR) formalism to evaluate the nuclear relaxation contribution to (hyper)polarizabilities. 57,58 In the present study, we confine the analysis solely to static properties, including polarizability and first hyperpolarizability.As far as the nuclear relaxation contribution is concerned, following BHK we define: 57 where m(F,R F ) is the dipole moment obtained at the field-relaxed geometry and m(0,R 0 ) is the same property for field-free conditions.The expansion coefficients yield the static properties: When applying these FF-NR formulas, one should not overlook that the geometry relaxation must not include the rotations of the molecule through the 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).Therefore, the field-dependent optimization must be performed strictly maintaining the Eckart conditions. 59Such optimizations can be carried out with the aid of the procedure developed by Luis et al. 60

Excess vibrational (hyper)polarizabilities and their partitioning
Let the equilibrium geometry of a complex composed of A and B in an applied electric field, F, be denoted by AB F , while E(F,A AB,F ) and E(F,B AB,F ) are the field-dependent energy of A and B at the geometry corresponding to the field-relaxed geometry of AB F .AB stands for the equilibrium geometry of the complex without the field present.The field-dependent interaction energy at field-relaxed geometry AB F is thus given by: The derivative of DE int (F,AB F ) with respect to electric field components along Cartesian directions i, j,. . . is given by: For n 4 1 the first term on the r.h.s. of eqn (6) can be recognized as the sum of the static electronic (hyper)polarizability (P e ) and the nuclear relaxation (P nr ) of a complex AB: The two remaining terms on the r.h.s of eqn (6) can also be split in the electronic and nuclear-relaxation contributions to a property P of subsystems A and B, respectively.Henceforth, we will refer to P nr (A AB ) and P nr (B AB ) as pseudo-nuclear relaxation (pseudo-NR) contributions to monomer properties.We add the prefix pseudo because these monomeric contributions are computed at geometries corresponding to interacting AB, instead of the equilibrium geometry of A and B, respectively: Combining eqn (6) with eqn (7-9) yields: where Since DE int can be represented as the sum of interaction energy components, E int,X , it is possible to re-write eqn (10): Finally, the above equation can be presented in a form that allows us to directly analyze the DP nr contribution in terms of interaction-energy components DE int,X :

Interaction energy partitioning schemes
In the preceding section we have presented a framework for the partitioning of nuclear-relaxation (hyper)polarizabilities.However, such partitioning is general and it still requires a particular choice of interaction energy decomposition scheme.Since the seminal work of Morokuma, who was the first to introduce a robust scheme for the decomposition of supermolecular interaction energy into physically meaningful components based on analysis and interpretation of the Fock matrix elements, 61 many similar approaches have been proposed including the extended transition state (ETS) method of Ziegler and Rauk, 62 the improved Kitaura-Morokuma (KM) scheme, 63 the constrained space orbital variation (CSOV), 64 the reduced variational space self-consistent field (RVS-SCF), 65 the natural energy decomposition analysis (NEDA), 66 the block-localized wavefunction approach (BLW-ED) 67 and the method recently proposed by Su and Li. 68here is also a variational-perturbational decomposition scheme (VP-EDS) [69][70][71][72] which, like the above schemes, is based on the partitioning of the supermolecular interaction energy.However, in this scheme the interpretation of the obtained components refers to the intermolecular perturbation theory. 73][71][72] In this paper we decompose electronic and vibrational contributions to electric excess properties using the ab initio VP-EDS scheme and the ETS method within the Kohn-Sham density functional theory (KS-DFT) framework.Although the VP-EDS and KS-DFT schemes were previously used for the partitioning of interaction-induced properties, 23,27,39,40 the reliability of the latter approach is largely determined by the choice of exchange-correlation potential.It can be argued that the excess property partitioning imposes additional requirements regarding the choice of density functional.Unfortunately, there is not much known about the performance of density functionals in computations of excess properties and available data are to some extent contradictory.Zawada et al. studied

the purely
This journal is © the Owner Societies 2016 electronic contributions to the interaction-induced electric properties by employing several exchange-correlation functionals, including B3LYP, LC-BLYP, PBE0, M06-2X and CAM-B3LYP, and concluded that their overall performance is satisfactory, at least for hydrogen-bonded chains consisting of HF, H 2 CO and NH 3 molecules. 29More pessimistic conclusions regarding the reliability of the most popular functionals were drawn by Baranowska-Ła ˛czkowska et al., based on the results of computations of excess properties for CO-and N 2 -(HF) n complexes with n = 1,. ..,4. 41 These authors reported that in the case of long chains, range-separated hybrid functionals or functionals with a large fraction of exact exchange are necessary to avoid errors exceeding 20%. 41Maroulis employed the B1LYP, B3LYP, B3PW91 and mPW1PW91 exchange-correlation functionals and determined the interaction-induced properties of the water dimer. 36Their comparison against the CCSD(T) reference values revealed that none of the functionals correctly predicts the sign (and magnitude) of the interaction-induced mean second hyperpolarizability. 36 Moreover, a recent study by Bulik et al. demonstrates that computations of vibrational contributions by means of density functional theory still pose serious challenges. 74o bring the discussion to a close, let us conclude that an insufficient amount of data regarding the performance of DFT was one of the factors that motivated us to employ two decomposition schemes for the excess property partitioning.We thus decompose electronic and vibrational contributions to electric excess properties, using the ab initio VP-EDS and ETS schemes with a selection of exchange-correlation functionals.
2.3.1 VP-EDS scheme.According to the VP-EDS scheme, the total interaction energy of a dimer, calculated by a supermolecular approach in the dimer-centered basis set (DCBS), 75 using the second-order Møller-Plesset perturbation theory (MP2) is partitioned into the Hartree-Fock (HF) and the electron correlation interaction energy components: The HF term can be further partitioned into the electrostatic interactions of unperturbed monomer charge densities, e (10) el , as well as the associated exchange repulsion (DE HL ex ), and the charge delocalization (DE HF del ) term encompassing the induction and the associated exchange effects due to the Pauli exclusion principle: The second-order electron correlation term, DE MP2 corr , DE MP2 corr = e (12) el,r + e (20) disp + DE (2) ex (17)   includes the second order dispersion interaction, e (20) disp , as well as the electron correlation correction to the first order electrostatic interaction, e (12) el,r , and the remaining electron correlation effects (DE (2) ex ).The latter term accounts mainly for the uncorrelated exchange-dispersion and electron correlation corrections to the Hartree-Fock exchange repulsion. 71,73The e (10) el and the e (20) disp terms are obtained in the standard polarization perturbation theory, whereas the e (12) el,r term is calculated using the formula proposed by Moszyn ´ski et al. 76 The indices in parentheses denote perturbation orders in the intermolecular interaction operator and intramonomer correlation operator, respectively.
2.3.2ETS scheme.The ETS approach, which is usually carried out in the framework of density functional theory, can be employed to analyze chemical bonds between the fragments of a molecule.The energy of chemical bonds (D e ) is defined as the difference between the total energy of a molecule E mol and the sum of the energies of the fragments E frag : The dissociation energy may be split into two terms, the preparation energy (DE prep ) and the interaction energy (DE int ): The first one is given by two contributions: (i) the energy required to distort the equilibrium geometry of the isolated fragments to the geometry that they have in the molecule; and (ii) the energy necessary to excite the electronic ground states of the fragments to the states with the electronic configuration suitable to form the bonds.Obviously, if the fragments are atoms, DE prep = 0.The interaction energy, which is the difference between the total energy of the molecule and the sum of the energies of the prepared fragments, can be split into three different contributions: (i) the electrostatic interaction energy between the frozen charge densities of the fragments (DE el ); (ii) the exchange repulsion energy due to the antisymmetrization and renormalization of the Hartree product of the fragment wavefunctions (DE Pauli ); and (iii) the orbital relaxation energy (DE orb ): It is interesting to note that if this scheme is applied to a weakly bound complex with the fragments chosen as interacting molecules, it becomes essentially equivalent to the partitioning of DE HF int in VP-EDS.However, if carried out within the KS-DFT framework, the obtained contributions also account for the intra-monomer correlation and to some extent the intermonomer correlation effects (i.e.dispersion interaction), which in VP-EDS are included via the Møller-Plesset perturbation theory.In this paper we compare the ab initio VP-EDS and KS-DFT ETS partition of the electronic and vibrational hyperpolarizabilities to analyze the performance of BLYP, LC-BLYP and LC-BLYP-dDsC functionals.

Software and computational details
Geometry optimizations and vibrational structure calculations using MP2, CCSD and CCSD(T) methods were performed using the Gaussian suite of programs 77 while property calculations were carried out using custom computer programs.VP-EDS calculations were carried out using a modified version of the Gamess (US) program. 78,79Finally, all DFT calculations were done using ADF program. 80

Results and discussion
In a series of papers, Baranowska and collaborators carried out the assessment of basis sets on the excess electric properties [24][25][26]38,41,42 and advocated the use of basis sets specifically tailored for such types of calculations (for the analysis of schemes for eliminating the basis set superposition error in calculations of the properties in question we refer to ref. 81 and 82). This is  line with the conclusions presented by Skwara et al. 22 However, as demonstrated recently by some of us, property-oriented basis sets are not always the best choice to predict the nuclear relaxation (hyper)polarizabilities. 83 On the other hand, it was reported in the same study that the average errors in these properties associated with the aug-cc-pVTZ basis set can be negligible (i.e. less than 1%) for small molecules.Considering these findings, in the present study we used Dunning's correlationconsistent basis sets aug-cc-pVXZ (X = D, T, Q) [84][85][86] for VP-EDS calculations.
The structure and geometrical parameters for the HCN monomer and the dimer are shown in Fig. 1.The most remarkable change in the monomer geometry upon the complex formation is the increase in the C-H bond length from 1.065 Å (isolated monomer) to 1.071 Å (monomer A in the complex AB).There is also an accompanying shortening of the C-N bond length by 0.002 Å in the case of monomer B. In comparison with the isolated monomer, upon the complex formation one thus expects a much more significant perturbation to the vibrational structure and properties of A than to B. For the exhaustive analysis of the energetics of (HCN) 2 isomer formation we refer to the study of Smith et al. 87 and to the ESI in ref. 27 and references cited therein.
MP2 and coupled-cluster electronic and vibrational contributions to a and b of the HCN monomer and the dimer are shown in Table 1.The electronic electric-dipole properties of the monomer and the dimer were also studied by other authors and there is a good agreement between their results and the values reported in Table 1. 26,27,88It should be highlighted that there is a satisfactory convergence of monomer and dimer properties with respect to the basis set size, and the aug-cc-pVTZ basis set already delivers accurate estimates.Then, although in the following analysis we will focus on the data obtained at the MP2/aug-cc-pVQZ level, the same conclusions are obtained from aug-cc-pVDZ and aug-cc-pVTZ basis set results.We can also see that the MP2 results are in fairly good agreement with those provided by more sophisticated CCSD and CCSD(T) methods.This is in line with results obtained for the polarizabilities of long polyynes. 89,90 comparison of the properties of the monomer with those of the dimer reveals a few interesting conclusions.For both systems the electronic polarizability is much larger than the nuclear relaxation polarizability, with the a nr /a e ratio increasing from 0.008 (monomer) up to 0.097 (dimer).In contrast to polarizability, b nr prevails over b e for both the monomer and the dimer, as the b nr /b e ratio is about 15 for the monomer and 5 for the dimer.Upon the complex formation there is a large increase in the nuclear relaxation contributions to polarizability (roughly 30-fold) and first hyperpolarizability (roughly 20-fold).In order to elucidate the origins of this enhancement we have performed normal mode analysis based on the Bishop-Kirtman perturbation theory.Fig. 2 shows that the major contribution of double harmonic diagonal terms (i.e.[m 2 ] (0,0) and [ma] (0,0) ) is due to the vibrational normal mode of frequency n 3 = 123 cm À1 , which corresponds to the intermolecular symmetric stretching along the principal symmetry axis (see also Fig. 3).
There is only a negligible net anharmonic contribution to nuclear hyperpolarizability of the monomer (3.28%).On the contrary, although the electrical and mechanical anharmonic contributions to b nr of the dimer cancel each other out to a Fig. 1 Geometrical parameters (given in [Å] units) for the monomer and the dimer of HCN and the orientation of molecules in the Cartesian coordinate system.Equilibrium geometries were obtained at the MP2/ aug-cc-pVQZ level of theory.
Table 1 Electronic and nuclear relaxation polarizability and first hyperpolarizability of the HCN monomer and dimer (given in a.u.) at equilibrium geometries corresponding to the indicated methods  In contrast to other problematic hydrogen-bonded systems, like for instance the water dimer or the hydrogen fluoride dimer studied by Eckart and Sadlej, 91 the data in Table 1 do not indicate that the BKPT expansion is likely to diverge.
We will now turn to the excess properties.Table 2 contains the values of electronic and nuclear relaxation interactioninduced polarizability and hyperpolarizability computed based on eqn (11) and (12).Although we noted earlier that a e is far larger than a nr for both the monomer and the dimer, the excess electronic and nuclear relaxation polarizabilities are similar in magnitude and they amount to 4.63 and 5.10 a.u., respectively.On the contrary, the excess Db nr is much larger than its electronic counterpart and the Db nr /Db e ratio is nearly 7. Table 2 also contains the values of pseudo-NR contributions, i.e.P nr (A AB ) and P nr (B AB ).These quantities arise from the definition of interaction energy which involves geometries of monomers corresponding to a dimer.Table 2 shows that there are large differences between b nr (A AB ) and b nr (B AB ).The comparison of a nr and b nr values listed in Tables 1 and 2 reveals that P nr (A AB ) deviates more substantially from the values for the isolated monomer than P nr (B AB ).Both a nr (A AB ) and b nr (A AB ) have opposite signs compared to the corresponding properties of the monomer at its equilibrium geometry.Presumably, a major factor contributing to observed discrepancies is the formerly mentioned elongation of the C-H bond in monomer A upon complex formation.
In order to shed more light on the pseudo-NR contributions, we employed an approach proposed by Ingamells et al. who studied vibrational corrections to electric properties at nonequilibrium geometries. 92,93We shall focus solely on a nr as the extension of this approach to nuclear relaxation hyperpolarizability is not straightforward.Following their treatment, we define the geometry correction (a nr,gc zz ) to the nuclear relaxation polarizability as: where k runs over all vibrational normal modes Q and E is the total energy at non-optimum geometry.The largest normalmode contributions to a nr,gc zz are shown in Table 3.The major contributions to a nr,gc are due to two vibrational modes (k = 3 and k = 4).The latter vibrational mode contribution prevails in the case of monomer A and it corresponds to C-H stretching.It comes as no surprise that this vibrational mode does not make a significant contribution to geometry correction for monomer B. As pointed out by Ingamells et al. the geometry correction is approximately equal to the difference between the electronic   polarizability at the non-optimum (a e (R neq )) and optimum geometries (a e (R eq )): The results in Table 3 show that there is an excellent agreement between a nr,gc and the difference given on the r.h.s. of the above equation.Finally, we briefly comment on the fact that the values reported in Tables 2 and 3 were obtained, for the sake of consistent analysis, employing the monomer-centered basis sets (MCBSs) 75 in the case of pseudo-NR contributions.The differences in excess properties between MCBS and DCBS results at the MP2/aug-cc-pVQZ level are not substantial as Da e + Da nr is 9.73 a.u.for MCBS and 9.66 a.u.for DCBS, while Db e + Db nr is À307.71a.u.for MCBS and À301.58 a.u.for DCBS.
We shall now analyze the results of electronic and vibrational excess property decomposition in terms of interaction energy components.Table 4 shows the results obtained based on the VP-EDS scheme.The excess electronic (hyper)polarizabilities of HCN aggregates have been thoroughly studied by some of us. 27ote, however, that in the present study, electronic (hyper)polarizabilities and equilibrium geometries are determined at the very same level of theory.This is the factor contributing to small differences between the present values and those given in ref. 27.With this in mind, we only briefly comment on our results of calculations of electronic properties.Electronic excess polarizability is dominated by the electrostatic contribution which prevails over the remaining components.Exchange-repulsion and delocalization components are similar in magnitude but of opposite signs while the sum of electron correlation corrections is negligible.
Much more complicated interplay of interaction energy terms is found in the case of electronic contributions to excess first hyperpolarizability, i.e. one finds that the first-order electrostatic component is quenched by the associated exchange repulsion, and the delocalization component makes the major contribution to excess properties.Except a significant contribution of electron correlation correction to the electrostatic component, the remaining terms computed at the correlated level are not substantial.
We now turn to the decomposition of excess nuclear relaxation polarizability and hyperpolarizability.There are several interaction energy components which are large in magnitude and make a substantial contribution to Da nr .These include: first-order electrostatic (+194%), exchange repulsion (À246%) and delocalization (+143%) components.As indicated by their signs given in parentheses, the first two terms cancel each other out to a large extent.The positive Hartree-Fock component is further slightly increased by the sum of dispersion contribution (+40%) and the electron correlation correction to exchange energy (À25%).In the case of excess nuclear relaxation first hyperpolarizability the pattern of interactions is similar (i.e.first-order electrostatic (+175%), exchange repulsion (À265%), delocalization (+143%), dispersion (+40%), and electron correlation correction to exchange repulsion (À25%)).Interestingly, there is an identical interference pattern between the negative and positive interaction components for Da nr and Db nr .Indeed, even the relative magnitude of the different components of Da nr and their counterparts for Db nr are strikingly similar.
Finally, we will discuss the results of calculations of electronic and nuclear relaxation contributions to a and b performed at the DFT/aug-TZ2P level of theory (cf.Table 5).The aug-TZ2P basis set is composed of Slater-type orbitals (STOs), which have the correct nuclear cusp and long-range behavior, as opposed to Gaussiantype orbitals (GTOs).Therefore, the properties converge much quicker with STOs than with GTOs with respect to the basis set size. 94he calculations were performed employing three functionals, namely: BLYP, LC-BLYP and dispersion corrected LC-BLYP-dDsC. DP HL ex DP HF del DP (12) el,r DP (20) disp DP (2)  The relative error with respect to the CCSD(T)/aug-cc-pVTZ value of a e is 6.13%, 2.62% and 1.92%, for BLYP, LC-BLYP and LC-BLYP-dDsC, respectively.Note the excellent performance of the latter two functionals in predicting electronic polarizability.In the case of electronic first hyperpolarizability, the corresponding relative errors with respect to CCSD(T)/aug-cc-pVTZ are much larger at the DFT/aug-TZ2P level of theory and they amount to 37.50%, 25.36% and 27.50% for BLYP, LC-BLYP and LC-BLYP-dDsC, respectively. Alsor a nr , the relative errors with respect to the CCSD(T)/aug-cc-pVTZ value are larger than for a e , i.e. 33.33%, 10.62% and 6.58% for BLYP, LC-BLYP and LC-BLYP-dDsC, respectively.Similar discrepancies are found in the case of b nr , i.e. the errors are equal to 30.83% and 0.02% for BLYP and LC-BLYP. Due to lnumerical instabilities, the value of this property is not shown for the dispersion corrected LC-BLYP-dDsC functional.Therefore, for the HCN dimer, the LC-BLYP functional not only delivers much more accurate estimates in comparison with the BLYP for the electronic contributions of a and b, but also for their nuclear relaxation counterparts.Table 5 also contains the results of decomposition of electronic and nuclear relaxation interaction-induced properties based on the ETS scheme.As already highlighted in the previous section, when the chosen fragments are interacting molecules, the ETS scheme is essentially equivalent to the partitioning of DE HF in VP-EDS.However, the DFT-ETS terms also partly account for electron correlation effects.Therefore, the DFT-ETS electrostatic term should be compared with VP-EDS DP (10) el + DP (12) el,r ; DFT-ETS exchange repulsions with VP-EDS DP HL ex + DP (2) ex ; and DFT-ETS orbital relaxation with VP-EDS DP HF del .In what follows we will compare the VP-EDS data determined at the MP2/aug-cc-pVQZ (GTOs) level of theory basis set with the results obtained based on the DFT-ETS scheme and using the aug-TZ2P (STOs) basis set.The absolute values of the DFT-ETS electrostatic, exchange repulsion and orbital relaxation terms contributing to Da e , Da nr , Db e and Db nr are larger than their VP-EDS counterparts.
The only exceptions are the electrostatic DFT-ETS terms of Db e and the electrostatic and exchange terms of Da nr determined by employing the LC-BLYP-dDsC functional.The discrepancy is larger for the orbital relaxation term than for the other two terms, with the exception of LC-BLYP and LC-BLYPdDsC Db e .Since the DP (2) disp term was not explicitly taken into account in the above comparison between VP-EDS and DFT-ETS terms, at least for BLYP and LC-BLYP, the observed discrepancies could partly be attributed to the missing (or non-representative) dispersion effects in DFT calculations.Our results show that the ability of DFT methods to reproduce the ab initio results is directly related to the magnitude of the orbital relaxation terms.That is, the origin of the too large DFT NLOP properties is the overshooting of the orbital relaxation term.Finally, we focus on the effect of the dispersion correction of LC-BLYP-dDsC results.The dDsC corrections to Da, and Db contributions are far smaller than the DP (2) disp term, and do not systematically improve LC-BLYP values.Thus, although dDsC correction improves the accuracy of the interaction energy, the same is not true for the excess properties.

Summary and conclusions
The concept of excess properties is a convenient framework to gain insight into the effect of intermolecular interactions on the electric properties of molecular complexes and clusters, which can also be relevant to set up appropriate models for calculations of bulk properties.Moreover, in combination with an appropriate interaction energy decomposition scheme this concept enables us to understand the physical origins of the excess properties.This paper reports the results of a pioneering investigation of nuclear relaxation contributions to the excess electric properties of a model molecular complex.Due to the presence of the hydrogen bond, which plays a pivotal role in supramolecular chemistry, the HCN dimer was chosen as a model system in our study.In order to analyze the interaction-induced nuclear relaxation (hyper)polarizabilities, the approach for the evaluation of purely electronic excess properties has been extended to the total (i.e.electronic and vibrational) excess properties.The new scheme relies on the computationally efficient Bishop-Hasan-Kirtman method for the partitioning of the total value into the electronic and nuclear relaxation contributions.Since in our extension both contributions to electric properties can be expressed in terms of the interaction energy derivatives with respect to the external field, it is straightforward to use an appropriate interaction-energy decomposition scheme and express the excess properties as the sums of physically meaningful contributions.In our study, we apply the VP-EDS scheme at the MP2/aug-cc-pVQZ level and the ETS scheme by employing selected DFT functionals (BLYP, LC-BLYP, and LC-BLYP-dDsC) with the aug-TZ2P basis set.A significant increase of both electronic and vibrational electric properties upon the complex formation is reflected in the magnitude of the corresponding excess quantities.Moreover, let us highlight the particularly large value of the excess nuclear relaxation first hyperpolarizability, which is roughly 7 times larger than its electronic counterpart.From the VP-EDS analysis it follows that: (i) while the electronic excess polarizability is strongly dominated by the electrostatic contribution, in the case of the electronic excess first hyperpolarizability, the electrostatic component is practically quenched by the exchange repulsion contribution, and then it is the delocalization term which finally plays the major role; (ii) there are several interaction energy components which are large in magnitude and make a substantial contribution to nuclear relaxation excess polarizability.It was observed that the first-order electrostatic contribution cancels out with the negative exchange repulsion term to a large extent resulting in a positive value of Da nr due to the contribution from the delocalization and dispersion terms.In the case of excess nuclear relaxation first hyperpolarizability, the pattern of interaction contributions is very similar to that for Da nr , both in terms of their sign as well as relative magnitude.
For the functionals explored in this work all three interaction energy terms within the DFT-ETS scheme contribute to the overestimation of excess properties.Nevertheless, the largest overshooting is always given by the orbital relaxation term.Finally, one may conclude that the LC-BLYP and LC-BLYP-dDsC functionals, which have a smaller orbital relaxation term than BLYP, are much more successful in predicting excess properties.Thus, comparison between VP-EDS and DFT-ETS partition schemes can be applied as a useful tool to analyze in depth the performance of DFT functionals for computing NLO properties and obtain helpful hints to pursue in the future design of more accurate functionals.

Fig. 2
Fig.2Convergence of longitudinal components of [m 2 ] (0,0) and [ma] (0,0) for the HCN dimer with respect to the vibrational normal modes (dotted curves).Individual contributions are represented by gray bars.All results were obtained at the MP2/aug-cc-pVQZ level of theory.

Fig. 3
Fig. 3 Atomic displacements for vibrational normal mode of frequency n 3 = 123 cm À1 for the HCN dimer.See the text for more details.

Table 2
Electronic and nuclear relaxation contributions to the polarizability and first polarizability of the HCN monomers and the HCN dimer.All data in this table correspond to MCBS calculations and were obtained at the MP2 level This journal is © the Owner Societies 2016 Phys.Chem.Chem.Phys., 2016, 18, 22467--22477 | 22473

Table 3
The geometry correction (a nr,gc zz ) to the diagonal nuclear relaxation polarizability together with non-vanishing normal mode (k) contributions.All results are given in a.u. and were obtained at the MP2/aug-cc-pVQZ level

Table 4
Results of decomposition of electronic and nuclear relaxation polarizability and first hyperpolarizability (given in a.u.) of the HCN dimer, performed at the MP2 level of theory.Note that each term, including DP HF and DP MP2 , was independently computed by numerical differentiation (corresponding numerical errors are given in the ESI)

Table 5
Results of decomposition of electronic and nuclear relaxation polarizability and first hyperpolarizability (given in a.u.) of the HCN dimer, performed at the DFT/aug-TZ2P level of theory.Corresponding numerical errors are given in the ESI This journal is © the Owner Societies 2016 Phys.Chem.Chem.Phys., 2016, 18, 22467--22477 | 22475