Polarizabilities of Ï-Conjugated Chains Revisited - ACS Publications

Subscriber access provided by - Access paid by the | UCSB Libraries

Article

Polarizabilities of #-Conjugated Chains Revisited: Improved Results from Broken-Symmetry Range-Separated DFT and New CCSD(T) Benchmarks M. Belen Oviedo, Niranjan V. Ilawe, and Bryan M. Wong J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00360 • Publication Date (Web): 22 Jun 2016 Downloaded from http://pubs.acs.org on June 28, 2016

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 36

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

Polarizabilities of π-Conjugated Chains Revisited: Improved Results from Broken-Symmetry RangeSeparated DFT and New CCSD(T) Benchmarks M. Belén Oviedo, Niranjan V. Ilawe and Bryan M. Wong* Department of Chemical & Environmental Engineering and Materials Science & Engineering Program, University of California-Riverside, Riverside, California 92521,United States *Corresponding author. E-mail: [email protected]. Homepage: http://www.bmwong-group.com

ABSTRACT. We present a detailed analysis of non-empirically tuned range-separated functionals, with both short- and long-range exchange, for calculating the static linear polarizability and second hyperpolarizabilities of various polydiacetylene (PDA) and polybutatriene (PBT) oligomers. Contrary to previous work on these systems, we find that the inclusion of some amount of short-range exchange does improve the accuracy of the computed polarizabilities and second hyperpolarizabilities. Most importantly in contrast to prior studies on these oligomers, we find that the lowest-energy electronic states for PBT are not closed-shell singlets, and enhanced accuracy with range-separated DFT can be obtained by allowing the system to relax to a lower-energy broken-symmetry solution. Both the computed polarizabilities

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

Page 2 of 36

and second hyperpolarizabilities for PBT are significantly improved with these broken-symmetry solutions when compared to previously published and current benchmarks. In addition to these new analyses, we provide new large-scale CCSD(T) and explicitly-correlated CCSD(T)-F12 benchmarks for the PDA and PBT systems, which comprise the most complete and accurate calculations of linear polarizabilities and second hyperpolarizabilities on these systems to date. These new CCSD(T) and CCSD(T)-F12 benchmarks confirm our DFT results and emphasize the importance of broken-symmetry effects when calculating polarizabilities and hyperpolarizabilties of π-conjugated chains.

Introduction Since the early 1960s,1-3 the linear and nonlinear optical (NLO) properties of πconjugated polymers have garnered immense interest from both theorists4,

5

and

experimentalists6, 7 for their use as novel optical materials. Specifically, recent developments in these NLO materials have led to a wide variety of technological advancements including optical memory, holography, optical computing, nonlinear microscopy, and electro-optic waveguide devices.8 Predictive computational design, particularly with quantum chemical methods, will play an important role in these advancements by providing a rational and guided path for accurately calculating the NLO properties in these conjugated materials. Among the various quantum chemical techniques currently in use, the most accurate calculations of polarizabilities and second hyperpolarizabilities of polymers have been obtained with wavefunction-based techniques (such as coupled cluster methods); however, these calculations have been limited to short polyenes due to their high computational costs.9 In contrast, the simple Hartree-Fock


Page 3 of 36

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


method is sometimes still employed for the calculation of hyperpolarizabilities of large molecules; unfortunately, the Hartree-Fock formalism by definition does not fully account for dynamic electron correlation, which plays an important role in the estimation of polarizabilities and hyperpolarizabilities (see Refs. 10, 11 and references therein). Kohn-Sham density functional theory (DFT), which includes an approximate treatment of electron correlation, has become an extremely popular and versatile method for its excellent balance between accuracy and computational cost;12 nevertheless, the selection of appropriate exchange-correlation (XC) functionals for calculating polarizabilities and hyperpolarizabilities still remains. Historically, the accurate calculation of polarizabilites in π-conjugated systems has presented an immense challenge for conventional DFT methods. In particular, previous work by others4, 5, 13-15 has demonstrated that DFT calculations with common XC-functionals, including LDA, GGA, and functionals constructed with a fixed percentage of exact exchange, dramatically fail to provide an accurate description of both the linear polarizability (α) and the second hyperpolarizability (γ) of one-dimensional polymers. These difficulties stem from the fact that polarizabilities are second-order electronic properties and, as mentioned previously, will be extremely sensitive to exchange-correlation approximations in Kohn-Sham DFT methods (in particular, for computing γ). In a very recent paper by Champagne and co-workers,16 the authors employed both standard and non-empirically tuned range-separated functionals (specifically CAM-B3LYP and LC-BLYP) to calculate the longitudinal linear polarizability and second hyperpolarizability of polydiacetylene (PDA) and polybutatriene (PBT) oligomers (Figure 1). To test the accuracy of these various range-separated methods, the authors utilized a limited set of coupled-cluster methods with single and double excitations and perturbative triple excitations,



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

