Partitioning of interaction-induced nonlinear optical properties of molecular complexes. I. Hydrogen- bonded systems†

Understanding the effects of different fundamental intermolecular interactions on the nonlinear optical properties is crucial for proposing efficient strategies to obtain new materials with tailored properties. In this study, we computed the electronic and vibrational (hyper)polarizabilities of ten hydrogen-bonded molecular complexes employing the MP2, CCSD and CCSD(T) methods combined with the aug-cc-pVTZ basis set. The vibrational contributions to hyperpolarizabilities included nuclear-relaxation anharmonic corrections. The effect of intermolecular interactions was analyzed in terms of excess properties, which are defined as the difference between a property of the complex and the net properties of the noninteracting subsystems. Considering systems covering a wide range of the hydrogen bond strengths, the electronic and vibrational excess (hyper)polarizabilities were decomposed into different interaction energy contributions (electrostatic, exchange, induction and dispersion). This systematic study, the very first of this kind, revealed that the physical origin of the electronic and vibrational excess properties is completely different. In the case of vibrational contributions, the decomposition pattern is very similar for the polarizability and first and second hyperpolarizabilities. The exchange contributions to excess vibrational properties are the largest and they have different sign than the electrostatic, induction and dispersion terms. On the other hand, no general patterns can be established for the electronic excess properties.


Introduction
It is now well established that bottom-up engineering of materials for nonlinear optics (NLO) applications is an effective strategy. 1,2At the molecular level, linear and nonlinear optical properties are governed by the electric dipole polarizability (α) and first (β ) and second (γ) hyperpolarizabilities, respectively. 3nderstanding the factors that determine the magnitude of the (hyper)polarizabilities is essential to model new molecules and supramolecular complexes characterized by large nonlinear optical responses, required to design novel materials for photonic applications.5][6][7] These effects can be conveniently discussed in terms of excess properties, which are defined as the difference between a property of the complex and the net properties of the noninteracting subsystems.Although there has been substantial progress recently in studies of the excess properties and related collision-induced spectroscopy of atomic pairs [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] and molecular complexes,  the interplay of fundamental interaction types in the context of high-order electrical properties of molecular complexes has not received appropriate attention yet. One aproach to analyze NLO properties of interacting species is the partitioning of total interacting system (hyper)polarizabilities into terms of different physical origins.In a pioneering work in 1992, Fowler and Sadlej analyzed the different contributions to the collisioninduced linear polarizabilities of He• • • F − and He• • • Cl − within the framework of the theory of intermolecular interactions.9 They found that the exchange-overlap interaction terms were indispensable for the quantitative description of the linear polarizability.Four years later Bishop and Dupuis analyzed the different terms of the linear polarizability and second hyperpolar-izability of He dimer using perturbation theory.46 In the same year, Heijmen et al. formulated a symmetry-adapted perturbation theory for interaction-induced polarizabilities of weakly bound complexes and reported expressions for the contributions up to the second order in the intermolecular potential based on the finite-field numerical differentiation of the corresponding energy terms.12 Some of the present authors have applied this method to small molecular dimers 25,42 and trimers, 29 as well as to 70 Watson-Crick B-DNA pairs, 41 analyzing the role of the electrostatic, exchange-repulsion, exchange-induction, dispersion and electron correlation terms of the purely electronic interactioninduced properties.It was observed that in many cases the electronic excess polarizabilities showed similar patterns, i.e., they were predominantly determined by the first-order electrostatic and the second-order induction (polarization) terms which were diminished by the exchange-repulsion effects (cf.HF dimer 25 , quasi-linear dimers of HCN, 29 urea, diformamide, 4-pyridone, 4nitroaniline, the complex of hydrogen fluoride with nitroacetylene 42 and stacked cytosine dimers 41 ).In the case of Watson-Crick guanine-cytosine and adenine-thymine pairs the exchange and induction terms also played dominant roles prevailing over the electrostatic contribution.41 On the other hand, a regular pattern of interaction types was not found for the excess first hyperpolarizability of hydrogen-bonded species, since two systems having very similar nature of interactions were often characterized by an entirely different origins of response properties.42 To avoid confusion it is worth pointing out that in the above discussion we refer to various contributions to response properties.For instance, by the first-order electrostatic term we understand the contribution to excess (hyper)polarizabilities due to the change of electrostatic interactions in a complex caused by the presence of an uniform external electric field.The details of our approach are summarized below.
][49][50][51][52][53][54][55][56] The studies summarized above share a common feature, which is the neglect of the vibrational contributions to molecular (hyper)polarizabilities. [57][58][59][60][61] The only exception is a pioneer study on the vibrational contribution of excess nonlinear properties of HCN dimer. 62As demonstrated by several studies, molecular vibrations play a key role in many nonlinear optical processes, especially in those involving static fields.There are no studies of excess vibrational hyperpolarizabilities analyzing their origins, although some authors indicated that these contributions may prevail over electronic ones. 63The present work, aiming at filling this gap, is the first systematic study revealing the interplay of various interaction types in excess electronic and vibrational (hyper)polarizabilities of weakly bound systems.As such, it contributes to the fundamental understanding of the physical origins of electronic and vibrational contribution to excess (non)linear optical properties.
As a case study we selected ten linear hydrogen-bonded complexes of unsaturated compounds.The choice of unsaturated hydrogen-bonded species is natural in view of the preliminary studies of such systems summarized above.Limiting this project to linear complexes allows us to compute highly accurate welldefined interaction components of their electronic and nuclearrelaxation (hyper)polarizabilities, which could serve as benchmarks in further studies.Let us underline that in all cases the linear configuration represents equilibrium geometry of the system.

