Partition of Optical Properties Into Orbital Contributions

Nonlinear optical properties (NLOPs) play a major role in photonics, electro-optics and optoelectronics, and other fields of modern optics. The design of new NLO molecules and materials has benefited from the development of computational tools to analyze the relationship between the electronic structure of molecules and its optical response. In this paper, we present a new means to analyze the response property through the partition of NLOPs in terms of orbital contributions (PNOC). This tool can be used to obtain a local representation of the NLOPs, providing a powerful visualization aid to connect the magnitude of the optical property with some parts of the molecule. Unlike other methods to analyze NLOPs, the PNOC decomposes the optical property into orbitals of the unperturbed system, furnishing this method with the ability to assess the performance of single- and multi-determinant electronic structure methods. PNOC can be also used to design small basis sets for an accurate description of large systems, saving a substantial amount of computer time for the calculation of optical properties.


Introduction
Nonlinear optical properties (NLOPs) have been extensively studied since the discover of the second-harmonic generation exhibited by a quartz sample irradiated by a ruby laser. 1 This discovery led to applications in optical communications, photonics, and optoelectronics, such as X-ray generation, the study of correlated photon pairs or biological imaging, spurring the quest for molecules and materials with high nonlinear optical (NLO) response. 2e design of new NLO materials and molecules has benefited from the evolution of quantum chemistry, which has provided the framework to compute NLOPs.In addition, the development of computational tools to analyze the relationship between the electronic structure of molecules and its optical response has greatly contributed to guide the design of nonlinear materials.Perhaps the simplest of such tools is the well-known two-level approach of Chemla and Oudar, 3 which simplifies the expression of the calculation of the hyperpolarizability, unveiling the most important molecular features that are responsible for the enhancement of the second-order nonlinear response.In the same vein, the sumover-states (SOS) technique, 4 which provides a perturbation theory approach to compute NLOPs, is currently still used as a way to assess the role of some electronic states in the magnitude of the NLOP.Additional works also considered the removal of states 5 or orbitals 6 in the SOS expression to estimate their importance.
One of the landmark papers on the analysis of NLOPs is due to Chopra and co-workers, 7 who studied orbital contributions to linear and NLOPs.Although the contributions sum up to the total NLOP value, they depend on the position of the origin of coordinates axis.Despite this unsettling downside, origin-dependent methods are widespread and the coordinate axis is usually fixed at the center of mass.Nakano et al. [8][9][10] defined the contributions to the n-order optical property as the n-order derivative of the density with respect to the electric field.Obviously, the sum of these quantities does not add up to the total value of the NLOP, but they are not origin dependent.The analysis has been extended to open-shells systems 11,12 using the density of effectively unpaired electrons, 13 which gives rise to pairwise contributions to the optical property.More recently, Geskin and Brédas have considered NLOP contributions based on the derivatives of the Mulliken charges with respect to the electric field, 14,15 providing origin-dependent contributions that do not add to the total values for all optical properties.The latter approach was originally proposed to study the static optical response and it was later extended to dynamic properties within the frameworks of time-dependent Hartree-Fock 16 and Kohn-Sham 17 levels.Localized orbitals have been also employed to decompose NLOPs. 17,18Other analyses provide a pairwise decomposition of the optical property.This is the case of the Hieringer and Baerends decomposition of the first hyperpolarizability within the framework of the Kohn-Sham density functional theory 19 and the recent work of Mandado and co-workers using field-induced orbitals. 20][23][24] Finally, one should also mention the decomposition of NLOPs in the real space.Most of these analyses [25][26][27][28][29][30] are framed in the quantum theory of atoms in molecule (QTAIM), 31 although Hirshfeld-based atomic partitions have been also used to analyze the polarizability. 32][39] Decomposing NLOPs into orbital contributions would help identifying pitfalls of electronic structure methods such as the wrong selection of active orbital spaces or the adequate size of orbital basis sets.
In this paper, we present the partition of NLOPs into orbital contributions (PNOC), a new means to analyze the response property and estimate its local contributions.PNOC provides a convenient way to relate linear and nonlinear optical properties to the structure of the molecule, adding a new tool for the construction of molecules with potential applications in photonics and optoelectronics.In addition, PNOC provides a feature that is absent in many other tools to analyze NLOPs: a decomposition in terms of the orbitals of the unperturbed system that can be also used to assess the performance of electronic structure methods and design conveniently small basis sets to study large molecules.
In the following we provide some background knowledge needed to introduce the PNOC and, afterwards, we give some examples that illustrate the potential of this tool.