CCSD(T), in conjunction with second-order and fourth-order Møller-Plesset perturbation theory (MP2 and MP4) as benchmark standards. Based on these benchmarks, the authors concluded with the following statements: (1) “it is not expected that [adjusting the fraction of short-range HF exchange] will improve the polarizability and the second polarizability,” (2) “for all levels of approximation, the overestimation of the α values [was observed],” and (3) “the bad performance of all levels of approximation to estimate γ of PBT chains [was observed].”

Figure 1: Molecular structures of the PBT and PDA oligomers (N = 1 – 6) studied in this work. The longitudinal linear polarizability and second hyperpolarizability of both structures are computed along the z-axis shown in the figure.

To shed additional light on these previous conclusions, we present a new, detailed investigation using (1) non-empirically tuned range-separated functionals that include a portion of short-range HF exchange, and (2) an extensive analysis of broken-symmetry effects in rangeseparated functionals for calculating polarizabilities and second hyperpolarizabilies in PDA and PBT chains. In regards to the first point, previous work by our group17 and others18-21 has suggested that the inclusion of some short-range HF exchange does improve the accuracy of computed excited-state properties, and we demonstrate that this also enhances the accuracy of


Page 4 of 36

Page 5 of 36

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


computed polarizabilities. Addressing the second point, broken-symmetry effects arise when the restricted (closed-shell) wavefunction becomes unstable towards an unrestricted (open-shell) solution.22,

23

In particular, we find that the lowest-energy electronic states for PBT are not

closed-shell singlets, and enhanced accuracy with range-separated DFT can be obtained by allowing the system to relax to a lower-energy broken-symmetry solution. To both supplement and verify our findings, we also provide new large-scale CCSD(T) and explicitly correlated CCSD(T)-F12 benchmarks for the PDA and PBT systems. It is worth mentioning that in the original work by Champagne and co-workers,16 only a limited set of CCSD(T) benchmark calculations were carried out due to the immense computational cost of these wavefunctionbased methods (i.e. the largest PDA and PBT chains with 5 and 6 oligomers were not computed). For this reason, the authors commented that their range-separated calculations could not be checked to assess if the extrapolated DFT trends would either degrade or improve as a function of oligomer size. In this work, we complete these computationally intense CCSD(T) calculations as well as provide a new set of explicitly-correlated CCSD(T)-F12 benchmarks which comprise the most complete and accurate calculations of polarizabilities and second hyperpolarizabilities on these systems to date. Taken together, our new broken-symmetry range-separated DFT calculations in conjunction with these high-level CCSD(T) and CCSD(T)-F12 benchmarks highlight the importance of broken-symmetry effects when calculating polarizabilities and hyperpolarizabilties of π-conjugated chains. Finally, we give a detailed analysis for all of these effects on various PDA and PBT oligomers and discuss the implications of both short-range exchange and broken-symmetry effects in calculating polarizabilities in these challenging systems.



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 36

Theory and Methodology Since one of the main purposes of this work is to assess the accuracy of various rangeseparated functionals for computing polarizabilities, we briefly review the underlying theory for these methods. In contrast to conventional hybrid functionals, the range-separated formalism24, 25 mixes short range density functional exchange with long range Hartree-Fock exchange by partitioning the electron repulsion operator into short and long range terms (i.e., the mixing parameter is a function of electron coordinates). In the most general form of the range-separated approach, the interelectronic Coulomb operator is given by:26, 27 1 1 − α + β ⋅ erf ( µ ⋅ r12 )  α + β ⋅ erf ( µ ⋅ r12 )  = + . r12 r12 r12

(1)