Theory
The properties of interest in this study are the electronic P e and nuclear-relaxation P nr contributions to a property P (P=α, β , γ).The total property P is the sum of these two terms and the curvature contribution P curv : 58 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 change of the equilibrium geometry and the shape of the potential energy surface, respectively.The nuclear relaxation contributions to β and γ contain all the lowest-order anharmonic corrections.The curvature contributions, which include the rest of the high-order anharmonicity, are far more computationally expensive and hence they are not computed here.A efficient method to compute the nuclear-relaxation hyperpolarizabilities is that proposed by Bishop, Hasan and Kirtman (BHK), hereafter denoted as FF-NR. 64,65In the present study, we confine the analysis solely to static properties, including polarizability and first and second hyperpolarizabilities.As far as the nuclear relaxation contribution is concerned, following BHK we define: 64 where µ i (F, R F ) is the i-th component of the dipole moment obtained at the relaxed equilibrium geometry in the presence of an external electric field F and µ i (0, R 0 ) is the same property for field-free conditions.The expansion coefficients yield the static properties: When applying these FF-NR formulas the geometry relaxation must not include the rotations of the molecule through the alignment of the permanent or induced dipole moment in the field direction (indeed this may be the easiest way for the system to decrease its energy).Therefore, the field-dependent optimization are performed strictly maintaining the Eckart conditions. 66Such optimizations can efficiently be carried out with the aid of a procedure developed by Luis et al. 67 In our recent study, we have shown how to partition the nuclear-relaxation (hyper)polarizabilities into contributions due to the various interactions types. 62In the remainder of this section, we will briefly outline our approach and the adopted intermolecular-interaction energy decomposition scheme.
The equilibrium geometry of a complex composed of A and B in an applied electric field, F, is denoted by AB F , while E(F, AB F ), E(F, A AB,F ) and E(F, B AB,F ) are the field-dependent energies of AB, A and B, respectively, at the geometries corresponding to the field-induced relaxed structure of AB F .AB stands for the fieldfree equilibrium geometry of the complex.The field-dependent interaction energy at field-relaxed geometry AB F is thus given by: and then the nth derivative of ∆E int (F, AB F ) with respect to electric field components along Cartesian directions i, j, . . . is given by: It can be easily shown that the nth-order derivative of an interaction energy component, E int,X , can be expressed as the sum of the interaction-induced electronic and pseudo-nuclear relaxation contributions to a corresponding property P: ∆P nr of the previous equation does not correspond to the rigorous expression of the nuclear relaxation contribution to P nr because for its definition we use E(F, A AB,F ) and E(F, B AB,F ) instead of E(F, A F ) and E(F, B F ), as would be appropriate for P nr .Finally, the above equation can be presented in a form which allows to analyze the ∆P nr contribution in terms of interaction-energy components ∆E int,X : 7][78][79] In this scheme the interpretation of individual components refers to the intermolecular perturbation theory. 807][78][79] According to VP-EDS, the total interaction energy of a dimer, calculated by a supermolecular approach in the dimer-centered-basis set (DCBS), 81 at 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, ε el , as well as the exchange repulsion (∆E HL ex ), and the charge delocalization (∆E HF del ) terms encompassing the exchange effects due to the Pauli antisymmetry principle and the induction, respectively.
The second-order electron correlation term, ∆E MP2 corr : includes the second order dispersion interaction, ε disp , as well as the electron correlation correction to the first order electrostatic interaction, ε (12)   el,r , and the remaining second order electron correlation effects (∆E ex ).The latter term accounts mainly for the uncorrelated exchange-dispersion and electron correlation corrections to the Hartree-Fock exchange repulsion. 78,80The ε (10)  el and the ε (20)   disp terms are obtained using standard polarization perturbation theory, whereas the ε el,r term is calculated using the formula proposed by Moszyński et al. 82 The indices in parentheses denote perturbation orders in the intermolecular interaction operator and intramonomer correlation operator, respectively.