Methodology
In the quantum-mechanical picture, under the perturbation of an electric field F F F, the Hamiltonian of a molecular system Ĥ is expressed as a sum of the unperturbed Hamiltonian Ĥ and the term that describes the interaction of the electrons (a) and nuclei (A) with the external electric field F F F: In the present study, the electronic nonrelativistic Hamiltonian and electronic part of the wavefunction are considered, with static and homogeneous electric fields taken into consideration.This prevents the contribution of higher multipole moments (like quadrupole or octupole) to the total energy. 40,43However, the present scheme could be easily generalized for such multipole moments.
The classic definition of the molecular static electric (hyper)polarizabilities comes from the Taylor expansion of the fielddependent dipole moment µ µ µ( ( (F F F) ) ): 7,40,41,43 where i, j, k and l can be any of the Cartesian components: x, y or z.In the above equation, the permanent dipole moment (µ i (0 0 0)) is altered by the linear response, quantified by the polarizability (α i j ), and the nonlinear response, represented by the first hyperpolarizability β i jk and second hyperpolarizability γ i jkl , which are consecutive partial derivatives of µ i (F F F): An equivalent description of the NLOPs can be obtained in the terms of the derivatives of the one-electron densities with respect to the external field utilizing the definition of the dipole moment µ i (F F F) in terms of the field-dependent one-electron density function ρ(r r r, F F F): [7][8][9] µ i (0 0 0) = − r i ρ(r r r, 0 0 0)dr r r Although the aforementioned density derivatives ρ ( j) (r r r), ρ ( jk) (r r r) and ρ ( jkl) (r r r) have been named after their NLOPs predecessors: α-, β -, γ-densities, respectively, 8,9 they are not rigorous property densities since their integration does not give the corresponding NLOP.
In a given basis (for instance atomic orbital (AO) basis), the induced dipole moment µ i (F F F) can also be defined in terms of the field-dependent one-particle reduced density matrix (1-RDM), D(F F F), and one-electron transition dipole moment matrix h (i) (with elements defined as h Further differentiation with respect to the external electric field leads to similar expressions for linear and NLOPs (note that the basis is unchanged for the perturbed and unperturbed systems, i.e., ∂ n h In the latter equations, D µν , D µν and D ( jkl) µν are the derivatives of 1-RDMs with respect to the external field (expressed in the same basis), which employ the equivalent notation to the one from Equations 5-7.

Partition of NLOPs into Orbital Contributions
In this section we describe the partition of NLOPs into orbital contributions (PNOC), which is based on the 1-RDM representation of the NLOPs.Only the PNOC expressions for the second hyperpolarizability component γ i jkl are given.An equivalent derivation and partition applied to α and β is available in the ESI.† Firstly, the electric field derivative of 1-RDM, D ( jkl) , is computed.In the current PNOC implementation, D ( jkl) are calculated with the finite field procedure using the AOs basis, because it is invariant to the applied electric perturbation.
Secondly, the h (i) and D ( jkl) matrices in the atomic orbital (AO) basis are projected to matrices M (i) and ∆ ∆ ∆ ( jkl) that are expressed in the basis of the natural orbitals (NOs) of the unperturbed (i.e., field-free) molecule * , giving rise to an equivalent representation of γ i jkl in the new basis, with the new matrices ) † , obtained through the transformation between two bases pq = χ p | r i | χ q is the i-th component of the transition dipole moment vector between the p-th and q-th natural orbitals.
where χ χ χ NO and φ φ φ AO are row vectors, and the LCAO coefficients are organized in columns of matrix C.
Finally, the selected property is partitioned into p NOs components, γ i jkl,p .In PNOC, each γ i jkl,p component is obtained under the assumption that the pair contribution of two NOs p and q is equally distributed between them: Although PNOC can use any arbitrary basis of orbitals, we employ NOs of the unpertubed system for the sake of convenience, as this orbital basis is often employed in the description of the 1-RDM of many-body wavefunctions.In the case of spin unrestricted computations, such as UHF (or broken spin-symmetry KS-DFT) we employ the unrestricted NOs defined by Pulay and Hamilton. 46ndeed, the PNOC can be applied to any quantum chemistry method provided a 1-RDM is available and, unlike other partition methods, it yields an exact decomposition of the NLOP.The main advantages of PNOC are the simplicity and consistency of the partition for the high-order NLOPs.
One of the drawbacks of most of the NLOPs decomposition schemes found in the literature is the so-called origindependency, which originates from the explicit dependency of dipole moment operator on the position vector. 8,14,15,17,18,25,47hile for the neutral system (total charge Q=0) with fixed orientation the total values of the properties (such as α zz , β zzz and γ zzzz ) do not depend on the relative position with respect to the center of the Cartesian system, the partitioned values in terms of 3D fragments or functions that depend on spatial coordinates usually do.Nevertheless, it is interesting to note that for systems that have the σ h symmetry plane along the target direction or systems that are centrosymmetric, the PNOC decomposition is free of the origin-dependency problem (see Section S.II in ESI † for further details).As previously done by many authors, we adopt here the convention to center the origin at the center of mass of the molecule. 17

Computational details and studied systems
In this work we analyze molecular (hyper)polarizabilities of the four different chemical models presented in Figure 1.In all calculations, the center of mass of the system was placed at the origin of the Cartesian coordinate system, and the system was rotated so that the main component of the diagonalized inertia tensor coincides with the z-axis.In this arrangement, only the longitudinal components of α zz and γ zzzz were studied.Unless stated otherwise, the highest possible symmetry constraints were applied.The electronic structure calculations were performed with the Gaussian09 48 computational package and the PNOC scheme was implemented within an in-house code. 49The Chemcraft program was used to generate all the orbital pictures. 50Our SOS computations include the first 20 ( 16) FCI (TDHF) Σ + u excited states (D 2h symmetry constraints were applied to the wavefunction).FCI computations were carried out with Molcas8, 51,52 whereas TDHF were done with Gaussian09.In all cases, excepting the (H 2 ) 3 chain (for which we have considered the geometry shown in Figure 1), we have studied the ground state geometry of the molecule.Benzene and p-benzyne were optimized using (U)B3LYP/aug-cc-pVDZ (the same basis set was employed for NLOP calculations), whereas for all-trans-hexatriene we employed the CAM-B3LYP/cc-pVDZ method.For the (H 2 ) 3 computations we employed the cc-pVDZ basis set of Dunning et al.. 53 This basis was chosen because its size is small enough to permit the FCI calculations of (H 2 ) 3 (6 electrons in 30 active orbitals), which is the only electronic structure method for which the SOS expressions yield the exact optical response.
To achieve the best stability in the 1-RDMs numerical differentiation (see Section S.III in ESI †), tight convergence criteria were used in the electronic structure computations (given in atomic units): 10 −11 in HF, 10 −11 in CASSCF total energies, 10 −10 and 10 −8 for the energy and the norm of the amplitudes vector in CCSD.Unrestricted coupled Hartree-Fock (UCHF) calculations were obtained from the data computed with Gaussian09 using the expressions given in section S.IV of the ESI †.

Comparison between SOS and PNOC schemelinear (H 2 ) 3 chain
In this section, we demonstrate the potential of the PNOC decomposition by comparing the NOs contributions to the ones provided by SOS (for a description of the SOS see Section S.IV of ESI †).To this end, we analyze α zz of a model system of three colinear hydrogen molecules, (H 2 ) 3 (see Figure 1).As expected, the SOS and the PNOC FCI values of the NLOPs are almost equal.A small negligible difference of 0.12 a.u.for α zz originates from the incompleteness of the set of FCI excited states used in the SOS formulation (see Table 1).For this system, α SOS zz using TDHF excited states is also in excellent agreement with the coupled HF value of α zz .Such good agreement between TDHF SOS and coupled NLOPs is in general not expected for larger molecular systems.On the other hand, UCHF α SOS zz is underestimated by 34.5% with respect to the coupled HF value.
Results of the SOS orbital partitioning and the PNOC scheme are compared in Table 1.In general, a semiquantitative agreement between SOS and PNOC is observed.The orbital contributions to the NLOP provided by the SOS and PNOC schemes are both in agreement with chemical intuition: the frontier orbitals, 2σ + g (HOMO) and 2σ + u (LUMO) in (H 2 ) 3 are the most important to the linear response (see Figure S.1 in the ESI † for a graphical representation of all NOs along with their occupancies).The relative contributions of the NOs to SOS and PNOC coupled α zz are very similar (differences smaller than 11 %) -especially for the HF wavefunction, for which exactly the same set of MOs of the ground state is used in both decompositions (differences smaller than 6 %).The differences between the SOS and PNOC contributions to FCI α zz originate from the approximate SOS orbital partitioning, which assumes that the NOs of the excited states are the same as those of the ground state.Conversely, in PNOC, this assumption is not required since the shape and occupancies of the perturbed NOs are projected to the NOs of the unperturbed system.It is interesting to notice that whereas the PNOC decompositions of static first and second hyperpolarizabilities are as simple as the PNOC splitting of the linear polarizability, the complexity of the SOS analysis rapidly increases with the order of property.
In our analysis, we choose the (U)CCSD orbital contributions as the reference values to assess the PNOC analysis at lower levels of theory.The total values of NLOPs computed using (U)CCSD(T) are also given for comparison.
According to PNOC, the difference between the overall contributions of π-and σ -type NOs to the linear response of benzene is smaller than 3% (see Table 2).As expected, the core orbitals are not involved in the static response (core electrons are highly confined near the nuclei).The largest contributions to CCSD α zz correspond to four frontier π NOs: 1e 1g , 2e 1g , 1e 2u and 2e 2u .
In the third-order response, the main role is shifted towards the π-type orbitals, for which the overall contribution to γ zzzz is about 60%.The most important contributions to α zz come from crossorbitals terms mixing occupied and virtual orbitals, whereas the dominant contributions to γ zzzz are due to terms involving only virtual orbitals.Indeed, the largest individual contributions belong to the frontier 2e 1g and 1e 2u orbitals, whereas the relative contribution of 1e 1g and 2e 2u is not as high as in the case of α zz .All occupied σ orbitals along with 1e 1g and 2e 2u π orbitals have a small negative contribution to the property.Indeed the largest contributions to γ zzzz come from the overall virtual σ and π orbitals.However, closer inspection reveals that individual higher virtual orbitals do not have large contributions, but there is a large number of small contributions adding to the total.
The formation of the singlet p-benzyne diradical changes the character of the electronic structure of the system, which has two unpaired electrons (of opposite spin) occupying 5b 1u and 6a g NOs with the occupancies 1.18 and 0.81, respectively (see Figure S.1 in ESI †).Such phenomenon is a direct manifestation of type a nondynamic correlation, 56 which results from an absolute neardegeneracy of HOMO (5b 1u ) and LUMO (6a g ) orbitals.In singlereference HF and post-HF methods (such as MP2 or CCSD), symmetry breaking allows to partially account for this type of electron correlation. 56][60][61][62] The formation of C 6 H 4 diradical from C 6 H 6 decreases the longitudinal CCSD α zz by approximately 5.3 a.u.However, the PNOC analysis reveals that in fact the contribution of the occupied σtype NOs increases by 6.0 a.u., and the lowering of the linear response is observed due to the virtual σ and occupied and virtual π-type orbitals.Within the latter group, the largest decrease is reported for the frontier π-type NOs, i.e., 2e 1g (1b 2g ) and 1e 2u (2b 3u ), respectively.
In contrast to the small decrease of α zz upon formation of the diradical, the total value of CCSD γ zzzz dramatically increases by 17148 a.u (140%).The PNOC decomposition shows that the nature of the response in p-benzyne (γ zzzz ) is completely dominated by the cross-orbital terms corresponding to two radical orbitals, 5b 1u and 6a g (the contribution of both NOs is 46% of the total γ zzzz ).Unlike C 6 H 6 , C 6 H 4 γ zzzz is almost fully dominated by the response of the σ -type NOs, and, although the larger overall part corresponds to the σ (vir), for C 6 H 4 the σ (occ) also plays a very important role.Conversely, π-type orbitals display a negative small contribution, coming mostly from the valence π NOs (see Table 2).
Next, we discuss the capability of PNOC to decompose the effect of the electron correlation in the description of static NLOPs.
On the examples of benzene and p-benzyne, we analyze the effect of the electronic correlation through the differences in the NO contributions to the properties obtained with the HF method and selected post-HF methods: MP2 and CCSD (which introduce dynamic correlation) and CASSCF (which introduces nondynamic correlation).HF, MP2 and CASSCF are compiled in Tables 3 and  4, whereas CCSD values, which serve as the reference, are taken from Table 2.
In benzene, the total α zz obtained from UHF is almost equal to the CCSD one.MP2 overestimates it by 2.0 a.u., whereas CASSCF underestimates it by 5.1 a.u.This suggests that UHF describes the linear response better than CASSCF, which by the variational principle yields energies closer to the exact value.PNOC reveals that CASSCF underestimates α zz due to a systematically smaller contribution of the π-type orbitals to α zz .On the other hand, MP2, which overestimates total α zz , actually describes the contributions for each of the σ and π valence orbitals much better than HF or CASSCF.The small error of HF total α zz is actually due to a fortunate compensation between an underestimation of the contribution of the σ orbitals and an overestimation of the contribution of the π orbitals.
In comparison to C 6 H 6 , the total value of α zz for C 6 H 4 is smaller by 15.6 a.u. in HF, 1.4 a.u. in MP2, and 7.4 a.u. in CASSCF (with respect to the CCSD value of 5.3 a.u.).Although the magnitude of the MP2 and CASSCF relative error of α zz with respect to UCCSD value is similar, changes in the nature of the linear response, such as the increase of the contribution of 5b 1u (and decrease of 6a g ) and overall decrease in contributions of the π-type NOs (especially frontier orbitals) are described better by CASSCF than UHF and MP2.UMP2 significantly overestimates the relative contribution of 1b 1g and 1a u NOs to the total α zz .The inaccuracy given by UHF and UMP2 partition of the linear response is in agreement with the important role of nondynamic correlation in p-benzyne.
Table 3 NOs contributions to longitudinal α zz (in a.u.) of benzene (C 6 H 6 ) and p-benzyne (C 6 H 4 ) obtained at different levels of theory.In the electronic structure computations, restricted and unrestricted HF/MP2 variants were used for C 6 H 6 and C 6 H 4 , whereas in the CASSCF/aug-cc-pVDZ computations active spaces of (6,6) and (8,8)  Differences in the description of the response are especially large for the second hyperpolarizability.In benzene, the γ zzzz total values obtained with HF and CASSCF are underestimated (with respect to CCSD ones).This drawback is found neither for MP2 or CCSD, which include dynamic correlation and provide a very similar value of γ zzzz .The PNOC also supports this conclusion: MP2 provides a well-balance description of the nonlinear response coming from σ , as well as π NOs (both occupied and virtual), Conversely, CASSCF strongly underestimates the overall contribution of the σ -type NOs and yields very inaccurate individual contributions of the valence π 1e 1g and 2e 2u NOs, revealing that the nonlinear response given by CAS (6,6) is still incorrect.
Upon formation of the diradical p-benzyne, dramatic discrepancies in the nature of γ zzzz are observed among these four methods.Interestingly, PNOC detailed analysis does not follow the trends expected from the total values of γ zzzz , i.e., in comparison to UCCSD, UMP2 provides the best total value for the nonlinear response, while UHF and CASSCF(8,8) perform unsatisfactory.However, PNOC analysis unveils that the small error of the MP2 total γ zzzz value of p-benzyne is accidental and comes from a fortunate cancellation of large errors.The largest discrepancy is observed for the response given by the 5b 1u and 6a g NOs (the only ones occupied by radical electrons), which incorrectly display an insignificant role in the nonlinear response.Similar errors are also observed for UHF decomposed NLOPs.Furthermore, according to the reference UCCSD results, relative contributions of π(all) should be -3.7% which UHF and UMP2 greatly overestimate (49.9% by UHF and 40.0% by UMP2).In contrast, the role of 5b 1u and 6a g NOs is correctly described by the CASSCF (8,8)  calculation.These two orbitals contribute approximately to 50% of the total CASSCF γ zzzz (against 46% in UCCSD).Moreover, at the CASSCF level of theory, the response of π NOs has dramatically decreased in the diradical (as compared to benzene) down to 8.3%, which is in much better agreement with UCCSD predictions.PNOC analysis shows that the error of CASSCF γ zzzz is again due to a strong (but systematic) underestimation of the overall contribution of the virtual σ -type NOs.This approach tacitly assumes a uniform contribution of the orbital to the property, as some of us have done in other contexts. 55ne of the most important properties of the latter definitions of α zz (r r r) and γ zzzz (r r r), is that the total values of α zz and γ zzzz are retrieved upon integration over the whole Cartesian space (likewise, the integration of the local p-th NO contribution yields the total NO contribution due to orbital p).This conclusion does not hold for density derivatives ρ (1) (r r r) and ρ (3) (r r r), which always integrate to zero.
In Figure 2, comparison of the information provided by α zz (r r r) and γ zzzz (r r r), and density derivatives ρ (1) (r r r) and ρ (3) (r r r) is shown for benzene and p-benzyne.Overall, the pictures provided by PNOC and density derivatives are in good agreement.However, PNOC representations are simpler and easier to analyze.For example, in C 6 H 6 , all carbons contribute equally to the polarizability and much more than hydrogens.Among those, the hydrogens in para position have slightly larger contribution to the longitudinal polarizability.γ zzzz (r r r) is mostly localized in the vicinity of the C1 and C4 atoms, and is showing a predominant π character.In C 6 H 4 , carbons with unpaired electrons have larger contribution to both the linear and nonlinear response, although the localization is far larger in γ zzzz .Additionally, for γ zzzz (r r r), one can see the small negative contributions of π type of orbitals at carbon atoms without unpaired electrons, in agreement with the results presented in Section 5.2.
Interestingly, a good qualitative correspondence is found between local NLOPs and the profiles of local nondynamic (I nd (r r r)) and dynamic (I d (r r r)) electron correlation of the unperturbed system. 55This can be easily seen in p-benzyne diradical, for which the regions with large I nd (r r r) approximately match γ zzzz (r r r).On the other hand, in benzene, γ zzzz (r r r) follows the profile of I d (r r r), but with larger contributions from the carbons aligned in the direction of the applied field.

