The First and Second Static Electronic Hyperpolarizabilities of Zigzag

ARTICLE pubs.acs.org/JPCA

The First and Second Static Electronic Hyperpolarizabilities of Zigzag Boron Nitride Nanotubes. An ab Initio Approach through the Coupled Perturbed Kohn
Sham Scheme Roberto Orlando,*,† Radovan Bast,‡ Kenneth Ruud,‡ Ulf Ekstr€om,§ Matteo Ferrabone,|| Bernard Kirtman,^ and Roberto Dovesi|| †

Dipartimento di Scienze e Tecnologie Avanzate, Universita del Piemonte Orientale, Viale T. Michel 11, 15121 Alessandria, Italy, Centre for Theoretical and Computational Chemistry (CTCC), Department of Chemistry, University of Tromsø, Tromsø N.9037, Norway, § Centre for Theoretical and Computational Chemistry (CTCC), Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N.0315 Oslo, Norway, Dipartimento di Chimica IFM, Universita di Torino and NIS - Nanostructured Interfaces and Surfaces - Centre of Excellence, http://www.nis.unito.it, Via P. Giuria 7, 10125 Torino, Italy, and ^ Department of Chemistry and Biochemistry, University of California, Santa Barbara, California 93106, United States

)

‡

ABSTRACT: The coupled perturbed Kohn
Sham (CPKS) computational scheme for the evaluation of electric susceptibility tensors in periodic systems, recently implemented in the CRYSTAL code, has been extended to third-order. It is, then, used to obtain static electronic hyperpolarizabilities of zigzag BN nanotubes for the first time. This procedure, which is fully analytic in all key steps, requires a double self-consistent treatment for taking into account the first- and second-order response of the system to the applied field. The performance of different functionals is compared and the B3LYP hybrid is ultimately chosen for calculations on nanotubes having radii as large as R = 20 Å (6
200 atoms in the unit cell). Such large radii are sufficient to give the pure longitudinal component of the (hyper)polarizability tensors to within 1% of the “exact” hexagonal BN monolayer limit. Other tensor components involving the transverse direction converge more slowly. They can, however, be extrapolated to the monolayer limit to within 4% accuracy except for the pure transverse second hyperpolarizability, which has an error of 13% in that limit.

1. INTRODUCTION The existence of boron nitride nanotubes (BNNT) was proposed theoretically1,2 in 1994, on the basis of the analogy between graphite and hexagonal boron nitride (h-BN), and was first synthetized soon after,3 both as multiwalled (MWBNNT) and as single-walled nanotubes (SWBNNT). Though h-BN and graphite are isoelectronic and have very similar structures, the corresponding nanotubes show significant differences in their properties. In particular, BNNTs have greater thermal stability4 and, thanks to a wide band gap (=5.9 eV), a less dramatic dependence of the electrical properties on the rolling direction and tube diameter. The stability of the properties with respect to tube size is an important advantage of BNNTs over carbon nanotubes. For the latter, poor control of tube chirality and tube size in their synthesis leads to poor control of tube properties. As a result, BNNTs are the object of intense experimental5
10 and theoretical11
15 analysis, associated with their possible applications in fields such as supertough composite materials, biosensors, and nanoelectronic devices.16 A SWBNNT (Figure 1) is formed by rolling an h-BN monolayer into a cylinder along a (n1,n2) lattice vector (see animations r 2011 American Chemical Society

at the CRYSTAL Web site).17 The n1 and n2 indices determine the diameter and chirality, which are the key parameters that characterize the structure.16
19 In recent years, much theoretical work has been devoted to the investigation of the response of nanotubes to the perturbation created by an electric field.11
15,20 In all these works, only the polarizability has been considered. In the present paper we show how the high-frequency (i.e., electronic) static polarizability (R) and first (β) and second (γ) hyperpolarizability tensors of zigzag (n,0)-SWBNNTs can be computed using a development version of the CRYSTAL code.21 The scheme is self-consistent (“local effects” are taken into account) and very efficient. Tubes as large as (50,0), containing 200 atoms in the unit cell, can be investigated with large basis sets using (in principle) any density functional theory (DFT) approximation. This permits us to study the dependence of the optical properties Special Issue: Richard F. W. Bader Festschrift Received: April 7, 2011 Revised: May 27, 2011 Published: June 23, 2011 12631