Software and Computational Details
Geometry optimizations and vibrational structure calculations using MP2, CCSD and CCSD(T) methods were performed with the GAUSSIAN 83 and ACES II 84 of programs, whereas property calculations were carried out employing custom computer programs.VP-EDS calculations were carried out using a modified version of the GAMESS (US) program. 85,86aranowska et al. performed the assessment of basis sets on the electronic excess properties [26][27][28]40,43,44 and suggested that the use of basis sets specifically tailored for such type of calculations (for the analysis of schemes for eliminating the basis set superposition error in calculations of the properties in question we refer to Refs 30,87,88). However,as demonstrated recently by some of the present authors, property-oriented basis sets are not always the best choice to predict the nuclear relaxation (hyper)polarizabilities. 89 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 were negligible (i.e. less than 1% compared to aug-cc-pVQZ results) for small molecules.Therefore, in this study we used Dunning's correlation-consistent augcc-pVTZ basis set [90][91][92] in all property calculations.

Results and Discussion
In order to systematically analyze the role of particular interaction types and possible common patterns in interactioninduced properties of selected hydrogen bonded systems, the VP-EDS approach was employed to partition the electronic and nuclear-relaxation (hyper)polarizabilities of 10 linear molecular complexes: HNC Due to the presence of a single lone pair on the hydrogen bond acceptor, the optimized geometries of all studied complexes belong to the C ∞v symmetry point group and were aligned along the Cartesian z-axis.Only the longitudinal components of all properties were calculated, i.e. µ z , α zz , β zzz and γ zzzz .As indicated in the preceding sections, the decomposi- tion of interaction-induced properties based on VP-EDS scheme is equivalent to the partitioning of excess properties evaluated based on the MP2 method.The rationale behind this choice is that interaction-energy decomposition, an essential ingredient of the approach followed here, is not yet feasible for all studied systems at higher levels of theory.However, in order to assess the performance of MP2 method (and CCSD as well), the electronic and nuclear-relaxation (hyper)polarizabilities of all 10 complexes were computed using the CCSD(T)/aug-cc-pVTZ approach.The summary of results is shown in Fig. 1.By and large, the CCSD method does not bring any systematical improvement upon the MP2 method as far as the average unsigned errors are concerned.In fact, except for β e and β nr , the errors associated with MP2 are smaller than those corresponding to CCSD.To sum up, our results indicate that the MP2 method used for interaction-induced property decomposition offers satisfactory accuracy in predicting total properties of the studied complexes, and the corresponding average (maximum) percentage errors do not exceed 12% (30%) and in most cases are much smaller.
Before we analyze the interplay of interaction types governing the excess properties, let us first discuss the results of interaction energy decomposition for the studied systems shown in Figure 2. The systems are ordered according to decreasing value of the total MP2 interaction energy, computed at the respective MP2/augcc-pVTZ equilibrium geometries.This partitioning shows that, in terms of absolute values, the exchange repulsion and electrostatic terms prevail over other interaction energy components.However, these two first-order perturbation theory terms cancel each other out to a large extent and the first order interaction energy is rather small.Note that the delocalization and dispersion components are always attractive (negative) and make a contribution to the overall stabilization much larger than the net first order interaction.Although the total interaction energy shows substantial variation on passing from HNC• • • HNC (the largest value) to N 2 • • • HF (the smallest value), the overall picture of relative contributions to the interaction energy is similar for the whole stud-Fig.2 Interaction energy partitioning for the studied species at the equilibrium geometries optimized using MP2/aug-cc-pVTZ method.The decomposition of interaction-induced electronic dipole moment (the corresponding nuclear relaxation contribution is by definition always zero) shown in Figure 3 reveals that this property is strongly dominated by the electrostatic contribution.This is an expected result observed in all previous studies.External electric fields induce electric multipole moments on interacting species which interact primarily via Coulomb forces even in systems such as helium dimer, whereas higher order effects are less important.The external field changes the electronic density distribution and influences also the induction and exchange-induction interactions described by the delocalization term, which also contributes noticeably to the excess dipole moment.The exchange repulsion contribution is in this case far smaller than its electrostatic counterpart.Indeed, for half of the molecular complexes studied, its contribution to the total excess dipole moment is practically negligible.Other interactions types do not play any qualitatively or quantitatively important role.
A similar pattern is found for the electronic excess polarizability (upper panel in Figure 4).The electrostatic contribution is again the largest one, but the exchange term is more important than for ∆µ el .Similarly, the delocalization contribution is also more pronounced than in the case of the excess dipole moment.Finally, it should not be overlooked that the Heitler-London exchange and the accompanying correlation correction have different sign from the remaining components, thus diminishing the overall interaction-induced polarizability.This holds for all the studied systems.The partitioning of the excess vibrational polarizability is shown in the lower panel of Figure 4. Similarly to what has been found for ∆α el , the exchange contributions have the opposite sign compared to the remaining components.However, there is a clear difference between electronic and vibrational excess polarizabilities.For the latter, the electrostatic contribution is in this case smaller than the corresponding exchange-repulsion for all studied systems (cf.Table S1 in ESI).This makes the firstorder contribution much less important and ∆α nr is to a large extent determined by the magnitude of the delocalization contribution.
Figure 5 shows the comparison of electronic (upper panel) and nuclear relaxation (lower panel) excess first hyperpolarizabilities.Note that due to numerical instabilities, the decomposition for HCN• • • HCCH complex was considered unreliable and it is not shown in Fig. 5.In the case of the electronic contribution, similarly as in the previously studied complexes, 42 it is difficult to point out any general trends holding for the whole set.Nevertheless, we note that with one exception (HCN• • • HCl), the Heitler-London exchange contribution has the opposite sign compared to the electrostatic and delocalization contributions.For several systems, e.g.HNC• • • HNC, HNC• • • HCN, FCN• • • HCN, the latter two terms make a similar contribution to excess electronic first hyperpolarizability.The overall contribution to ∆β el arising from the dispersion term is very small and its sign depends on the particular complex.Interestingly, for ∆β nr , the interplay of interaction types is much more regular and very similar to that observed for ∆α nr .Firstly the relative weight of the VP-EDS contributions are similar for the whole set (with the largest deviation for HCN• • • HCl; see Table S2 in ESI) and secondly the exchange contributions have opposite sign to the remaining components.Finally, the electrostatic contribution is less important than the exchange and delocalization terms.Based on the comparison of electronic and nuclear-relaxation terms we may conclude that for excess first hyperpolarizability, similarly to excess polarizability, the vibrational contributions are larger than the electronic counterparts.However, the differences in these two contributions to ∆β are much more pronounced than in the case of excess polarizability.
Let us now turn to the excess second hyperpolarizability (see Figure 6).This property involves fourth-order derivatives of interaction energy components and due to the instabilities of numerical differentiation, the results are shown only for five out of ten molecular complexes.However, even for this smaller subset it is possible to make very interesting observations.In the case of electronic excess second hyperpolarizability, the pattern regarding contributing terms is quite systematic.The signs of all terms remain constant for all studied complexes.The firstorder exchange-repulsion contributions are comparable in absolute magnitudes to their electrostatic counterparts and generally slightly larger, thus making the total first-order contribution less important than the higher-order terms.However, although the second-order electron-correlation terms ∆γ (12)   el,r and ∆γ (2)   ex have much greater relative weight, they also cancel each other to a large extent, similarly to the first-order terms.Then, even though the delocalization and dispersion terms appear to be much smaller, they importance is clear considering cancellation of correlated and uncorrelated electrostatic and exchange-repulsion terms.
Interestingly, in the case of ∆γ nr , one finds a very similar pattern to that already found for ∆α nr and ∆β nr .Again, the exchange contributions have opposite sign to the remaining components and the electrostatic contribution is less important than the exchange and delocalization terms.Similarly to lower order excess properties, the vibrational counterpart of ∆γ is much larger than the electronic one.Finally, for HCN• • • HCl complex, the relative dispersion contribution to ∆γ nr is unusually large comparing to other systems.
The HCN• • • HCl complex is interesting in its own right and deserves a deeper analysis, as the largest vibrational contributions to excess polarizability and first and second hyperpolarizabilities are found for this complex.We have already mentioned that the relative VP-EDS contributions to ∆α nr and ∆β nr do not fully follow the patterns observed for other complexes (Table S2 in ESI).For instance, the total correlation contribution for both quantities is about 30%, whereas it is typically about 15% for the other cases.Electrostatic, exchange, and delocalization contributions are also significantly larger compared to other dimers.In order to shed light on whether this complex exhibits any other unusual features (e.g.exceptionally large anharmonicity contribution), we provide a comparative analysis between HCN• • • HCl and HCN• • • HNC complexes.In the case of HCN• • • HCl we find that β nr zzz amounts to -760.49a.u.while ∆β nr zzz is nearly twice as large (-1217.59a.u.), showing a significant effect of intermolecular interactions, whereas for HCN• • • HNC, the corresponding values are -345.17a.u. and -635.79 a.u.Comparison of these two complexes shows that there is a two-fold increase in both proper- These modes correspond to intermolecular stretching vibrations (ν 3 , ν 5 ) and to intramolecular H-Cl (ν 9 ) and H-N (ν 13 ) stretching vibrations.To sum up, considering the anharmonicity the nature of ∆β nr zzz in the HCN• • • HCl complex does not differ significantly from other hydrogen bonded complexes.Its large magnitude can be rationalized in the harmonic approximation and attributed to lower harmonic frequencies and larger value of ( Finally, let us briefly comment on the overall effect of intermolecular interactions on vibrational properties.To this end, we will take a look at ∆P nr /P nr ratio.It turns out that this ratio increases with the order of property, i.e.: ∆α nr /α nr < ∆β nr /β nr < ∆γ nr /γ nr and such trend holds for vibrational contributions of all studied complexes.For example, for HCN• • • HCl it is: ∆α nr /α nr = 1.2 < ∆β nr /β nr = 1.6 < ∆γ nr /γ nr = 2.3 In the case of the very same complex, the corresponding electronic terms are: ∆α el /α el = 0.1 < ∆β el /β el = 1.4 > ∆γ el /γ el = 0.5 This comparison shows that the effect of intermolecular interactions is significant and, by and large, the ∆P el /P el is much smaller than the corresponding ratio for vibrational counterpart, as demonstrated in the case of HCN• • • HCl complex.

Conclusions
We studied the effects of different fundamental intermolecular interactions on nonlinear electrooptic properties of ten hydrogenbonded molecular complexes: For all these weaklybound systems we computed electronic and vibrational (hyper)polarizabilities employing the MP2, CCSD and CCSD(T) methods combined with the aug-cc-pVTZ basis set.As far as the vibrational counterpart is concerned, the nuclear relaxation terms were determined, which include all the different types of lower-order anharmonic corrections.The electronic and vibrational excess properties, defined as the difference between properties of the complex and its noninteracting constituents, were decomposed into contributions with different physical origin (electrostatic, exchange, induction and dispersion).The partition in question was performed at the MP2/aug-cc-pVTZ level and its reliability was confirmed by CCSD(T) calculations.
Our results clearly demonstrate qualitative differences between vibrational and electronic excess properties.While the origins of nuclear relaxation contributions are consistent in all the studied systems, regardless of the analyzed property, their electronic counterparts do not show any apparent patterns.This may indicate that systems with similar structural motifs and nature of interactions should also have similar excess vibrational (hyper)polarizabilities.Generally, the exchange-repulsion contribution to the latter quenches the corresponding electrostatic, induction and dispersion terms.These conclusions were drawn for a series of hydrogen-bonded complexes.It is thus highly compelling to pursue, if the same holds for other structural motifs, e.g. for complexes in stacked alignments.This investigation is currently in progress.

Fig. 1
Fig. 1 Percentage errors (minimum, maximum, average) in properties w.r.t.CCSD(T) reference values for the studied set of 10 dimers.Average errors were computed based on absolute values.The numerical data can be found in ESI.

Fig. 3
Fig. 3 Breakdown of interaction-induced electronic dipole moment into various interaction types.

Fig. 4
Fig.4 Breakdown of electronic and nuclear-relaxation contributions to interaction-induced polarizability into interaction types.

Fig. 5
Fig. 5 Breakdown of electronic and nuclear-relaxation contributions to interaction-induced first hyperpolarizability into interaction types.

Fig. 6
Fig.6 Breakdown of electronic and nuclear-relaxation contributions to interaction-induced second hyperpolarizability into interaction types.