2545
J . Phys. Chem. 1994,98, 2545-2554
Density Functional Theory Predictions of Polarizabilities and First- and Second-Order Hyperpolarizabilities for Molecular Systems Nobuyuki Matsuzawa* SONY Corporation Research Center, 174 Fujitsuka-cho, Hodogaya-ku, Yokohama 240, Japan
David A. Dixod DuPont Central Research and Development, Experimental Station,? P.O. Box 80328 Wilmington. Delaware 19880-0328 Received: September 1 , 1993'
First- (8) and second-order (y) hyperpolarizabilities of substituted benzenes and polyene oligomers are calculated by applying density functional theory (DFT) at the local (LDFT) and nonlocal (NLDFT) levels with field induced polarization (FIP) functions augmenting the basis sets. It is found that the LDFT method without the FIP functions gives 8 and y values which are too small as compared to experimental ones, whereas inclusion of the FIP functions at the LDFT or NLDFT level leads to an increase in the 8 and y values giving good agreement with the experimental values. It is also found that inclusion of nonlocal corrections improves the accuracy of the calculated 8 and y values although the improvement is not large. From comparison with a b initio molecular orbital theory values obtained by others, it is found that the NLDFT/FIP calculations yield values which are comparable to those at the MP-2 level.
Introduction The prediction of molecular propertiesof materials is becoming more important as the cost to design such materials continues to rise. A variety of theoretical methods have been used to predict molecular structure, and from these structures, the associated wavefunctions and properties can then be calculated. The most accurate predictions have, in general, come from ab initio molecular orbital theory. However, ab initio molecular orbital theory does not scale well with the size of the molecular system. For N basis functions, the Hartree-Fock method scales formally as N4, and if correlation effects are included, the scaling goes as Nm, m 1 5 . Although many molecular chemical problems can be solved at the Hartree-Fock level, the inclusion of correlation effects is often required for quantitative predictions. There is considerable interest in exploiting new computational methods that scale with a smaller value of m. The method that we have been testing for chemical systems along with other research groups is thatofdensityfunctional theory (DFT) in both thelocal ( L D m ) and nonlocal (NLDFT) approximations. This method has been extensively used to predict the propertiesof solids but only recently has it been used to predict the properties of molecular systems.' We have had a significant amount of success in using LDFT for the prediction of the structuresand vibrational spectraof molecular systems that require a high-level correlation energy treatment in order to obtain even qualitative agreement with experiment.Q.2 There is real interest in discovering new materials for use in nonlinear optical applications such as optical switches and high density optical recording.3 We are interested in predicting the nonlinear optical properties of molecules by using theoreticall computational methods. Our initial studies focused on using a semiempirical method based on the MNDO Hamiltonian to predict the NLO properties of a wide range of conjugated molecules for direct comparison to experiment.' A semiempirical method was employed because we wanted to validate its use as a survey tool for large molecules, as traditional ab initio molecular orbital calculationsincorporatingcorrelationeffects and adequate basis sets are simply too expensive computationally. In line with our other studies of the applications of DFT to chemical t Contribution No. 6638.
Abstract published in Advance ACS Absrracrs. December 15, 1993.
QQ22-3654/94/2098-2545%04.5Q/Q
systems,1r.g.2we decided to investigate its use for the prediction of NLO properties ( p , a!, 8, and y) for some model molecular systems. The LDFT method has been used to calculatethe secondorder hyperpolarizability (y) of benzene and (260.5 Consistent with calculations of other molecular properties, LDFT theory was found to give reasonable results for the prediction of a! and y for these two compounds. In other previous work, Jasien and Fitzgerald6 reported a study of molecular polarizabilities at the LDFT level. As we were completing this study, Chong and coworkers published their LDFT calculations for b for some simple molecules where they included field induced polarizationfunctions (FIP) (diffuse functions) in the basis set.' The need for diffuse functions for calculating such properties has been discussed extensively in terms of traditional ab initio molecular orbital approaches.* Our original LDFT study without diffuse functions has been extended to include such functions and the effect of nonlocal corrections was also examined.
Calculations The calculations described below were done with the program system DMol.9 The atomic basis functionsare given numerically on an atom-centered, spherical-polar mesh. The radial portion of the grid is obtained from the solution of the atomic DFT equations by numerical methods. Since the basis sets are numerical, the various integrals arising from the expression for the energy need to be evaluated over a grid in terms of radial functions and spherical harmonics. The number of radial points NR is given as
NR= AR,
1.4(2
+ 2)'13
(1) where Z is the atomic number and R ,, is the maximum distance for any function given in au. The angular integration points Ne are generated at the NRradial points to form shells around each nucleus. The value of Ne ranges from 14 to N- depending on the behavior of thedensity and a maximum 1value for the spherical harmonic, These quantities for the two meshes used in this study are given in Table 1. The Coulomb potential corresponding to the electron repulsion term is determined directly from the electron density by solving Poisson's equation. The form for the exchange-correlation energy of the uniform electron gas for the 0 1994 American Chemical Society
2546 The Journal of Physical Chemistry, Vol. 98,No. 10, 1994
Matsuzawa and Dixon
TABLE 1: Mesh Parameters parameter
FINEu
~
A R*X
1.2 12.0
thresholdb
L* NaU a
0.00001 29 302
XFINEa 1.5 15.0 0.000001 35 434
Mesh parameter in DMol. Angular sampling threshold.
LDFTcalculationsis that derived byvon Barth and Hedin (BH)loa or that derived by Volk, Wilks, and Nusair (VWN).Iob Calculations including nonlocal corrections (NLDFT) were also performed, The gradient corrections to the exchange potential derived by Lee,Yang,andParr (LYP)I1and that to thecorrelation potential derived by Beckelf were used (BLYP). The LYP correction also requires the use of a different local potential as described in ref 11. The DMol calculations were done starting with a double numerical basis set augmented by polarization functions (DNP). Calculations were also performed with a basis set obtained by augmenting the DNP basis set with the FIP basis functionsgiven by Chong and co-~orkers.~ Calculations were performed with the FIP d functions for first row atoms and the FIP p function on H (DNP d p). As recommended by Chong and coworkers, the exponents of the tight d polarization were changed to their recommended values. For some cases, FIP s and p functions were also added to the first row atoms (DNP + spd + p). The former augmented basis set is denoted as DNP+ 1, and the latter as DNP+2. The multipolar fitting functions for the model density used to fit the effective potential have angular momentum numbers, I, one greater than that of the polarization functions. Geometrieswere optimizedby using analytic gradient methods at the LDFT level with the von Barth and Hedin potential and the DNP basis ~ e t . * + ~The J 3 density was converged to 1od for the geometry optimizations, and the FINE mesh was used. The calculation of the nonlinear optical properties was done with a finite field approach.3 In this approach, the response of the ground state charge distribution to an applied electric field at zero frequency is investigated. A molecule in an applied electric field will exhibit an induced dipole moment. This induced moment can be expanded in a Taylor series in powers of the applied electris field as shown in eq 2
+ +
X
H
\ /c=c,
H
-
P. CI. Br. I
/H H
Figure 1. Orientation of the calculated molecules and the numbering systems of the atoms (amino hydrogens are above the xy-plane). as follows (4)
where g1 in eq 5 is defined as
BI = A similar expression is available for the energy. A dipole expansion was used as opposed to an energy expansion because the field occurs as the third power in the expression for y based on the dipole expansion and the fourth power in the energy expansion.s Part of this work is to study the dependence of the NLO properties on various computational parameters. In order to use the perturbation expansion, the value of the applied field must be small enough so that the expansion given in eq 2 is valid but large enough to incorporate higher order nonlincarities; thus, values for the field ranging from 0.001 to 0.01 au were used. Because of the small values of the field, the evaluation of y requires 2 to 3 orders of magnitude more accuracy in the energy expansion as compared to the dipole expansion. The initial calculations employed the equations given by Kurtz et al.I4 for p, a,8, and y in terms of dipole moments calculated in an applied electric field. Sim et al.lS have shown that there is a small error in the equation for 8 , due ~ to the assumption made by Kurtz et al. that p determined from the finite field energy expansion is the same as evaluated from the wavefunction. Values from both expressions are given below and the value based on the expression of Sim et al. is the preferred one.15 The quantities of p, a,8, and y calculated from their vector or tensor components are defined
BIjj
(7)
I
In order to account for the presence of the field terms in the expansion for direct comparison to experiment, B and y need to be multiplied by the appropriate correction factors. If the experimental hyperpolarizabilitiesare defined as shown in eq 8
=
+ f f l f+ ] BIjpfk + yIjklFfp, + "*
(8) the experimental numbers were multiplied by factors of 2 and 6 for B and y,respectively. The electric field Fis sometimes defined as
+
F = F(dO' e4"") (9) and the experimental numbers for 8 and y need to be multiplied by factors of 1/2 and 1/4, respectively, for direct comparison to our calculated values. Results and Discussions Structures. The orientation of the calculated molecules with respect to the Cartesiancoordinates are shown in Figure 1,together with the numbering system of the atoms for each molecule. The optimized geometries are given in Table 2 and compared to the experimental values.16
Predictions of Polarizabilities The LDFT calculated value for r(C-C) for benzene is in good agreement with the experimental value, but for r(C-H), the calculated value is longer than the experimental value by 0.014 A. For nitrobenzene, the calculated C-C and C-H bond lengths are in good agreement with the measured ones as is r(N-O), whereas the calculated value for r(C-N) is too short by 0.036 A. The calculated angles are, in general, in good agreement with the largest difference (2.4O) being found for B(CIN0). For aniline, the calculated bond lengths are in good agreement with the experimental values except for r(N-H), where the calculated value is 0.027 A longer than the measured one. The next largest difference is for r(Ca-H3) where the difference is 0.018 A. For the bond angles, the largest difference is found on the amino group (B(HNH)) with the calculated value 2.9O smaller than the measured angle. The angle 7 (angle between the phenyl ring and NHN plane) is calculated to be 6.1O larger than the experimental value. For p-nitroaniline, a complete gas-phase experimental structure is not available, so a comparison to both the partial gas-phase structure and to the crystal structure was made. The calculated value for r(c1-N~)is too short by 0.039 A, similar to what was found for nitrobenzene. There is better agreement between theory and experiment for this value determined from the crystal data. The amino group is clearly nonplanar as found experimentally. The calculated value for r(C,-N1) is longer than the experimental values by 0.01-0.02 A. For fluorobenzene, the calculated C-C and C-F bond lengths are in good agreement with the measured ones with differences less than 0.006 A. Larger differencescan be found for the C-H bond lengths with the largest difference being 0.017 A for r(C4H4). The angles are in good agreement with differences less than 0.8O. For chlorobenzene, the same tendency exists, although the difference between the calculated and measured values is somewhat larger. The calculated C-C and C-CI bond lengths are in agreement with the experimental ones with the largest difference being 0.015 A for r(C&). The differences for the C-H bond lengths are less than 0.018 A. The angles are in good agreement with differences less than 0.9O. For bromobenzene, again, the calculated C-C and C-Br bond lengths are in good agreement with the measured ones with differences less than 0.008 A. The differences are slightly more pronounced for C-H bond lengths, although they are still less than 0.012 A. The calculated angles are in good agreement with the measured ones with differences less than 0.9O except for B(C&-HI), where a difference of 2.3O is found. For iodobenzene, the gas-phase parameters were obtained by assuming that the benzene ring is a regular hexagon. The calculated value for r(C-I) is 0.03 A longer than the experimental value and the calculated ring C-C bond lengths are shorter than the experimentalones. Considering the approximationsin the experimental structure, the agreement is quite good. For cyanobenzene, the differences in bond lengths between the calculated and measured ones are slightly larger as compared to the halobenzenes. The largest difference is for ~ ( C I - C ~with ) the calculated value being smaller by 0.033 A than the measured one. The calculated angles are in good agreement with the measured ones with the largest difference being 1.8O for B ( C r CI-CS). For ethylene, the agreement between theory and experiment is good with differences in bond lengths and angles less than 0.006 A and 0.6O, respectively. For tram-butadiene, the agreement in the bond lengths is not as good as found for the case of ethylene or for the substituted benzenes. The calculated r ( C r C3) single bond is shorter by 0.046 A than the measured one, whereas that for r(C1-C2) double bond is longer by 0.003 A. The calculated C-H bond lengths are all longer by 0.016-0.020 A as compared to the experimentalvalues. The calculated bond angles are in agreement with the measured angles with the largest difference being 1.4O. For tram-hexatriene, the agreement
The Journal of Physical Chemistry, Vol. 98, No.10, 1994 2547 between the calculated and measured bond lengths is improved as compared to tram-butadiene. The calculated r(C&) single bond is shorter by 0.029 A as compared to the measured one, whereas the r(CI-C2) and r ( C 4 4 ) double bonds are longer by 0.007 and 0.015 A, respectively. The calculated C-H bond len@hs are all longer by 0.001-0).007A as compared to the measured ones. The difference in the bond angles between the calculated and measured is somewhat larger as compared to ethylene and tram-butadiene. The largest differences are for ~ ( C I - C ~ H ~ ) , e(Cl-C&3), and O(C,-C3-C1) with differences of 2.4O, 3.7O, and 3.2O, respectively. Hyperpolarizrbilitie-s. The initial plan was to use the FINE mesh and the DNP basis set for evaluating thevarious properties. This mesh and basis set has been shown to give good predictions of geometries, frequencies, and reaction energies.2 The dependence of the various polarizabilities on the quality of the mesh, the level of convergence of the SCF quations, and the strength of the field were tested. The initial results for benzene showed that within a few percent, the value of y can be calculated with the FINE grid as compared to the use of a finer grid, if a field strength of 0.005 au and a convergence level of 1 V is used.’ The calculation of a was found to be independent of the parameters, if at least a FINE grid, a field strength of 20.001 au, and a convergencelevel of 10-7 are used. The initial studies on benzene gave a = 1.012 X 10-23 cm395 in good agreement with an experimental value of 1.032 X l e z 3cm3.I7 The value of y was calculated to be 6.01 X 10-36 esu,16 in good agreement with an experimental value of 5.9 X 10-36 esu.18 The computational parameters were further tested with calculations on nitrobenzene for which p and Bare nonzero with the results shown in Table 3. The values for p and a are converged at the levels shown in Table 3, although a slight increase in p can be seen for the results with a field strength of 0.010 au. The value of B is converged at an SCF convergenceof lP7. However, there is a significant difference in the value of B obtained from Kurtz et ale’sexpression and that from Sim et 81,’s expres.sion, if the applied field strength is 20.005 au. Thevalueof @obtained from Kurtz et al.’s expression begins to decrease in magnitude. The two expressions for calculating /3 are given in q s 10 and 1la14J5
B,,, = [-(2/5)P,(O) + (4/3)(P,(4) + P,(-F,)) ( l / w ( P , ( 2 q + r,(-24))I/F? (11) In these quations, Fdenotes the applied field, and i andf denote Cartesian axes. Following Kurtz et al., q 10 is used both for the diagonal and off-diagonal terms of 8, whereas following Sim et a1.h expression, q 10 is used for the off-diagonal terms and q 11 is used for the diagonal terms. As described above, p does depend on the field strength especially for higher field strengths, which can be attributed to the effect of higher order terms (P, I”,...) in the expansion. BecauseSim et al. use pdO) for calculating 8, the effect of the change in p due to higher fields is accounted for properly. Thus, the value of B obtained from q 10 and 11 was used below. For such ps, the quantity is converged if the applied field 20.005 au. The values of /3 are not strongly dependent on the quality of the grid for the FINE and XFINE grids. The value of y is approximately converged at 10-7 convergence for the SCF but not converged if the applied field