The erf term denotes the standard error function, r12 is the interelectronic distance between electrons 1 and 2, and µ is the range-separation parameter in units of Bohr-1. The parameters, α and β, satisfy the following relations: 0 ≤ α + β ≤ 1, 0 ≤ α ≤ 1, and 0 ≤ β ≤ 1. The parameter α in the partitioning allows a contribution of HF exchange over the entire range by a factor of α, and the parameter β allows us to incorporate long-range asymptotic HF exchange by a factor of (α + β). For example, the CAM-B3LYP functional of Yanai and co-workers28 uses α = 0.19, α + β = 0.65, and µ = 0.33; however, the CAM-B3LYP functional does not incorporate a “full” rangeseparation as it only has 65% HF exchange at long-range (instead of the correct 100% asymptotic HF exchange). In our previous work on range-separated functionals,29-33 we have used and parameterized “full” range-separation schemes that correspond to setting α = 0.0 and β = 1.0. In particular, we31 and others

34, 35

have previously shown that maintaining a full 100%

contribution of asymptotic HF exchange is essential for accurately describing valence excitations


Page 7 of 36

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


in even relatively simple molecular systems. For the two range-separated LC-BLYP methods used in this work, we fix α + β = 1.0 (with different values of α) in conjunction with selfconsistently tuning the range-separation parameter µ by satisfying DFT-Koopmans’ theorem.36-38 In summary, this theorem states that the energy of the highest occupied molecular orbital (HOMO) equals the negative of the ionization potential (IP), which is defined as the (∆SCF) difference between the ground-state energy of the N electron and the N − 1 electron systems. Within the Kohn-Sham DFT formalism, this condition is fulfilled for the exact XC-functional; therefore, adjusting the range-separation parameter in this self-consistent manner provides a theoretical justification for this procedure. Although several numerical schemes exist, one practical approach for self-consistently tuning the range-separation parameter µ is to numerically minimize the following function: µ µ J 2 ( µ ) = ε HOMO ( N ) + IPµ ( N )  + ε HOMO ( N + 1) + IPµ ( N + 1)  , 2

2

(2)

µ where ε HOMO ( N ) is the HOMO energy of the N-electron system, and IPµ ( N ) is the ground state

energy difference between the N and N – 1 electron systems with the same range-separation parameter. The second term in this equation takes into account the N + 1 system to indirectly tune the LUMO energy of the N electron system. The LUMO energy cannot be directly incorporated in this equation since DFT-Koopmans’ theorem does not explicitly relate the electron affinity (EA) to the negative of the LUMO energy. To obtain the optimal µ values for each oligomer, several single-point energy calculations on fixed geometries (discussed further below) were carried out by varying µ from 0.0 to 0.4 (in increments of 0.02) for each of the N, N + 1, and N – 1 electron states. These calculations were computed using two non-empirically tuned LC-BLYP methods: the first method does not include



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

any short-range exchange (i.e., α = 0.0, β = 1.0), and the second method contains 20% exchange over the entire range (i.e., α = 0.2, β = 0.8). It is important to emphasize that each of the different LC-BLYP parametrizations used in this work still recover the full 100% exchange at asymptotic distance (α + β = 1.0) even though each parametrization has a different exchange contribution at short range. In addition, we also tested the performance of the LC-BLYP (µ = 0.47) method since there has been recent work demonstrating that hyperpolarizabilities of various chromophores are more accurately described without utilizing the non-empirically tuned procedure.39 For all of the oligomers, we carried out a DFT stability analysis at all LC-BLYP levels of theory to converge (if possible) towards a lower-energy, broken-symmetry solution, which allows for an unrestricted spin state as well as a reduction in symmetry of the orbitals. With the broken-symmetry solutions in hand, was computed (Eq. 2) as a function of µ for each polyene and all of the various monomers. Spline interpolation was subsequently used to refine the minimum for each individual system. All DFT calculations were carried out with the Gaussian 09 package

40

using default SCF convergence criteria (density matrix converged to at

least 10-8) and the default DFT integration grid (75 radial and 302 angular quadrature points). In order to maintain a consistent comparison with the previous study of Champagne and co-workers, identical molecular geometries obtained from Ref. 16 were used throughout this work, and the Cartesian coordinates for all the systems are listed in the Supporting Information for completeness. Similarly, we utilized the same 6-31+G(d) basis set to compute the longitudinal static polarizability and second hyperpolarizability of the various PDA and PBT chains, ranging from one to six monomer units. Following the same approach by Champagne et al., all of the polarizabilities obtained with all DFT methods were calculated using analytical