Basis set dependency
PNOC proves useful also in the benchmark studies of the basis set effect in the NLOPs computations.In this section, we study the influence of diffuse functions in the character of the response using the PNOC analysis.We analyze a small π-conjugated system, all-trans-hexatriene (C 6 H 8 ), calculated with the CCSD method.Starting from the cc-pVDZ basis set, we separately add diffuse functions of selected angular momentum to the basis functions of carbon (C) or/and hydrogen (H), (see Section 4 for the description of these basis sets).We analyze the overall contributions of NOs belonging to two defined subspaces -full valence (FV) and higher virtual (HV) orbitals.The FV subspace consists of first 19 occupied (including core 1s C orbitals) and first 16 virtual σ -type NOs, and 3 occupied and first 3 virtual π-type NOs.All remaining virtual σ and π NOs are classified as the HV orbitals (and their number depends on the size of the basis set).
Results of calculated and decomposed NLOPs are compiled in Tables 5 and 6.PNOC reveals that the addition of diffuse functions to the basis sets increase the values of polarizability and second hyperpolarizability in a different manner.Changing from VDZ/VDZ to AVDZ/AVDZ, α zz rises only by 10%, 52% of this in-Fig.2 Comparison of functions representing the NLOPs in the Cartesian space for benzene (left) and p-benzyne (right): α zz (r r r) and γ zzzz (r r r) obtained from PNOC contributions, and density derivatives ρ (1) (r r r) and ρ (3) (r r r).Orange and green colors represent positive and negative values of the selected isosurfaces.The last row presents local nondynamic and dynamic correlation indices I nd (r r r) (light gray) and I d (r r r) (dark gray), 55 respectively.All the representations were obtained at the (U)CCSD/augcc-pVDZ level of theory.We have used the following isosurface values (in atomic units): α zz (r r r) = 0.40, γ zzzz (r r r) = 40.0,ρ (1) (r r r) = 0.085, ρ (3) (r r r) = 5.00, and I nd (r r r) = I d (r r r) = 0.006.crement coming from the contributions of the valence NOs (and the rest from the HV NOs).On the other hand, large changes of γ zzzz in hexatriene are mainly due to the larger contributions of the HV orbitals (which account for 72% of the property value).As expected, the second-order hyperpolarizability is more basis-set dependent than the polarizability.
The PNOC analysis can be used to make a judicious selection of the smallest basis set providing an accurate description of NLOPs.In the following example we aim at reproducing the PNOC analysis and the total values of the response calculated with AVDZ/AVDZ basis, with the smallest number of diffuse functions.The data compiled is in Tables 5 and 6.The addition of diffuse functions to carbon atoms is much more important than to the hydrogen centers for the correct description of both α zz and γ zzzz .Using AVDZ/VDZ basis, a very good agreement with AVDZ/AVDZ is found not only in the total values of α zz and γ zzzz , but also in the particular orbital components, namely σ (FV), σ (HV), π(FV) and π(HV).On the other hand, VDZ/AVDZ suffers mostly from the incorrect description of the response of the πtype NOs, which is especially important for the accurate description of γ zzzz .
Among the set of diffuse functions of C, p-type diffuse functions seem to be the most important for the description of both α zz and γ zzzz .They substantially enhance the response of the π-type NOs, especially the part coming from the HV NOs.On the other hand, s-type diffuse functions marginally affect the optical response in hexatriene (compare columns 2 and 3 or 4 and 6 in Tables 5 and  6).Interestingly, d-type diffuse functions enhance the total value of α zz by almost the same value as the p-type functions.However, PNOC unveils that the π(all) contribution to α zz is not so large, and this drawback is even magnified in γ zzzz .
Since the diffuse functions of p angular momentum in carbon seem to be crucial, we analyze the effect of adding diffuse functions to H centers too.The VDZ+p/VDZ+p basis set has substantially better performance than VDZ+p/VDZ and yields results comparable to AVDZ/VDZ.In the description of α zz , VDZ+p/VDZ+p corrects all orbital contributions to the response that were insufficiently described by VDZ+p/VDZ.The VDZ+p/VDZ+p basis correctly describes the particular γ zzzz components of both FV and HV NOs.Needles to say, the reduced size of this basis set could become in handy for the description of NLOPs in large molecules.