dx.doi.org/10.1021/jp203237m | J. Phys. Chem. A 2011, 115, 12631–12637

The Journal of Physical Chemistry A

ARTICLE

on the nanotube radius with reference to the limiting case of a BN monolayer. The paper is organized as follows. In section 2, the most relevant formulas are provided, and in section 3, the computational conditions are defined. In section 4 our results for the static linear and nonlinear polarizabilities of SWBNNTs are discussed with reference to the choice of the Hamiltonian, including the local density (LDA), generalized gradient (GGA) and Hartree
Fock (HF) approximations, as well as the “hybrid” functionals, such as B3LYP and PBE0. Finally, in section 5, some conclusions are drawn.

2. METHOD In our previous papers,22
26 we have shown how to carry out local basis set crystal orbital (CO) calculations of the static polarizability tensor (R), as well as first and second hyperpolarizability tensors (β and γ), for a periodic system based on a fully coupled perturbation theory formulation. The original coupledperturbed Hartree
Fock (CPHF) treatment was extended in ref 25 to the coupled perturbed Kohn
Sham (CPKS) R and the 2n + 1 formulation of β given by (see below eq 3 for the definition of the quantities appearing in these equations) ! BZ 2 occ virt ðtÞ ðuÞ ð1Þ Rtu ¼
Ptu R ∑ ∑ ∑ Ξjl Ulj k nk j l