Page 8 of 36

Page 9 of 36

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


derivatives of the energy with respect to field strength within the coupled-perturbed Kohn-Sham (CPKS) method. Second hyperpolarizabilities for all DFT methods were evaluated as secondorder numerical derivatives of the polarizability with respect to the applied external electric field. At the wavefunction-based CCSD(T) and explicitly-correlated CCSD(T)-F12 levels of theory, polarizabilities and hyperpolarizabilities were manually calculated with a custom-developed code as the second- or fourth-order numerical derivatives of the energy with respect to the applied external electric field. In these finite field approaches, the following field amplitudes in atomic units (1 a.u. = 5.142206 10 V/m) were chosen: 0.0, 1 10 , 2

10 , 4 10 , 8 10 , 16 10 , and 32 10 a. u. To maintain a consistent comparison with the various DFT methods, the CCSD(T) calculations were also performed with the same 6-31+G(d) basis set. In addition to the CCSD(T)/6-31+G(d) calculations, new CCSD(T)-F12 calculations were also carried out to verify both the overall trends and quality of the CCSD(T) benchmarks. The CCSD(T)-F12 methods have attracted recent attention for their ability to calculate extremely accurate electronic energies (typically at a higher level of accuracy than conventional CCSD(T) with the same basis), by constructing a wavefunction that depends explicitly on the interelectronic coordinates. As such, the explicitly-correlated CCSD(T)-F12 methods exhibit dramatic improvements in basis set convergence, and results of quintuple-zeta quality that were obtained with smaller triple-zeta basis sets have previously been shown.41 To this end, we carried out our CCSD(T)-F12 calculations with the cc-pVDZ basis, which is the largest basis set available commensurate with both the density-fitting algorithm in Molpro and our computational resources. As a side note, both the CCSD(T) and CCSD(T)-F12 calculations were extremely computationally intensive, especially for the large PDA[5], PDA[6], PBT[5], and



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 36

PBT[6] structures. For example, the largest of these structures, PDA[6], took up to 18 continuous days (for each of the finite field F amplitudes) on 16 × 2.3 GHz AMD Opteron CPUs, and each calculation consumed up to 356 GB of disk space on rapid-access solid state drive storage. Finally, all of the CCSD(T)-F12 polarizabilities and second hyperpolarizabilites reported in this work were obtained from self-consistent CCSD(T)-F12a energies. The CCSD(T)-F12a energies were chosen over the CCSD(T)-F12b results since extensive benchmarks have shown that the CCSD(T)-F12a method gives better results for smaller basis sets (such as the cc-pVDZ basis set used in this work) than CCSD(T)-F12b.42 The F12b variant differs from the F12a method by the inclusion of an additional energy correction which approximately doubles the magnitude of the coupling between the conventional and explicitly correlated pieces of the calculation.42, 43 Taken together, the CCSD(T)-F12 polarizabilities and second hyperpolarizabilities offer a second check on both the DFT and CCSD(T) results as well as provide new, additional high-quality benchmarks for the PDA and PBT systems.

Results and Discussion Figure 2 shows the smooth curves that result from computing J2 as a function of µ for (a) PDA without short-range exchange (LC-BLYPα=0.0,β=1.0), (b) PDA including short-range exchange (LC-BLYPα=0.2,β=0.8), (c) PBT without short-range exchange (LC-BLYPα=0.0,β=1.0), and (d) PBT including short-range exchange (LC-BLYPα=0.2,β=0.8). As discussed in the Theory and Methodology Section, we carried out a full DFT stability analysis at all LC-BLYP levels of theory for both PDA and PBT to allow (if possible) for a lower-energy, broken-symmetry solution. As shown in Figure 2, only closed-shell solutions were obtained for PDA (regardless of

µ value), whereas a broken-symmetry configuration was obtained for large values of µ in PBT


Page 11 of 36

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


