Computation of Second Harmonic Generation for Crystalline Urea and

Subscriber access provided by University of Sussex Library

Article

Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham scheme Michel Rérat, Lorenzo Maschio, Bernard Kirtman, Bartolomeo Civalleri, and Roberto Dovesi J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b00791 • Publication Date (Web): 01 Dec 2015 Downloaded from http://pubs.acs.org on December 2, 2015

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

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

Page 1 of 24

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

Journal of Chemical Theory and Computation

Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham scheme. Michel Rérat1 ,∗ Lorenzo Maschio2 ,∗ Bernard Kirtman3 , Bartolomeo Civalleri2 , and Roberto Dovesi2 1

Equipe de Chimie Physique, IPREM UMR5254, Université de Pau et des Pays de l’Adour, 64000 Pau, France 2

Dipartimento di Chimica and NIS (Nanostructured Interfaces and Surfaces) Centre, Università di Torino, via Giuria 5, I-10125 Torino, Italy

3

Department of Chemistry and Biochemistry, University of California, Santa Barbara, California 93106 E-mail: [email protected]; [email protected]

∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


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

November 9, 2015 Abstract The electronic second harmonic generation (SHG) tensor, d, of crystalline urea and potassium dihydrogen phosphate (KDP) is evaluated as a function of frequency using a gaussian type basis set and the Coupled Perturbed Hartree-Fock (CPHF) and Kohn Sham (CPKS) schemes as implemented in the CRYSTAL code. The results of various functionals, including LDA, GGA (PBE), global and range-separated hybrids (B3LYP, PBE0, LC-BLYP), as well as Hartree-Fock, are compared. It is found that the calculated SHG intensity always decreases as the percentage of exact exchange increases. The hybrid functionals turn out to provide results that agree well with experiment. For urea and KDP the percentage of exact exchange determined by the inverse dielectric constant is too large. At 1064 nm the vibrational contribution for urea is found to be less than 5% of the total value. To the authors knowledge, this is the first coupled (self-consistent) calculation of SHG for any periodic system. Keywords: Urea and KDP, Second Harmonic Generation (SHG), non linear optics (NLO), Ab initio calculations, coupled-perturbed Hartree-Fock and Kohn-Sham method (CPHF/CPKS), DFT, Hartree-Fock hybrid functionals, CRYSTAL code

Introduction Second Harmonic Generation (SHG), also called frequency doubling, is one of the most technologically interesting nonlinear optical processes; it is exploited, for instance, in the generation of visible light lasers from infrared radiation, and in the SHG microscopy, widely used in medical and biological science. In SHG two phonons of frequency ω interact with matter to generate one single phonon of doubled frequency (2ω). The intensity of the effect is expressed by the third-order frequency-dependent nonlinear response of the material, χ(2) (−2ω; ω, ω).

2


Page 2 of 24

Page 3 of 24

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


Solid state SHG measurements are, typically, recorded with respect to crystalline urea or potassium dihydrogen phosphate (KDP) as a reference. For absolute values, or even relative values at different frequencies (ω), it is then necessary to know the absolute values for urea and KDP. However, only very limited results are available from experiment. In order to remedy this deficiency we have undertaken ab initio calculations using the analytical coupled perturbed Hartree-Fock and Kohn-Sham (CPHF/CPKS) method implemented in the CRYSTAL code (see Refs 1–4), which has very recently been generalized to include frequency-dependent first hyperpolarizabilities (β in Ref. 5). Both urea 6 and KDP 7 have been considered previously, although only in the static limit. Thus, we build on what has already been learned, particularly with regard to structure, basis set and computational parameters, as discussed in the following section. The choice of method is especially important. In that connection, calculations were carried out using Hartree-Fock (HF) as well as pure DFT and hybrid DFT/HF functionals from which the role of exact exchange, among other things, can be assessed. In the UV-visible range of energy, the vibrational (ionic) contribution to the SHG electric susceptibility is negligible. However, the Nd:YAG laser with a wavelength in the far infrared at 1064 nm is often used experimentally. For that case, the most important ionic contribution to SHG, i.e. the double harmonic term, 8 has also been taken into account in order to improve the agreement with experiment at low photon energy. In Section II, we summarize the frequency-dependent formulas along with the computational conditions. Then the coupled perturbed results are considered in Section III as well as a comparison of these results with simpler Sum-Over-States (SOS) calculations. Finally, Sec. IV contains a summary of our findings.

3



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

Page 4 of 24