(uv) are the first and second derivatives of the i,j • G(u) ij and Gij KS Hamiltonian matrix element with respect to the components of the applied electric field indicated by u and v. G(u) ij consists of two terms, namely: Ξ(u) ij , the first derivative of the twoelectron interaction i,j matrix element with respect to the u component of the applied electric field, and the analogous ). The latter is derivative of the DFT XC potential (G(u)XC ij calculated in the atomic orbital (AO) basis (see eq 19 in ref 25) and transformed to the CO basis. consists of the corresponding second derivatives except G(uv) ij for the dipole moment term which vanishes. In the AO basis (denoted by jμ, jν), the second derivative of the XC term in the KS matrix is (" ∂3 f XC t u XC gtu FF ½F μν ¼ wi ∂F3 i

∑

∂3 f XC ∂3 f XC t u u t ðF ∇F ∇F + F ∇F ∇F Þ + 4 ∇F 3 3 ∂F2 ∂j∇Fj2 ∂F∂ðj∇Fj2 Þ2 2 XC ∂2 f XC t u ∂ f tu t u 3 ∇F ∇F 3 ∇F + ∂F2 F + 2 ∂F∂j∇Fj2 ð∇F 3 ∇F # " ∂3 f XC tu 0 g + ∇F 3 ∇F Þ jμ jν + 2 Ft Fu ∇F ∂F2 ∂j∇Fj2 +2

∂3 f XC ðFt ∇F 3 ∇Fu + Fu ∇F 3 ∇Ft Þ∇F ∂F∂ðj∇Fj2 Þ2 ∂3 f XC +4 ∇F 3 ∇Ft ∇F 3 ∇Fu ∇F ∂ðj∇Fj2 Þ3 ∂2 f XC + ðFtu ∇F + Ft ∇Fu + Fu ∇Ft Þ ∂F∂j∇Fj2 ∂2 f XC +2 ð∇F 3 ∇Ftu ∇F + ∇Ft 3 ∇Fu ∇F + ∇F ∂ðj∇Fj2 Þ2 )  ∂f XC t u u t tu 0 g 3 ∇F ∇F + ∇F 3 ∇F ∇F Þ + ∂j∇Fj2 ∇F 3 ∇ðjμ jν Þ r +2

βtuv ¼

8 2 0 < BZ 2 occ virt virt ðtÞ @ ðuÞ ðvÞ 4
Ptuv R Ulj Glm Umj
: k nk j l m

∑ ∑∑

∑

occ

∑m

ðvÞ ðuÞ Ulm Gmj

ðvÞ

+i

∂Ulj

∂ku

13 A5 + βXC tuv

9 = ;

ð2Þ The CPKS procedure is more difficult to implement than CPHF because of the need to evaluate terms involving up to third-order derivatives with respect to F and/or |3F|2 in the case of a GGA exchange
correlation (XC) functional. Here we extend the CPKS treatment one step further to calculate γ using the 2n + 1 formulation:24 γtuvw ¼
Ptuvw

(

BZ

1 occ R nk j

∑k ∑

"

virt

all

∑

i

where fXC denotes an XC functional, g is a lattice vector and Ft, Ftu, 3Ft, and 3Ftu are the first and second derivatives of the electron density and density gradient, with respect to the is field components, evaluated at the grid point ri. G(uv)XC ij obtained from this matrix upon transformation to the CO basis. is an element of the block (occupied,virtual) off• U(u) ij diagonal matrix that determines the first-order CO coefficients as a linear combination of AOs:

occ

ðvwÞ ðvwÞ ðuÞ ∑l UljðtÞ ∑m GðuÞ lm Umj
∑ Ulm Gmj m

1 virt ðuvÞ ðwÞ 1 occ ðwÞ ðuvÞ + G U
U E 2 m lm mj 2 m lm mj

∑



) ðtÞ # ðvwÞ ∂Ulj Ulj + γXC tuvw ∂ku l

virt

+i

ð4Þ

∑

ð3Þ

Cμj ¼

∑i C0μi UijðuÞ

ðuÞ

Gij

ðuÞ

All matrices in eqs 1
3 are in the unperturbed CO basis; they depend upon the reciprocal space k point but, for brevity, this dependence has been omitted. In the above relations: • Ptu, Ptuv, and Ptuvw indicate the sum of all permutations of the subscripted field components. • nk denotes the number of k points in the first Brillouin zone (BZ). • R is the real part of the expression in parentheses that follows. • Ξ(t) ij is the t component of the appropriate dipole moment operator for periodic systems.

ð5Þ

It is given by Uij

ðuÞ

¼

½0

½0

Ej
Ei

ð6Þ

where E[0] and E[0] are the eigenvalues associated with the i j (occupied) ith and (virtual) jth unperturbed COs, respectively. U(u) is the result of a self-consistent first-order CPKS cycle (since G depends on U). 12632

dx.doi.org/10.1021/jp203237m |J. Phys. Chem. A 2011, 115, 12631–12637

The Journal of Physical Chemistry A

ARTICLE

provides arbitrary-order derivatives of a variety of functionals including GGAs, hybrids, and meta-GGAs.

Figure 1. Single-walled zigzag (12,0) BN nanotube.

• U(uv) is an element of the matrix that determines the secondij order CO coefficients as a linear combination of AOs: ðuvÞ

Cμj

¼

∑i C0μi UijðuvÞ

ð7Þ

The calculation of the block off-diagonal elements of U(uv) ij (e.g., virtual-occupied) implies a second-order self-consisthrough tent CPKS cycle as they are related to G(uv) ij ðtuÞ

ðtuÞ

Uij

Gij + Ptu ð ¼

virt

occ

∑l GðtÞil UljðuÞ
∑m UimðuÞ GðtÞmj Þ ½0

½0

Ej
Ei

ð8Þ

whereas the block-diagonal elements (e.g., occupied
occupied) are obtained from a product of first-order U matrices: ðtuÞ

Uij

1 virt ðtÞ ðuÞ ¼
Ptu Umi Umj 2 m

∑

ð9Þ

(v) • (∂U(v) lj )/(∂ku) is the derivative of Ulj with respect to the u component of the reciprocal lattice vector k, which is computed analytically.27,28 XC • βXC tuv and γtuvw contain third- and fourth-order derivatives, respectively, of fXC with respect to the applied field compoXC nents. βXC tuv is reported in eq 20 of ref 25; γtuvw is obtained by with respect to another field further differentiation of βXC tuv component. The derivatives of the XC functional fXC appearing in eq 4 are provided on the fly by the XCFun library,29,30 and we will therefore not consider them in any detail here. The XCFun library29,30 has been interfaced to a development version of the CRYSTAL program using (forward mode) automatic differentiation and

3. COMPUTATIONAL DETAILS The present CPKS calculations of electronic (hyper)polarizabilities were performed with a development version of the periodic ab initio CRYSTAL09 code21 using the HF approximation and several XC functionals in conjunction with a split valence 6-31G* basis set, which has been shown to provide accurate polarizabilities.31 The XC contributions were evaluated by numerical integration over the unit cell volume. In CRYSTAL, radial and angular points of the integration grid are generated through Gauss
Legendre radial quadrature and Lebedev twodimensional angular point distributions. In our calculations, a (75,974) pruned grid (XLGRID keyword in the CRYSTAL09 manual),32 corresponding to 75 radial and 974 angular points, was employed. The integration accuracy can be estimated by the error in the electronic charge of the unit cell (Δe = 2.8  10
3|e| out a total of 288 electrons for the (12,0) BNNT). Other details on the grid generation and its influence on the accuracy and cost can be found in refs 33, 34, and 35. Evaluation of the Coulomb and exact exchange infinite series is controlled by five parameters32 (T1, T2, T3, T4, T5), whose values are set to T1 = T2 = T3 = T4 = 1/2T5 = 8. This setting corresponds to an uncertainty in the determination of the (hyper)polarizability tensor components below 0.4%. The selfconsistent-field calculation of the wave function for the unperturbed systems and the CPKS cycles for the various tensor components were all well converged, i.e., to 10
10 for the former and 10
4 au for the latter. Convergence of the CPKS cycles benefited from mixing Fock matrix derivatives of successive cycles to damp oscillatory behavior. A mixing factor of 50% (FMIXING parameter)32 turned out to be effective. The unit cell of the BN nanotubes (or a monolayer) may be considered small because the lattice parameter along the tube axis is short and is independent of the tube diameter. For small unit cell compounds, it has been shown25,26 that the number of reciprocal points (nk in eqs 2 and 3) must be relatively large to obtain accurate values of β and γ. The dependence of the nonnull components of the (hyper)polarizabilities on the shrinking factor S has been examined using the B3LYP Hamiltonian for the (6,0)-BN nanotube and the h-BN monolayer (all other nanotubes considered here may be thought of as lying between these two limiting cases). It was found that the calculation of the (hyper)polarizability tensors for the monolayer requires a higher value of S than for the SWBNNT to provide accurate results, in particular in the case of components involving the periodic directions x and y. However, the asymptotic values are reached within six significant figures in all cases for S g 24. With reference to the asymptotic limits, S = 12 leads to an error that is far below 1% in γxxxx or γxxyy and is reduced to 0.014% with S = 16, whereas the error in βxxx and the “nonperiodic” components of γ are smaller by 2 orders of magnitudes. Thus, S = 16 was used for all calculations because it provides very accurate results for all tensor components through fourth order. 4. RESULTS AND DISCUSSION A flat h-BN sheet corrugates slightly upon rolling and forms a thin double layer with B atoms in the inner layer, i.e., nearer the tube axis (Figure 1). With B3LYP the separation between the two layers is 0.2 Å for the (3,0)-SWBNNT and decreases rapidly with 12633

dx.doi.org/10.1021/jp203237m |J. Phys. Chem. A 2011, 115, 12631–12637

The Journal of Physical Chemistry A

ARTICLE

Table 1. Static (Hyper)Polarizabilities per BN Unit of the Zigzag (12,0)-BN Nanotube and the h-BN Monolayer with Different Hamiltonians (au) structure

Hamiltonian

(12,0)-SWBNNT

h-BN monolayer

Rxx

Rzz

βxxx

βxzz

γxxxx

γxxzz

γzzzz

SVWN

36.23

11.93

501.2


50.94

109506

3073

2002

PBE B3LYP

36.00 31.34

12.03 11.34

494.5 358.7


51.01
42.91

102038 48124

2982 1668

1987 1337

1483

1234

PBE0

30.52

11.18

327.3


40.61

40778

HF

22.32

9.76

143.1


24.21

9014

568.1

SVWN

35.30

5.566

460.8

99994

408.9

68.02

PBE

35.65

5.634

460.2

95687

429.0

79.26

B3LYP

30.98

5.504

341.3

44784

366.5

67.61

PBE0

29.97

5.475

312.4

38023

357.0

76.03

HF

21.90

5.254

143.8

8840

242.4

61.82

increasing nanotube radius to 0.1 Å for n = 6 and 0.04 Å for n = 12. Before calculating the dielectric properties, equilibrium structures of the SWBNNT have been determined for every value of n considered by minimizing the total energy through a conjugate-gradient algorithm, as implemented in CRYSTAL.32 Because the large radius behavior as well as the trend of R, β, and γ toward the monolayer limit are of interest here, nanotubes of different sizes were investigated. We adopt the following orientation for the monolayer and the nanotubes: x and y are periodic directions in the monolayer plane, whereas SWBNNTs are periodic along the x direction. With this orientation, the following inequivalent components of the static (hyper)polarizability tensors are nonzero for the h-BN monolayer: Rxx, Rzz, βxxx, βxyy, γxxxx, γxxyy, γxxzz, and γzzzz. The nonzero inequivalent components for the SWBNNTs are the same except for the replacement of γxxyy by γyyzz. Even though the lists of nonvanishing components are very similar for the nanotube and the monolayer, the correspondence between components in the very large radius limit is not straightforward. This is because the equivalent transverse directions in the monolayer (x and y) are periodic whereas the perpendicular direction is nonperiodic; in contrast, for the SWBNNT the equivalent y and z directions are nonperiodic and there is periodicity only in the x direction. However, several relations between the components of the two structures, in the limit of an infinite nanotube radius, can be established. For R, as already noticed by Guo et al.,12 through circular averaging one can obtain (the superscipts n and l refer to the nanotube and monolayer, respectively) 0 1 l n l B Rzz ¼ 2 ðRzz + Ryy Þ ð10Þ @ Rnxx ¼ Rlxx Again, by circular averaging, the following relations may be derived for β and γ: 0 1 l n B βxyy ¼ 2 βxyy ð11Þ @ n βxxx ¼ βlxxx 0

3 γnzzzz ¼ ðγlzzzz + γlyyyy + 2γlyyzz Þ B 8 B B γn ¼ γl xxxx B xxxx @ 1 l n γxxzz ¼ ðγxxyy + γlxxzz Þ 2

ð12Þ

590.4

We remind that Rlyy in eq 10 is equal to Rlxx and that in eq 12 γlyyyy = γlxxxx, γlyyzz = γlxxzz, and γnxxzz = γnxxyy. There are other relations that hold separately for the nanotube and the monolayer due to symmetry, and for this reason results for βlxyy, γnyyzz, and γlxxyy are not reported because βlxyy =
βlxxx, γnyyzz = 1/3γnzzzz, and γlxxyy = 1/3γlxxxx. The effect of the choice of functional on the determination of the (hyper)polarizability tensors for one representative nanotube, the (12,0)-SWBNNT, and the h-BN monolayer is shown in Table 1 (all values in au per BN unit, which are equivalent to 1.649  10
41 (C2 m2)/J, 3.206  10
53 (C3 m3)/J2, and 6.235  10
65 (C4 m4)/J3 for R, β and γ, respectively). The following functionals have been examined: the SVWN representation of LDA; PBE as an example of a GGA; the hybrids B3LYP and PBE0; and for comparison purposes also HF. As already observed for three-dimensional systems,25,26 LDA and GGA provide very similar estimates and are generally believed to overshoot the (hyper)polarizabilities, whereas it is also wellknown that HF can undershoot them significantly. As expected, the two hybrid functionals provide similar results that lie within the range spanned by the LDA and HF values, with the PBE0 results being slightly closer to HF than B3LYP as a consequence of the larger amount of HF exchange incorporated in the former (25%, as compared to 20% in the latter). Our calculated values vary significantly with the choice of functional, particularly for the γ tensor. Unfortunately, due to difficulties in the purification of synthesized SWBNNT samples,16 no experimental measurements of the optical properties are available. However, it is well-known that the (hyper)polarizabilities are strongly influenced by the band gap. Arenal et al.,5 Jaffrennou et al.,6,7 and Lee et al.8,9 have reported an experimental band gap of 5.8
5.9 eV. The band gap predicted by hybrid functionals for the (12,0) nanotube is 5.9 and 6.3 eV with B3LYP and PBE0, respectively, whereas it is much lower (4.1 eV) according to LDA (the PBE band gap of 4.2 eV is quite similar) and much larger (13.2 eV) in the HF approximation. Even considering that real samples contain much larger radius nanotubes than the (12,0) case and that the band gap calculated with B3LYP for very large SWBNNTs reaches about 6.3 eV, we expect the hybrid functionals to perform the best as far as optical properties are concerned. Although B3LYP calculations result in a large overshoot of longitudinal hyperpolarizabilities for small-gap conjugated polymers,36 this problem does not exist for wide band gap systems such as BN nanotubes. Thus, the B3LYP functional is here used to analyze the linear and nonlinear optical properties of SWBNNTs and to verify their 12634

dx.doi.org/10.1021/jp203237m |J. Phys. Chem. A 2011, 115, 12631–12637

The Journal of Physical Chemistry A

ARTICLE

Figure 3. Rzz and γzzzz per BN unit as a function of the (n,0)-BN nanotube size. Due to the slow convergence of γzzzz data have been fitted to the function c0 + c1/n + c2/n2 + c3/n3. Rzz in au and γzzzz in 1000 au.

Figure 2. Rxx, βxxx, and γxxxx per BN unit as a function of the (n,0)-BN nanotube size. The parameters in the inset are obtained by fitting values for the nine largest-radius nanotubes (n = 20
50) to the polynomial function c0 + c1/n + c2/n2. Asymptotes represent the h-BN monolayer limit. Values in atomic units (10 000 au for γxxxx).

trend toward the asymptotic limit given by the h-BN monolayer values. To study the variation of the (hyper)polarizabilities as a function of radius, we considered SWBNNT zigzag nanotubes in the range (3,0) to (50,0). Our results are shown in Figures 2
5. It may be seen that some of the curves exhibit a maximum or minimum at a small radius in the range n = 3
6. In

particular, the (3,0)-SWBNNT possesses a quasi triangular section with interatomic distances similar to those in the bulk and from the numerical values there appears to be a discontinuity between n = 3 and n = 4. Monotonic behavior is observed for R within the entire range of n values, as well as for β and γ when n g 6. For n g 6 all curves are in addition in the asymptotic region, i.e., beyond the point of inflection. For a given tensor the components differ significantly in magnitude from one another, with the largest values corresponding to the direction of periodicity. Although the components involving nonperiodic directions would likely increase significantly by enriching the basis set for those directions,37 this general behavior is expected to remain. All components of R, β, and γ within the (6,0)
(50,0) range tend asymptotically to a limiting value in a smooth and regular way. However, while the pure periodic components (Figure 2) converge rapidly, all the other components converge much more slowly. To assess the limiting values, extrapolations were carried out using a truncated power series in 1/n, as indicated in the figure captions. For the slowly convergent cases (Figures 3
5) a relatively small set of large n nanotubes was used to obtain a reasonable fit. For those tensor components involving the periodicity direction only, the extrapolated values are Rxx = 30.95, βxxx = 341, and γxxxx = 4.475  104 au, respectively, which closely approximate 12635

dx.doi.org/10.1021/jp203237m |J. Phys. Chem. A 2011, 115, 12631–12637

The Journal of Physical Chemistry A

ARTICLE

pure nonperiodic component, γzzzz, is more problematic. In this case the value we obtain by extrapolation undershoots the correct layer limit of 1.709  104 au obtained from eq 12 by about 13%.

Figure 4. βxyy per BN unit as a function of the nanotube size n. Other details as in Figure 3.

5. CONCLUSIONS The coupled perturbed Kohn
Sham (CPKS) scheme for calculating linear and nonlinear electric susceptibilities recently implemented in the CRYSTAL code has been extended to thirdorder and used to obtain static (hyper)polarizability tensors of zigzag SWBNNTs. It is seen that values determined for LDA/ GGA functionals are very different from those obtained from hybrid functionals. Using the B3LYP hybrid, which yields a good value for the bandgap, we have carried out calculations on (n,0) nanotubes with n up to 50, corresponding to radii as large as 20 Å. For the larger radius tubes, the pure longitudinal component of each tensor has a magnitude that is 2
6 times larger than any other component. This value saturates rapidly with increasing tube radius. It is within 1% of the infinite radius limit for n = 17, 22, and 42, for R, β, and γ respectively. Extrapolation gives values that are virtually identical, as they should be, to those are predicted on the basis of h-BN layer calculations. All other nonvanishing components saturate more slowly. For the largest nanotubes examined here (n = 50) the fraction of the infinite radius limiting value that is obtained ranges from a high of 84% for Rzz to a low of 63% for γzzzz. Nonetheless, extrapolation leads to results that lie within 4% of the “exact” h-BN layer value except for γzzzz where the discrepancy is 13%. Considering the very slow convergence in this case, we regard that level of agreement as satisfactory. ’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

Figure 5. γxxzz per BN unit as a function of the nanotube size n. Other details as in Figure 3.

the BN monolayer values (Rxx = 30.98, βxxx = 341.3, and γxxxx = 4.479  104 au). As expected, the convergence is most rapid for Rxx and slowest for γ. To achieve a value within 1% of the asymptotic limit, the size of the nanotube required is n = 17, 22, and 42 for R, β, and γ. All the remaining nonzero components of the hyperpolarizability tensors involve at least two electric field components along nonperiodic directions. The limiting values are approached from below in all cases and convergence is slow. γzzzz, which does not involve any component of the electric field along the periodicity direction, also shows a point of inflection for n e 20. Because of this behavior the fitting polynomial has been extended to the third power in these cases, and the best-fitting procedure has been applied only to the subset of data on the right of the flexus. The asymptotic values so determined are as follows: Rzz = 17.8, βxyy =
164, γxxzz = 737  10, and γzzzz = 1.49  104 au. Despite some difficulty in fitting, the extrapolated values are close to the h-BN monolayer limit given by eqs 10
12. In particular, Rzz undershoots the value of 18.24 au, based on layer calculations, by about 2% and γxxzz undershoots the layer predicted value of 7647 au by about 3.6%. The fitting of the

’ ACKNOWLEDGMENT The authors thank CINECA for granting access to supercomputing resources (ISCRA Award No. HP10AC4ZGA.2010). R.B., K.R. and U.E. have been supported by the Research Council of Norway through Research Grant No. 179568/V30. ’ REFERENCES (1) Rubio, A.; Corkill, J.; Cohen, M. Phys. Rev. B 1994, 49, 5081. (2) Blase, X.; Rubio, A.; Louie, S.; Cohen, M. Europhys. Lett. 1994, 28, 335. (3) Chopra, N.; Luyken, R.; Crespi, V.; Cohen, M.; Louie, S.; Zettl, A. Science 1995, 269, 966. (4) Chen, Y.; C., S. J.; Zou, J.; Caer, G. L. Appl. Phys. Lett. 2004, 84, 2430. (5) Arenal, R.; Stephan, O.; Kociak, M.; Taverna, D.; Loiseau, A.; Colliex, C. Phys. Rev. Lett. 2005, 95, 127601. (6) Jaffrennou, P.; Donatini, F.; Barjon, J.; Lauret, J. S.; Maguer, A.; Attal-Tretout, B.; Ducastelle, F.; Loiseau, A. Chem. Phys. Lett. 2007, 442, 372. (7) Jaffrennou, P.; Barjon, J.; Lauret, J. S.; Maguer, A.; Golberg, D.; Attal-Tretout, B.; Ducastelle, F.; Loiseau, A. Phys. Stat. Sol. b 2007, 244, 4147. (8) Lee, C. H.; Wang, J.; Kayashta, V. K.; Huang, J. Y.; Yap, Y. K. Nanotechnology 2008, 19, 455605. (9) Lee, C. H.; Xie, M.; Kayashta, V.; Wang, J.; Yap, Y. K. Patterned Growth of Boron Nitride Nanotubes by Catalytic Chemical Vapor Deposition. Chem. Mater. 2010, 22, 1782. 12636

dx.doi.org/10.1021/jp203237m |J. Phys. Chem. A 2011, 115, 12631–12637

The Journal of Physical Chemistry A

ARTICLE

(10) Sun, C.; Yu, H.; Xu, L.; Ma, Q.; Qian, Y. J. Nanomater. 2010, 2010, 163561. (11) Hao, S.; Zhou, G.; Duan, W.; Wu, J.; Gu, B. L. J. Am. Chem. Soc. 2006, 128, 8453. (12) Guo, G. Y.; Ishibashi, S.; Tamura, T.; Terakura, K. Phys. Rev. B 2007, 75, 245403. (13) Guo, G. Y.; Lin, J. C. Phys. Rev. B 2005, 71, 165402. (14) Lan, H.-P.; Ye, L.-H.; Zhang, S.; Peng, L.-M. Appl. Phys. Lett. 2009, 94, 183110. (15) Margulis, V. A.; Muryumin, E. E.; Gaiduk, E. A. Phys. Rev. B 2010, 82, 235426. (16) Golberg, D.; Bando, Y.; Huang, Y.; Terao, T.; Mitome, M.; Tang, C.; Zhi, C. ACS Nano 2010, 4, 2979. (17) http://www.theochem.unito.it/crystal_tuto/mssc2008_cd/ tutorials/nanotube/nanotube.html. (18) Dai, H. Acc. Chem. Res. 2002, 35, 1035. (19) White, C. T.; Robertson, D. H.; Mintmire, J. W. Phys. Rev. B 1993, 47, 5485. (20) Wang, L.; Lu, J.; Lai, L.; Song, W.; Ni, M.; Gao, Z.; Mei, W. N. J. Phys. Chem. C 2007, 111, 3285. (21) Dovesi, R.; Orlando, R.; Civalleri, B.; Roetti, C.; Saunders, V. R.; Zicovich-Wilson, C. M. Z. Kristallogr. 2005, 220, 571. (22) Ferrero, M.; Rerat, M.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2008, 29, 1450. (23) Ferrero, M.; Rerat, M.; Orlando, R.; Dovesi, R. J. Chem. Phys. 2008, 128, 014100. (24) Ferrero, M.; Rerat, M.; Kirtman, B.; Dovesi, R. J. Chem. Phys. 2008, 129, 244110. (25) Orlando, R.; Lacivita, V.; Bast, R.; Ruud, K. J. Chem. Phys. 2010, 132, 244106. (26) Orlando, R.; Ferrero, M.; Rerat, M.; Kirtman, B.; Dovesi, R. J. Chem. Phys. 2009, 131, 184105. (27) Kirtman, B.; Gu, F. L.; Bishop, D. M. J. Chem. Phys. 2000, 113, 1294. (28) Bishop, D. M.; Gu, F. L.; Kirtman, B. J. Chem. Phys. 2001, 114, 7633. (29) Ekstr€om, U. XCFun library, http://www.admol.org/xcfun, 2010 (30) Ekstr€om, U.; Visscher, L.; Bast, R.; Thorvaldsen, A. J.; Ruud, K. arbitrary-order density functional response theory from automatic differentiation. J. Chem. Theory Comput. 2010, 6, 1971–1980. (31) Ferrabone, M.; Kirtman, B.; Rerat, M.; Orlando, R.; Dovesi, R. Phys. Rev. B 2011, 83, 235421. (32) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M. CRYSTAL 2009 User’s Manual; University of Torino: Torino, 2009. (33) Pascale, F.; Zicovich-Wilson, C. M.; Orlando, R.; Roetti, C.; Ugliengo, P.; Dovesi, R. J. Phys. Chem. B 2005, 109, 6146. (34) Prencipe, M.; Pascale, F.; Zicovich-Wilson, C.; Saunders, V.; Orlando, R.; Dovesi, R. Phys. Chem. Min. 2004, 31, 559. (35) Tosoni, S.; Pascale, F.; Ugliengo, P.; Orlando, R.; Saunders, V. R.; Dovesi, R. Mol. Phys. 2005, 103, 2549. (36) Champagne, B.; Perpete, E. A.; van Gisbergen, S. J. A.; Baerends, E. J.; Snijders, J. G.; Soubra-Ghaoui, C.; Robins, K.; Kirtman, B. J. Chem. Phys. 1998, 109, 10489. (37) Champagne, B.; Bulat, F. A.; Yang, W.; Bonness, S.; Kirtman, B. J. Chem. Phys. 2006, 125, 194114.

12637

dx.doi.org/10.1021/jp203237m |J. Phys. Chem. A 2011, 115, 12631–12637