for both LC-BLYPα=0.0,β=1.0 and LC-BLYPα=0.2,β=0.8. These results can be rationalized from the chemical structures of these systems since PDA is composed of successive double-single-triplesingle CC bonds (cf. Figure 1) and is more alternant than PBT, which exhibits a single-doubledouble-double bond CC bond pattern (i.e., the π orbitals are much more strongly conjugated along the backbone of PBT compared to PDA). From a more theoretical viewpoint, rangeseparated functionals with higher values of µ inherently contain larger contributions of HF exchange in the XC potential, and it is well-known44-48 that DFT methods containing a large percentage of HF exchange will favor a lower-energy, broken-symmetry configuration.44-48 It is interesting to note that for sufficiently large values of µ, all PBT oligomers (even the smallest PBT[2] structure) will exhibit a broken-symmetry solution where the alpha and beta spin densities alternate through the whole backbone of the molecule, and a long-range ordering of the spin density persists as the length of the oligomer increases. This phenomenon corresponds to electrons localizing in the p orbitals of the carbon atoms in an antiferromagnetic pattern, which can be visualized as the spin density differences in Figure 3. Note that we do not obtain an antiferromagnetic pattern for PDA since the ground states for all of the PDA structures (regardless of µ) are closed-shell singlets, and the alpha and beta spatial distributions are the same.



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 36

Figure 2: Plots of J2 as a function of mu for (a) PDA without short-range exchange (LCBLYPα=0.0,β=1.0), (b) PDA including short-range exchange (LC-BLYPα=0.2,β=0.8), (c) PBT without short-range exchange (LC-BLYPα=0.0,β=1.0), and (d) PBT including short-range exchange (LCBLYPα=0.2,β=0.8). The shaded regions in (c) and (d) denote the values of µ where a symmetrybroken solution is obtained.


Page 13 of 36

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


Figure 3: Spin density difference (blue = positive spin density and red = negative spin density) obtained with LC-BLYP (µ = 0.47) for the various PBT oligomers.

The optimally-tuned µ values for PDA and PBT as obtained by both the LCBLYPα=0.0,β=1.0 and LC-BLYPα=0.2,β=0.8 functionals are summarized in Table 1. It is interesting to note that the optimal µ values for both PDA and PBT are not affected by the broken-symmetry solutions since the shaded regions (where a broken-symmetry solution is obtained) in Figure 2 lie to the right of the minima of all the J2 curves. In other words, the ground-state wavefunction at the optimal µ values for PBT are closed-shell singlets and are not within the shaded regions in



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 36

Figure 2. From these tabulated data entries, we observe an inverse correlation between the optimal µ values and the length of the oligomers, where µ decreases as the number of monomers increases. This general size-dependence in other chemical systems has been previously reported by several groups20,

49-52

and also in the previous study on PDA and PBT in Ref. 16. Our

optimally-tuned µ values are in excellent agreement with the LC-BLYPα=0.0,β=1.0 calculations given in the previous study by Champagne et al., shown in parentheses in Table 1. Specifically, we observe a negligible difference of ~0.005 Bohr-1, and this deviation can be attributed to the difference between both tuning schemes: Ref. 16 directly employs the DFT-Koopmans’ theorem, where µ is tuned such that the energy of the –HOMO is equal to the IP, whereas we obtain µ by adjusting both the energy of the –HOMO to the IP and, in an indirect manner, the LUMO energy via the N + 1 electron system. All subsequent linear polarizability and second hyperpolarizability calculations were carried out using all of the optimal µ values listed in Table 1.

PDA N

PBT

LC-BLYP

LC-BLYP

LC-BLYP

LC-BLYP

(α = 0.0, β = 1.0)

(α = 0.2, β = 0.8)

(α = 0.0, β = 1.0)

(α = 0.2, β = 0.8)

1

0.289 (0.294)

0.240

2

0.241 (0.246)

0.196

0.254 (0.260)

0.205

3

0.212 (0.216)

0.170

0.213 (0.214)

0.168

4

0.194 (0.194)

0.153

0.185 (0.183)

0.142

5

0.181 (0.178)

0.142

0.166 (0.160)

0.125


Page 15 of 36

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


6

0.172 (0.166)

0.134

0.151 (0.145)

0.111

Table 1: Non-empirically-tuned µ values for all oligomers of PDA and PBT at the LC-BLYP/631+G(d) level of theory. The numbers in parentheses are µ values obtained by Champagne et al.16