Conclusions
In this paper, we have presented an orbital decomposition scheme for optical properties, PNOC.This tool can be used to obtain a local representation of the NLOPs, providing a powerful visualization aid to connect the magnitude of the optical property with some parts of the molecule.This feature can be employed in the quest for nonlinear optical materials.
Perhaps the most interesting property of PNOC is the fact that the optical properties are decomposed in terms of the natural orbitals of the unperturbed system, making it a convenient method to detect flaws in NLOP evaluation of electronic structure methods.In particular, we have seen how the orbital contributions are more sensible to the choice of the electronic structure method than the total value of the optical property.In this way, PNOC unveils compensations of errors between orbital contributions and helps identifying orbitals that have a prominent role in the description of the property and, therefore, should be included in the active space of orbitals in a correlated calculation.For p-benzyne, PNOC has been also used to show that CASSCF provides a poor description of the total value γ zzzz because it underestimates the virtual σ orbital contribution, whereas MP2 calculations yield inaccurate values due to the lack of nondynamic correlation that causes the underestimation of the occupied σ orbitals.
Optical properties (especially of high order) are particularly susceptible to the inclusion of diffuse basis sets that significantly increase the computational cost of the calculation.The PNOC can be also used to identify the proper number and type of atomic basis sets required, saving substantial amounts of computer time for the calculation of optical properties in large molecules.
PNOC decomposes exactly the optical property regardless the approximation employed in the electronic structure method.Unlike other analysis methods for NLOPs, PNOC does not pose a large computational cost and it can be straightforwardly applied to any response order.For instance, the exact SOS analysis of the second-order hyperpolarizability involves large expressions and it is rarely applied.PNOC only requires the first-order reduced density matrix, which is easily available in most electronic structure methods.PNOC formulae can be easily generalized to evaluate any static or dynamic one-electron property, such as multipole moments and diamagnetic susceptibilities.
Table 1 Comparison among the orbital contributions to α zz in (H 2 ) 3 obtained from the SOS and the PNOC schemes using various wave functions.For each particular method, the absolute values (in a.u.) of the NOs contributions to α zz are presented (relative contributions given in parentheses).The full table and the orbitals are included in the ESI †(Tables S.I and Figure S1 J o u r n a l N a me , [ y e a r ] , [ v o l .] , 1-12 | 1

5. 3
Cartesian representation of NLOPs using PNOC Another useful feature of PNOC is its capability of representing the NLOPs in the Cartesian space without any additional computational cost.Within the PNOC scheme, a local representation of the NLOPs can be obtained through multiplication of each NO contribution to the property, i.e.. α zz,p or γ zzzz,p by its NO local amplitude, |χ p (r r r)| 2 : α zz (r r r) = ∑ p α zz,p χ p (r r r) 2 (16) γ zzzz (r r r) = ∑ p γ zzzz,p χ p (r r r) r n a l N a me , [ y e a r ] , [ v o l .] , J o u r n a l N a me , [ y e a r ] , [ v o l .] , 1-12 | 7 J o u r n a l N a me , [ y e a r ] , [ v o l .] , 1-12 | 11 7,[40][41][42][43]

Table 4
NOs contributions to longitudinal γ zzzz (in a.u.) of benzene (C 6 H 6 ) and p-benzyne (C 6 H 4 ) obtained at different levels of theory.In the electronic structure calculations, restricted and unrestricted HF/MP2 variants were used for C 6 H 6 and C 6 H 4 , whereas in the CASSCF/aug-cc-pVDZ computations active spaces of(6,6)and(8,8)were selected, respectively.

Table 2
Contributions of selected NOs to longitudinal α zz and γ zzzz of benzene (C 6 H 6 ) and p-benzyne (C 6 H 4 ), obtained at the (U)CCSD/aug-cc-pVDZ level of theory.Units are a.u. and the relative contributions are given in parenthesis.These NOs along with their occupancies are represented in the Cartesian space on Figure S.2 in ESI.†

Table 6
Basis set dependency of the NOs contributions to longitudinal γ zzzz of hexatriene.Decomposed values were obtained at the CCSD level of theory.Additionally, CCSD(T) total values are given.