Coupled Perturbed HF/KS calculation of the dynamic first hyperpolarizability We have very recently 5 implemented the fully analytical calculation of coupled perturbed Hartree-Fock and Kohn-Sham (CPHF/CPKS) dynamic electronic (clamped nuclei) hyperpolarizabilities in the CRYSTAL code. 9 Our treatment is based on the 2n + 1-rule, which requires solution of only first-order perturbation equations. In addition to the electronic value, calculated by the CPHF/CPKS method, there is also an ionic contribution to SHG. It vanishes in the infinite frequency limit, but may become important as the frequency is lowered. We evaluate this contribution in the so-called double harmonic approximation following the prescription of Bishop and Kirtman. 8 For sake of completeness an abbreviated version of the expressions employed for both contributions is provided in this section.

The CPHF 2n + 1-rule expression for the electronic SHG of closed shell periodic systems may be written as

dtuv(−ωσ ;ω1 ,ω2 ) = −

1 2π 2 X X X ℜ P −ωt , +ωu , +ωv × σ 1 2 2 V nk a i ~k

×

  



(t)∗ Uai(−ωσ ) (~k) 

X b

(v) (u) Gab(+ω1 ) (~k)Ubi(+ω2 ) (~k) −

X j

 (v) ~k)  ( ∂U ai(+ω2 ) (u) (v)  Uaj(+ω2 ) (~k)Gji(+ω1 ) (~k) + ı  ∂ku

(1)

where i, j and a, b run over occupied and virtual crystalline orbitals (CO) per ~k-point, respectively, and V is the unit cell volume. Here nk indicates the number of k points sampled in the first Brillouin zone, ℜ indicates the real part and P is the permutation operator. The U-matrix in Eq. (1) determines the linear response of the crystal orbitals to the external electric field(s) which, in turn, determines χ(2) . In Eq. (1) the off-diagonal virt-occ (u) blocks of the U matrix, e.g. Uai(+ω1 ) (~k) (the diagonal blocks are chosen to vanish in the

4


Page 5 of 24

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


non-canonical treatment employed), are defined by: Gai(+ω1 ) (~k) =− Ea (~k) − Ei (~k) + ω1 (u)

(u) Uai(+ω1 ) (~k)

(2)

(u) Here Gai(+ω1 ) (~k) is obtained from the first-order atomic orbital (AO) Fock matrix in reciprocal

space (u) Gai(+ω1 ) (~k)

Xh

= Ξai (~k) +

(u) ∗ ~ Cµa (k)Fµ,ν(+ω1 ) (~k)Cνi (~k)

µ,ν

i

(3)

in which Ξai (~k) = ha(~k)|r + ı∇~k |i(~k)i is the transition moment matrix element between (u) occupied i(~k) and virtual a(~k) crystalline orbitals, 10,11 whereas Fµ,ν(+ω1 ) (~k) is calculated by a

Fourier transformation of the 2-electron AO integrals after contraction with the perturbed density matrix,i.e. Fµ,ν(+ω1 ) (~k) = (u)

X ~g

~

e−ık·~g

XX ~h,~ n ρ,τ

~ ~ (u) µν ~g ||ρh τ h+~n Dρ,τ ~n (+ω1 )

(4)

In this last expression ~g , ~h and ~n are lattice vectors, so that the notation τ ~n indicates that the atomic orbital τ is centered in cell ~n. On the other hand, the density matrix in Eq.(4) u Dρ,τ ~ n (+ω ) is determined through a back-Fourier transform of its reciprocal space image, which 1

is given for closed shells by

(u)

Dµ,ν~g (+ω1 ) = 2

X ~k

~



eık·~g 

X

(u) ∗ ~ (k) + Cµa (~k) Uai(+ω1 ) (~k) Cνi

ia

X ia



(u)∗ ∗ ~  (5) (k) Cµi (~k) Uai(+ω1 ) (~k)Cνa

This procedure is designed so as to avoid the costly direct transformation of 2-electron integrals from AOs to crystalline orbitals.

The extension to KS-DFT has also been carried out in the previous work (Ref. 5) and implemented in the CRYSTAL code, so that hybrid DFT/HF Hamiltonians can be, and are, used in the following application. We do not repeat the DFT expression here since it is quite 5



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

Page 6 of 24

messy (see Ref. 5), even though straightforward to evaluate. The double harmonic vibrational [µα](0,0) (or ionic) contribution, although not previously evaluated for any periodic system, is obtained by a simple generalization of the molecular formula given in Ref. 8:

dvib tuv (−ωσ ; ω1 , ω2 ) =

1 2π 2V

3N −6 X i=1

P −ωt

, u , v σ +ω1 +ω2

∂µt ∂Qi

∂αuv ∂Qi

ωi2 − ωσ2

(6)

where µt and αuv are dipole moment and polarizability components, while ωi and Qi are phonon frequencies and normal modes of the i-th vibration at the Γ-point. The terms and

∂αuv ∂Qi

∂µt ∂Qi

are the same tensors needed for the evaluation of infrared and Raman intensities

and are computed as described in Refs. 12 and 13, respectively. There are other vibrational contributions which depend upon electrical and mechanical anharmonicity. However, the double harmonic term is expected to be dominant, especially for urea. For KDP, anharmonicity may play a significant role in connection with H-bond tunneling but we have used an averaged room temperature structure and, thus, a much more sophisticated treatment would be required to treat that case.

Computational details Structure and symmetry Urea is a non-centrosymmetric tetragonal molecular crystal that has P¯421 m space group symmetry and contains two molecules per unit cell (16 atoms/cell). At room temperature KDP is tetragonal (I4d2 space group), with two molecules in the unit cell as well (also in this case there are 16 atoms/cell). Due to symmetry, both urea and KDP have only two independent non-vanishing components of the dielectric constant ǫ (ǫxx and ǫzz ) and just one such component of the first-order (2)

nonlinear susceptibility tensor χ(2) (χxyz ) at the static limit. On the other hand, for SHG 6


Page 7 of 24

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


there are two independent dynamic susceptibility (d = 12 χ(2) ) components: dxyz (−2ω; ω, ω) and dzxy (−2ω; ω, ω). It is well known that DFT functionals overestimate the lattice parameters of molecular crystals, 14 unless empirical dispersion corrections are used, in which case the cell volume is underestimated. 15 For that reason we used the experimental values for urea: 16 a=5.565 Å and c=4.684 Å at 12K, and for KDP: 17 a=7.444 Å and c=6.967 Å at room temperature. In the above symmetry the hydrogen atoms of KDP are at the midpoint between the O1 and O2 oxygen atoms of the H-bonds. The internal coordinates of both crystals were then optimized at the B3LYP level of theory, with a triple-ζ basis as specified in the following.

Hamiltonians and basis set A selection of LDA, GGA and global hybrid exchange-correlation functionals were used in this study, namely: S-VWN 18,19 (hereafter referred to as LDA), PBE, 20 PBE0 21 and B3LYP. 22–24 In addition, to assess the complete replacement of DFT exchange by HF exchange at longrange, the long-range corrected (LC) hybrid functional LC-BLYP proposed by Hirao and co-workers 25,26 was considered. A value of µ=0.33 Å−1 for the range separation parameter was used since it has been proved to give good performance for the prediction of response properties. 27 The same basis set was used as in our previous studies of urea 6 and KDP 7 to which we refer for the effect of various basis sets on the static hyperpolarizability tensor. It is an all-electron basis of triple-ζ quality enriched with a double set of polarization functions. For urea, it corresponds to the basis set denoted as TZPP in ref. 6, while for KDP it is the one labeled as basis set D in ref. 7. A careful detailed discussion of the rationale for this choice of basis set is given in the cited papers.

7



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

Parameters of the calculation The present calculations were performed with a development version of the periodic ab initio CRYSTAL14 code, 9 that uses a Gaussian basis set. Evaluation of the Coulomb and exact exchange infinite series is controlled by five parameters 28 (T1 , T2 , T3 , T4 , T5 ), whose values are set to T1 = T2 = T3 = T4 = TI and T5 = 2 ∗ T4 , with TI = 8 (in the HF case tighter exchange tolerances have been used: T4 = 14 and T5 = 28). The shrinking factor used to generate the commensurate net of points at which the Fock matrix is diagonalized and the CPHF/CPKS equations are solved is S = 4. Convergence of the field-free and first-order SCF cycles is controlled by the parameters TE = 10−10 (total energy difference between two subsequent field-free cycles, in Hartree) and TCP = 10−4 (polarizability difference between two subsequent first-order SCF cycles in a.u.). Various schemes are used for damping oscillations along the SCF cycles.

Results Comparison with experiment Very few experimental absolute values of the SHG susceptibility are found in the literature. Generally, the NLO d tensor is compared to that of urea or KDP. Based on our previous success with B3LYP 6,7 at the static limit, we first compare experimental data with values obtained using that functional, while a comparison with other Hamiltonians will be made in the following sub-section. Computed results for urea and KDP (see Table 2) show that B3LYP reproduces rather well the experimental values measured at 1064 and 600 nm. For urea, the quasi-isotropic electronic contribution is slightly smaller than the experimental d14 -value 29 at λ = 1064 nm. At this wavelength the vibrational contribution may not be completely negligible. For that reason the double harmonic vibrational term was calculated (see Eq. 6); its value is

8


Page 8 of 24

Page 9 of 24

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


0.041 and 0.035 pmV −1 for the xyz and zxy-component, respectively, and has the same sign as the electronic contribution. The total B3LYP value is, then, equal to 1.027 and 1.018 pmV −1 which falls essentially at the outer limit of the error bars for the experimental value: d14 = 1.2 ± 0.1 pmV −1 of Ref. 29 (dxyz ≃ dzxy for wavelength larger than 600 nm). Vibrational anharmonicity and/or temperature effects, which would increase the magnitude of this term, could be part of the reason for the small discrepancy. At 600 nm the calculated vibrational contribution is four times smaller than at 1064 nm and, thus, can be neglected. The B3LYP value (1.371 and 1.361) in this case is well within the experimental window: 1.3 ± 0.3 pmV −1 of Ref. 30. Finally, for KDP the B3LYP electronic value of 0.41 pmV −1 is in good agreement with the available experimental datum (0.41 pmV −1 ). 31 Although there may be a significant vibrational contribution that has not been accounted for in the case of KDP, the agreement with experiment for both crystals gives us some confidence that the B3LYP functional provides good predictions for the frequency-dependent SHG.

Influence of the Hamiltonian In Table 1 we have reported the values for the electronic dielectric tensor obtained with different Hamiltonians. Overall, DFT methods are far better than HF. Quite surprisingly, for this quantity, PBE and LDA seem to outperform hybrid functionals in the case of urea, although a better agreement is observed for ǫxx than ǫzz . It was shown in ref. 6 that the prediction of ǫxx is more delicate than for ǫzz , with former being more sensitive to the basis set. Indeed, results obtained at the static limit with the TZPP basis set were slightly underestimated with respect to even larger basis sets. 6 For instance, ǫxx =2.142 and ǫzz =2.465 are predicted at the B3LYP/QZVPPP level. Therefore, we expect that LDA and PBE would overestimate both components of the dielectric tensor if a more complete basis set were utilized. For the KDP crystal, the situation is quite different whether one looks at ǫxx or ǫzz : the former seems to be better described by hybrid functionals, in particular B3LYP, while 9



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

Page 10 of 24

the latter is better matched by PBE and LDA values. In Figure 1, the SHG dxyz component of urea has been plotted as a function of the wavelength for different hamiltonians. The other non-vanishing component, dzxy , has a similar dispersion behavior and has not been plotted. We can see in the figure that the curves for the various methods remain practically parallel except as one approaches the first resonance.

Results in Table 2 and Figure 1 clearly show the well-known tendency of LDA and GGA functionals to grossly overestimate high-order electric susceptibilities as far as the value of the wavelength approaches the resonance. It is noteworthy that dxyz increases at each frequency when the percentage of HF-exchange decreases. This is not unexpected because it correlates with the predicted band gap for the different Hamiltonians which decreases from HF (14.0 eV, 89 nm) to PBE0 (7.4 eV, 168 nm) and B3LYP (6.9 eV, 180 nm) to PBE (5.2 eV, 239 nm) and LDA (4.8 eV, 257 nm) It is interesting that LC-BLYP predicts a band gap of 7.7 eV (161 nm) and shows a slightly different trend with respect to global hybrids. This may be related to the fact that LC hybrids fully recover the correct decay of the exchange potential at long-range, while for global hybrids the long-range decay depends on the amount of HF exchange included.

In passing, we point out that if the inverse of the average static dielectric constant ǫ is used to evaluate the amount of HF-exchange in global hybrid functionals, as proposed in Ref. 32, it would lead to an apparently too large SHG value. Indeed, our calculated 1/ǫ for urea would imply a coefficient for HF-exchange that is larger than 0.45. This, in turn, would result in an SHG well below PBE0 and B3LYP, which would be much too small compared to experiment. For instance, in the case of urea crystal, the dxyz values obtained using a modified PBE0 with 45% HF exchange are 0.778, 0.868 and 1.133 pm/V at ω = ∞, 1064 and 600 nm, respectively. For the same frequencies, the corresponding results for KDP are 0.292, 0.309 and 0.350 pm/V.

10


Page 11 of 24

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


The frequency of the lowest allowed transition in urea (see Fig.1 caption for method of estimation) also increases as the percent HF-exchange decreases with estimated values of 121, 154, 162, 192 and 198 nm for HF, PBE0, B3LYP, PBE and LDA, respectively. LC-BLYP does not have a well-defined percent of HF-exchange; the corresponding estimate for the first transition is 150, ranking between HF and PBE0. Virtually identical pole values are obtained for the dzxy component: 122 (HF), 150 (LC-BLYP), 154 (PBE0), 163 (B3LYP), 192 (PBE) and 198 (LDA) nm. Not unexpectedly, an almost linear correlation exists between the fundamental direct band gap computed as the difference between the HOCO and LUCO electronic levels at ~k = ~0 (see above) and the value of the pole, which is related to the optical band gap. The corresponding dispersion results for KDP are also reported in Table 2. Again, both the dxyz and dzxy components increase at each frequency when the percentage of HF-exchange in the Hamiltonian decreases. In this case LC-BLYP ranks in between PBE0 and B3LYP. In fact, all the conclusions drawn above concerning urea apply equally as well to KDP.

Comparison of CPKS and Sum-Over-States (SOS) for B3LYP Figures 2 and 3 show a comparison between the CPKS and SOS frequency dispersion curves, obtained with the B3LYP functional, for the dxyz and dzxy components of urea and KDP, respectively. The horizontal line is the static limit. Note that the two curves (CPKS and SOS) are nearly parallel (except near the first resonance) even though the SOS values of d are unsatisfactory.

For both urea and KDP, the dxyz and dzxy values are quite similar with dxyz increasing slightly more rapidly near the first resonance. Thus, the calculated transition frequency 2ω associated with the x component of the dipole moment (in dxyz ) occurs at a slightly lower frequency than that associated with the z component (in dzxy ), just as the 11



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

Page 12 of 24

frequency associated with the resonance of ǫxx is less than the corresponding frequency of the resonance of ǫzz (see below). A fit of the SHG frequencies to d(λ = 2π/ω) = 1 1 1 a + b (ω0 −2ω)(ω for each component and for each system + + (ω0 +ω)(ω0 −ω) (ω0 +ω)(ω0 +2ω) 0 −ω) allows us to obtain the half frequency of the first resonance. Here the constant term is introduced to take into account higher excited states. An accurate fit of the B3LYP values leads to: λ0 = 2π/ω0 = 165 (154) and 163 (146) nm for the xyz- and zxy-components of urea (the values in parentheses are the SOS results). For KDP the corresponding results are 114 (103) and 112 (105) nm. For comparison, we have also carried out a fit for the relative dielectric constant of urea based on the formula ǫr (λ) = a + b/(ω02 − ω 2 ), which gives: λ0 = 168 (159) and 165 (155) nm for the xx- and zz-components (Fig. 4). The differences from the results obtained for SHG are relatively small: 3 (5) and 2 (9) nm, respectively. Although a more complete study near resonance would undoubtedly reduce these differences, the predicted values are not far from the observed experimental data for λ0 = 173 nm from ref. 33 and λ0 = 176 nm from polarized UV spectra of urea crystal. 34

Conclusions A very recent implementation in the CRYSTAL code of the analytical CPHF/CPKS method for computing dynamic electronic hyperpolarizabilities of periodic systems has been used to obtain the frequency-dependent SHG tensor of the reference compounds urea and KDP. All-electron calculations in a local basis set have been carried out with the Hartree-Fock Hamiltonian, as well as the LDA, PBE, PBE0, LC-BLYP and B3LYP functionals of density functional theory. The main conclusions of this work are the following: • There is good agreement between the few measured SHG values and global hybrid as well as long-range corrected functionals for urea (two frequencies) and KDP (one frequency). 12


Page 13 of 24

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


• The vibrational (or ionic) contribution to SHG is expected to be negligible at short wavelengths, but not necessarily at long wavelengths. A B3LYP test calculation for urea at 1064 nm, using the double harmonic approximation, gave a relative contribution (as compared to the pure electronic term) of less than 5%. • The calculated electronic SHG intensity always decreases as the percentage of exact exchange increases. We considered the recent proposal 32 to use the inverse dielectric constant as the optimum percentage to determine the optical gap of condensed systems (and the dielectric constant itself through an iterative procedure) for prediction of χ(2) as well. In the case of urea and KDP the results are not as good as those obtained with conventional hybrids. We believe, however, that more comprehensive tests, which we plan to undertake are needed to draw a firm conclusion. • In the case of B3LYP for example, the dispersion curves obtained by the simpler uncoupled, or sum-over-states, procedure parallel those of the coupled perturbed treatment at frequencies sufficiently below the first resonance. However, the absolute SHG values and the position of the resonance are significantly different. A similar behavior is observed for the dielectric constant of urea. • Overall, the CPKS method combined with both global and long-range corrected hybrid functionals provides a reliable SHG dispersion curve, at least for the two reference compounds considered here. The present validation of the frequency-dependent CPKS method implemented in the CRYSTAL code, as applied to two well-known 3D-systems allows us to consider promising new NLO nanomaterials of varying dimension and nature in future work.

13



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

Page 14 of 24

Acknowledgments References (1) Ferrero, M.; Rérat, M.; Kirtman, B.; Dovesi, R. Calculation of First and Second Static Hyperpolarizabilities of One- to Three-Dimensional Periodic Compounds. Implementation in the Crystal Code. J. Chem. Phys. 2008, 129, 244110. (2) Ferrero, M.; Rérat, M.; Orlando, R.; Dovesi, R. The Calculation of Static Polarizabilities of 1-3D Periodic Compounds. the Implementation in the Crystal Code. J. Comp. Chem. 2008, 29, 1450. (3) Ferrero, M.; Rérat, M.; Orlando, R.; Dovesi, R. Coupled Perturbed Hartree-Fock for Periodic Systems: The role of Symmetry and Related Computational Aspects. J. Chem. Phys. 2008, 128, 014110. (4) Orlando, R.; Lacivita, V.; Bast, R.; Ruud, K. Calculation of the First Static Hyperpolarizability Tensor of Three-Dimensional Periodic Compounds with a Local Basis Set: A Comparison of LDA, PBE, PBE0, B3LYP, and HF Results. J. Chem. Phys. 2010, 132, 244106. (5) Maschio, L.; Rérat, M.; Kirtman, B.; Dovesi, R. Calculation of the Dynamic First Hyperpolarizabilty β(−ωσ ; ω1 , ω2 ) of Periodic Systems. Implementation in the CRYSTAL Code. J. Chem. Phys. 2015, submitted . (6) Ferrero, M.; Civalleri, B.; Rérat, M.; Orlando, R.; Dovesi, R. The Calculation of the Static First and Second Susceptibilities of Crystalline Urea: A Comparison of HartreeFock and Density Functional Theory Results Obtained with the Periodic Coupled Perturbed Hartree-Fock/Kohn-Sham Scheme. J. Chem. Phys. 2009, 131, 214704. (7) Lacivita, V.; Rérat, M.; Kirtman, B.; Ferrero, M.; Orlando, R.; Dovesi, R. Calculation of the dielectric constant ǫ and first nonlinear susceptibility χ(2) of Crystalline KDP by 14


Page 15 of 24

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


the CPHF and CPKS schemes as implemented in the CRYSTAL code. J. Chem. Phys. 2009, 131, 204509. (8) Bishop, D. M.; Kirtman, B. A Perturbation Method for Calculating Vibrational Dynamic Dipole Polarizabilities and Hyperpolarizabilities. J. Chem. Phys. 1991, 95, 2646. (9) Dovesi, R.; Orlando, R.; Zicovich-Wilson, C. M.; Civalleri, B.; Maschio, L.; Erba, A.; Casassa, S.; Ferrabone, M.; Pierre, M. D. L.; D’Arco, P.; Noël, Y.; Causà, M.; Rérat, M.; Kirtman, B. Int. J. Quant. Chem. 2014, 114, 1287. (10) Otto, P. Calculation of the Polarizability and Hyperpolarizabilities of Periodic QuasiOne-Dimensional Systems. Phys. Rev. B 1992, 45, 10876. (11) Otto, P.; Gu, F. L.; Ladik, J. Calculation of Ab Initio Dynamic Hyperpolarizabilities of Polymers. J. Chem. Phys. 1999, 110, 2717. (12) Maschio, L.; Kirtman, B.; Orlando, R.; Rérat, M. Ab Initio Analytical Infrared Intensities for Periodic Systems through a Coupled Perturbed Hartree-Fock/Kohn-Sham method. J. Chem. Phys. 2012, 137, 204113. (13) Maschio, L.; Kirtman, B.; R´erat, M.; Orlando, R.; Dovesi, R. Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis. I. Theory. J. Chem. Phys. 2013, 139, 164101. (14) Civalleri, B.; Doll, K.; Zicovich-Wilson, C. M. Ab-initio investigation of structure and cohesive energy of crystalline urea. J. Phys. Chem. B 2007, 111, 26. (15) Civalleri, B.; Zicovich-Wilson, C.; Valenzano, L.; Ugliengo, P. B3LYP Augmented with an Empirical Dispersion Term (B3LYP-D*) as Applied to Molecular Crystals. Cryst. Eng. Comm 2008, 10, 405. (16) Swaminathan, S.; Craven, B. N.; McMullan, R. K. The Crystal Structure and Molecular

15



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

Page 16 of 24

Thermal Motion of Urea at 12, 60 and 123 K from Neutron Diffraction. Acta Cryst. B 1984, 40, 300. (17) See http://icsd.fiz-karlsruhe.de/ accessed July 2015. (18) Dirac, P. A. M. Note on Exchange Phenomena in the Thomas Atom. Proc. Cambridge Philos. Soc. 1930, 26, 376. (19) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis. Can. J. Phys. 1980, 58, 1200. (20) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. 1996, 77, 3865. (21) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: the PBE0 model. J. Chem. Phys. 1999, 110, 6158–6170. (22) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098. (23) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. (24) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648. (25) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A long-range correction scheme for generalized-gradient approximation exchange functionals. J. Chem. Phys. 2001, 115, 3540–3544. (26) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A long-range-corrected time-dependent density functional theory. J. Chem. Phys. 2004, 120, 8425–8433. 16


Page 17 of 24

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


(27) Tsuneda, T.; Hirao, K. Long-range correction for density functional theory. WIREs Comput. Mol. Sci. 2014, 4, 375–390. (28) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; Zicovich-Wilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M.; Causà, M.; Noël, Y. CRYSTAL14 User’s Manual, University of Torino 2014, (29) Halbout, J. M.; Blit, S.; Donaldson, W.; Tang, C. L. Efficient Phase-Matched SecondHarmonic Generation and Sum-Frequency Mixing in Urea. IEEE J. Quantum Electron. 1979, 15, 1176. (30) Bäuerele, D.; Beltzer, K.; Hesse, H.; Kaplan, S.; Loose, P. Phase-Matched SecondHarmonic-Generation in Urea. Phys. Stat. Sol. 1977, 42, K119. (31) Singh, S. Handbook of Laser Science; Ed. M. J. Weber, CRC, Boca Raton, Florida, 1986; pp 3–228. (32) Skone, J. H.; Govoni, M.; Galli, G. Self-Consistent Hybrid Functional for Condensed Systems. Phys. Rev. B 2014, 89, 195112–12. (33) Rosker, M. J.; Cheng, K.; Tang, C. L. Practical urea optical parametric oscillator for tunable generation throughout the visible and near-infrared. 1985, QE-21, 1600. (34) Campbell, B. F.; Clark, L. B. Polarized Vacuum Ultraviolet Spectra of Crystalline Urea. J. Am. Chem. Soc. 1989, 111, 8131–8136. (35) Rosker, M. J.; Cheng, K.; Tang, C. L. Practical Urea Optical Parametric Oscillator for Tunable Generation Throughout the Visible and Near Infrared. IEEE J. Quant. Electron. 1985, 21, 1600–1606. (36) See http://refractiveindex.info/?shelf=organic&book=urea&page=Rosker-o accessed June 2015.

17



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

Page 18 of 24

(37) Levine, Z. H.; Allan, D. C. Large Local-Field Effects in the Second-Harmonic Susceptibility of Crystalline Urea. Phys. Rev. B 1993, 48, 7783. (38) Eckardt, R. C.; Byer, R. L. Measurement of Nonlinear Optical Coefficients by PhaseMatched Harmonic Generation; SPIE, Inorganic Crystals for Optics, Electro-Optics and Frequency Conversion, 1991; Vol. 1561; p 119.

18


Page 19 of 24

4

HF PBE0 B3LYP LC-BLYP PBE LDA Experiment

3 dxyz (pm/V)

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


2

1

300

500

700

900

1100

1300

λ (nm)

Figure 1: SHG dxyz tensor component of bulk urea obtained with various functiond(λ) = a + als. The computed data are fit to the three parameter function: b 1/(ω0 − 2ω)(ω0 − ω) + 1/(ω0 + ω)(ω0 − ω) + 1/(ω0 + ω)(ω0 + 2ω) .

19



4.5 CPKS static SOS static CPKS xyz SOS xyz CPKS zxy SOS zxy

4 3.5 3 d (pm/V)

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

Page 20 of 24

2.5 2 1.5 1 0.5 300

500

700

900

1100

1300

λ (nm)

Figure 2: Dispersion of the SHG dxyz and dzxy tensor components of bulk urea obtained with the B3LYP functional. The four curves refer to the uncoupled (SOS=Sum-Over-States) and coupled perturbed Kohn-Sham (CPKS) results for the two components (cfr. Table 2).

20


Page 21 of 24

0.59 CPKS static SOS static CPKS xyz SOS xyz CPKS zxy SOS zxy

0.54

d (pm/V)

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


0.49

0.44

0.39

0.34 300

500

700

900

1100

1300

λ (nm)

Figure 3: Dispersion of the SHG dxyz and dzxy tensor components of bulk KDP obtained with the B3LYP functional. The four curves refer to the uncoupled (or SOS=Sum-Over-States) and coupled (CPKS) results for the two components (cfr. Table 2).

21



CPKS xx SOS xx CPKS zz SOS zz static CPKS xx static SOS xx static CPKS zz static SOS zz exp. o exp. e

3.4 3.2 3 2.8 ε

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

Page 22 of 24

2.6 2.4 2.2 2 100

300

500

700 λ (nm)

900

1100

1300

Figure 4: Dispersion of the dielectric ǫxx and ǫzz tensor components of bulk urea obtained with the B3LYP functional as compared to the ordinary (o) and extraordinary (e) experimental values. The four calculated curves refer to the electronic uncoupled (SOS) and Coupled Perturbed Kohn-Sham (CPKS) results for the xx and zz components. We fit the computed values using the function ǫ(λ) = a + b/(ω02 − ω 2 ) .

22


Page 23 of 24

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


Table 1: Relative dielectric ǫxx /ǫzz tensor components of urea and KDP (SOS values in parentheses) computed using various hamiltonians. λ (nm) urea ∞

HF

PBE0

B3LYP

LC-BLYP

PBE

LDA

Exp.

1.901/2.187

2.059/2.433

2.057/2.430

2.149/2.567

2.186/2.599

2.182/2.515a,b

1064

1.907/2.195

2.070/2.448

2.067/2.446

2.163/2.588

2.202/2.621

2.194/2.529a,b

600

1.920/2.214

2.093/2.483

2.070/2.451 (2.128/2.455) 2.081/2.467 (2.138/2.468) 2.106/2.504 (2.160/2.498)

2.091/2.481

2.195/2.635

2.237/2.671

2.220/2.577a,b

KDP ∞

2.025/1.868

2.223/2.039

2.219/2.049

2.327/2.132

2.341/2.160

2.259/2.133b

1064

2.030/1.872

2.231/2.046

2.227/2.055

2.339/2.141

2.352/2.169

2.231/2.131b

600

2.041/1.880

2.250/2.060

2.230/2.046 (2.244/2.086) 2.239/2.053 (2.251/2.091) 2.258/2.068 (2.284/2.117)

2.245/2.069

2.363/2.160

2.377/2.189

2.277/2.155b

a: Ref. 35 b: Ref. 36

Table 2: SHG electric susceptibility dxyz /dzxy (in pm/V) of urea and KDP (SOS values in parentheses) computed using various hamiltonians. λ (nm) urea ∞

HF

PBE0

B3LYP

LC-BLYP

PBE

LDA

0.680/0.680

0.823/0.823

0.876/0.876

0.949/0.949

1.106/1.106

1064

0.738/0.737

0.936/0.934

0.989/0.988

1.131/1.128

1.333/1.329

1.2 ± 0.1a,b

600

0.894/0.889

1.279/1.263

0.863/0.863 (1.401/1.401) 0.986/0.983 (1.564/1.555) 1.371/1.361 (2.042/1.997)

1.333/1.323

1.824/1.793

2.243/2.191

1.3 ± 0.3a,c

KDP ∞

0.198/0.198

0.354/0.354

0.342/0.342

0.467/0.467

0.487/0.487

1064

0.207/0.207

0.378/0.378

0.363/0.363

0.505/0.504

0.530/0.529

600

0.228/0.227

0.438/0.435

0.373/0.373 (0.341/0.341) 0.397/0.396 (0.361/0.361) 0.464/0.460 (0.412/0.412)

0.416/0.415

0.607/0.603

0.642/0.636

a: b: c: d: e:

Ref. Ref. Ref. Ref. Ref.

37 29 30 38 and references therein 31

23


Exp.

0.38d 0.41e


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

Page 24 of 24

Computation of Second Harmonic Generation for Crystalline Urea and

Recommend Documents