Tables 2 and 3 summarize the static longitudinal linear polarizability and second hyperpolarizability for the PDA and PBT oligomers, respectively, computed by the CAMB3LYP, LC-BLYPα=0.0,β=1.0, LC-BLYPα=0.2,β=0.8, LC-BLYP(µ=0.47), CCSD(T), and CCSD(T)F12 methods. The broken-symmetry (BS) results for CAM-B3LYP and LC-BLYP(µ=0.47) are also given in the tables for direct comparison to their closed-shell counterparts – recall that the optimally-tuned LC-BLYP methods for both PDA and PBT were not affected by the brokensymmetry solutions. In addition, we also tabulated the values for both CAM-B3LYP (BS) and LC-BLYP (BS) (µ=0.47) (cf. Table SI-12 in the Supporting Information) and found that while increases with size for PBT[2] – PBT[6], these systems are more accurately characterized as diradicals, without higher-lying spin states contributing to the trends in polarizability (cf. Figure SI-1, Tables SI-13, and SI-14 for LC- BLYP (BS) (µ=0.47) calculations of the triplet state in the Supporting Information). All mean absolute errors (MAE) were computed with respect to the CCSD(T)/6-31+G(d) benchmarks to allow for a consistent comparison with the DFT calculations that were computed with the same 6-31+G(d) basis set. It is worth mentioning that the discrepancy between our results and previous CCSD(T) calculations of the hyperpolarizability16 arise from the different numerical methodology and computational hardware used in these previous studies. Specifically, Champagne and co-workers utilized a Romberg differentiation procedure53 for their calculations, whereas we directly fitted the CCSD(T) energies (as a function of the applied external field) to second- and fourth-order



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 36

polynomials to obtain the linear polarizability and second hyperpolarizability, respectively. We found that the Romberg procedure was extremely sensitive to very small energy differences (even as small as 10-8 Hartrees), which can easily arise from different versions of software and machine architectures used in these previous studies. As such, all of the CCSD(T) results reported in this paper utilize the polynomial fitting procedure, which we found to be more numerically stable, for obtaining the linear polarizability and second hyperpolarizability. Furthermore, to check for possible non-dynamical correlation effects in our CCSD(T) calculations, we also computed the T1 diagnostic for all of the oligomers and obtained T1 values ranging from 0.0125 – 0.0131 and 0.0156 – 0.0164 for PDA[1] – PDA[6] and PBT[2] – PBT[6], respectively (T1 values greater than 0.02 indicate that a multireference electron correlation method is necessary54). In addition, we also carried out broken-symmetry unrestricted CCSD(T) (UCCSD(T)) calculations for all of the PBT oligomers (cf. Fig. SI-2 in the Supporting Information) and found that all of the UCCSD(T) energies using the broken-symmetry HF reference determinant were larger than their restricted CCSD(T) counterparts, further verifying (in conjunction with the T1 diagnostic discussed previously) that correlation effects are properly handled at the single-reference restricted CCSD(T) level of theory. Finally, we also did additional CCSD(T) calculations using the non-diffuse triple-zeta 6-311G(d,p) basis (cf. Table SI-15 in the Supporting Information) to assess the convergence of the 6-31+G(d,p) results used as benchmarks in both our study and the previous study by Champagne. We find that the presence of diffuse functions in the 6-31+G(d,p) basis plays a larger role than the additional valence basis functions in the 6-311G(d,p) basis, which is consistent with previous studies by Champagne and co-workers.5,

14

Furthermore, the explicitly-correlated CCSD(T)-F12


Page 17 of 36

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


calculations closely mirror the overall trends in the CCSD(T)/6-31+G(d) methods, giving additional indication of the basis-set convergence of our calculations.



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

N

CAM-B3LYP

CAM-B3LYP (BS)

Page 18 of 36

LC-BLYP

LC-BLYP

(α = 0.0, β = 1.0)

(α = 0.2, β = 0.8)

LC-BLYP

LC-BLYP (BS)

µ = adj.

µ = adj.

µ = 0.47

µ = 0.47

CCSD(T)

CCSD(T)-F12

α (a.u.)

1

140.14

140.14

139.14

139.50

131.73

131.73

123.33

123.92

(123.33) 2

323.26

323.26

327.87

327.03

290.62

290.62

267.20

272.71

(267.03) 3

558.04

558.04

584.22

577.84

483.61

483.61

439.70

454.19

(439.50) 4

825.12

825.12

893.36

875.60

695.05

695.05

627.67

654.38

(627.39) 5

1111.36

1111.36

1241.75

1206.20

916.30

916.30

824.27

865.75

6

1408.90

1408.90

1618.19

1558.40

1142.92

1142.92

1024.23

1081.85

MAE

176.64

176.64

249.69

229.70

58.97

58.97

116

116

115

112

93

93

110

104

γ ( × 103 a.u.)

1

(105) 2

964

964

982

951

666

666

522

624

(693) 3

3,785

3785

4131

3952

2228

2228

2657

2302

(2344) 4

9,670

9670

11592

10894

4949

4949

5366

5231

(5146) 5

18,913

18913

25091

23088

8658

8658

8468

7909

6

31,133

31133

45508

40903

13086

13086

13866

11220

5599

5599

9405

8152

330

330

MAE

Table 2: Longitudinal linear polarizability and second hyperpolarizability for increasingly large PDA oligomers at various levels of theory. All DFT and CCSD(T) calculations utilized the 631+G(d) basis with the CCSD(T)-F12 calculations using the cc-pVDZ basis and density-fitting approach described in the main text. The abbreviations “BS” and “µ = adj.” indicate a brokensymmetry calculation and a non-empirically tuned value of µ, respectively. Values in parentheses denote the CCSD(T) values obtained by Champagne et al.16


Page 19 of 36


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


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

N

CAM-B3LYP

LC-BLYP

LC-BLYP

CAM-B3LYP

(α = 0.0, β = 1.0)

(α = 0.2, β = 0.8)

LC-BLYP

LC-BLYP (BS)

(BS)

µ = adj.

µ = adj.

µ = 0.47

µ = 0.47

Page 20 of 36

CCSD(T)

CCSD(T)-F12

α (a.u.)

2

322.62

322.62

320.13

323.19

309.09

306.00

271.53 (271.51)

271.39

3

688.28

684.84

693.53

698.06

636.98

596.55

542.98 (542.83)

550.57

4

1196.72

1143.86

1235.13

1236.96

1069.64

955.85

894.94 (894.28)

918.93

5

1836.27

1693.92

1951.31

1942.37

1587.36

1363.60

1311.41

1364.21

6

2589.62

2315.23

2842.01

2811.11

2170.35

1804.02

1777.08

1870.96

MAE

367.11

272.51

448.83

442.75

195.10

45.62

γ ( × 103 a.u.)

2

363

363

403

374

341

552

543 (542)

303

3

2652

4387

2902

2698

2385

3469

3272 (3181)

1909

4

11535

17960

12623

11793

9361

11804

10874 (11042)

9608

5

35815

50818

40298

37788

25697

28113

30172

25531

6

87968

112605

104413

98032

55568

53443

68510

59783

MAE

5312

14624

9658

7760

4004

3652

Table 3: Longitudinal linear polarizability and second hyperpolarizability for increasingly large PBT oligomers at various levels of theory. All DFT and CCSD(T) calculations utilized the 631+G(d) basis with the CCSD(T)-F12 calculations using the cc-pVDZ basis and density-fitting approach described in the main text. The abbreviations “BS” and “µ = adj.” indicate a brokensymmetry calculation and a non-empirically tuned value of µ, respectively. Values in parentheses denote the CCSD(T) values obtained by Champagne et al.16

Figure 4 presents a graphical summary of Tables 2 and 3 by plotting the percent relative error for each of the various DFT methods using the CCSD(T)/6-31+G(d) calculations as benchmarks. Our calculations of the PDA linear polarizabilities (α) in Figure 4(a) using CAMB3LYP, LC-BLYP(µ=0.47), and CCSD(T) are in accordance with the values computed in Ref. 16. As mentioned in the Introduction, the previous work by Champagne and co-workers did not compute the CCSD(T) polarizabilities for the largest PDA and PBT chains with 5 and 6


Page 21 of 36

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


oligomers and, therefore, the overall trends in their range-separated calculations could not be checked to see if the extrapolated DFT trends would either degrade or improve as a function of oligomer size. By completing these computationally expensive benchmarks, we can now state that the same dramatic overestimation of the longitudinal linear polarizability with the nonempirically tuned LC-BLYPα=0.0,β=1.0 functional persists for large oligomers. However, in contrast to previous findings,16 we find that the accuracy of the linear polarizability does improve when a small amount of HF exchange is included at short-range (LC-BLYPα=0.2,β=0.8) compared to the “base” LC-BLYPα=0.0,β=1.0 approach (a reduction in the MAE from 249.69 to 229.70 a.u. is observed). Moreover, the accuracy of LC-BLYPα=0.2,β=0.8 appears to further improve as a function of the number monomer units, N, particularly for N ≥ 3. The CAM-B3LYP functional gives more accurate predictions of the linear polarizability than either of the non-empirically tuned LCBLYP methods, but the LC-BLYP(µ=0.47) functional gives the best agreement (with the lowest MAE values of 58.97 a.u.) compared to CCSD(T) benchmarks, which is consistent with Ref. 16 and previous work on polarizabilities by other groups.39 Turning to the second hyperpolarizabilities (γ) of PDA, Figure 4(b) shows a similar trend in accuracy compared to Figure 4(a) for the linear polarizabilities. As before, incorporating some portion of short-range HF exchange does improve the overall accuracy as a function of oligomer size; however, the most accurate second hyperpolarizabilities are still obtained with the LC-BLYP(µ=0.47) functional (MAE = 330 a.u.). Moreover, as shown in Figures 4(a) and 4(b), there is no change in either α or γ when the stability of the wavefunction is taken into consideration since the ground states of all the PDA oligomers have a closed-shell solution regardless of µ value.



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 36

Figure 4: Percent relative error (compared to CCSD(T)) in α and γ as a function of number of monomer units in PDA (upper panel) and PBT (lower panel) for different levels of theory.

Turning to the PBT oligomers, we find that the overestimation of the static linear polarizability (α) is more severe than for PDA. Similar to Figure 4(a), we also find that the polarizabilities are improved by including some short-range HF exchange, although the difference between these methods is smaller in the PBT system (MAEs for LC-BLYPα=0.0,β=1.0 vs. LC-BLYPα=0.0,β=1.0 in PBT are 448.83 and 442.75 a.u., respectively). However, in contrast to the PDA chains, the PBT oligomers are much more strongly conjugated along their backbone and


Page 23 of 36

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


can converge towards a lower-energy broken-symmetry solution. As mentioned previously, DFT methods that contain a large percentage of HF exchange will favor an unrestricted open-shell configuration, and we find that a lower-energy broken-symmetry solution is actually preferred in CAM-B3LYP and LC-BLYP(µ=0.47) for PBT. Interestingly, as the amount of HF exchange is increased in the XC-functional, the slope in the error of the static linear polarizability (as a function of oligomer size) decreases, as shown in Figure 4(c). However, it is worth noting that increasing the amount of HF exchange to 100% (i.e., pure Hartree-Fock) will lead to severe overestimations for both α and γ as a function of size (cf. Fig. SI-2 and Tables SI-13 and SI-14 in the Supporting Information). In the case of PBT, the characteristic plateau in the asymptotic limit, which still corresponds to a fairly large ~ 20% percent error, is reached for LCBLYP(µ=0.47). However, for the PBT linear polarizabilities, the broken-symmetry solutions give a more accurate result compared to their restricted, closed-shell counterparts in all cases. Specifically, allowing the system to relax to a lower-energy broken-symmetry solution with CAM-B3LYP leads to a constant value for the linear polarizability as N increases, where the associated error is less than 30%. This is in stark contrast to the growing error of ~45% in the restricted, closed-shell CAM-B3LYP method. A remarkably different behavior is obtained with the broken-symmetry LC-BLYP(µ=0.47) method where a negative slope as a function of oligomer size is obtained, resulting in a relative error of less than 5% when N > 4. As mentioned previously, we also carried out LC- BLYP (BS) (µ=0.47) calculations of the triplet state (cf. Fig. SI-1 in the Supporting Information) and found that the errors in the linear polarizability were larger than their LC-BLYP (BS) singlet state counterparts, further confirming that these systems are more accurately characterized as diradicals, without higher-lying spin states contributing to



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 24 of 36

the trends in polarizability. In summary, the MAEs for the linear polarizabilities in PBT can be summarized as follows: α[LC-BLYP(µ=0.47) (BS)]

Polarizabilities of Ï-Conjugated Chains Revisited - ACS Publications

Recommend Documents

Polarizabilities of Ï-Conjugated Chains Revisited - ACS Publications