Article pubs.acs.org/JCTC
Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Reduced Basis Set Dependence in Anharmonic Frequency Calculations Involving Localized Coordinates Magnus W. D. Hanson-Heine* School of Chemistry, University of Nottingham, University Park, Nottingham NG7 2RD, U.K. S Supporting Information *
ABSTRACT: Localized normal coordinates are known to be effective in speeding up anharmonic frequency calculations by reducing the complexity of the nuclear Hamiltonian and wave function. Displacing atoms in localized coordinates can also cause relatively small changes in the electronic structure, which can be exploited for further computational efficiency improvements during ab initio electronic structure calculations of the potential energy surface by reducing the electronic basis set dependence. Three different schemes for reducing the basis set dependence have been investigated in this work. These include combining localized coordinate schemes with general mixed basis sets, distance based force-field reductions, and using coordinate specific basis sets. The importance of accurately describing electronic interactions is found to diminish both for multicoordinate terms involving the displacement of remote atoms and when describing the interactions between more remote atoms within specific coordinates.
■
INTRODUCTION Three main computational bottlenecks exist when solving the nuclear Schrö dinger equation and calculating vibrational frequencies. The dimensions of the potential energy surface (PES) and the nuclear wave function scale exponentially with the number of nuclear degrees of freedom, while the cost of the underlying electronic calculations used to describe the PES scales similarly with the number of electrons present. While it is the aim of this work to address this third bottleneck, concerning the electronic structure calculations that underlie the PES, the methodology employed here is significantly informed by the advancements that have made in solving the first two. The most common approach for reducing the impact of these first two problems is to use the harmonic approximation.1 In this case the nuclear potential energy surface is assumed to be an exact quadratic function around an optimized geometry, reducing the necessary PES evaluation to second derivatives of the energy with respect to nuclear displacements. The nuclear vibrational Schrödinger equation then factorizes into 3N-6 (or 3N-5 for a linear molecule) noninteracting one-dimensional quantum harmonic oscillators for which exact wave function solutions are known, each contained along a different orthogonal normal coordinate. Vibrational frequencies, wave functions, and normal coordinates can then be computed routinely for large molecules from the eigenvalues and eigenvectors of the mass-weighted Hessian matrix, (H), corresponding to the nuclear Cartesian displacements for each atom Hij =
1 ∂ 2E(x) mimj ∂xi ∂xj
where xi represents a Cartesian displacement coordinate of an atom with mass mi. When higher accuracy anharmonic corrections are needed, the PES is often approximated as a truncated Taylor series expansion in normal coordinates taken around an optimized molecular geometry where derivatives of the energy with respect to nuclear displacement are evaluated up to the fourth order, referred to as the quartic force field (QFF).2 As a further approximation, energy derivatives involving more than a specific number of normal modes can also be excluded, which is referred to as an n-mode representation, or nMR.2−4 Similarly, additional approximations are made to the nuclear wave function, including taking harmonic solutions in normal coordinates and adding anharmonic effects either as perturbations4−10 or by using configuration interaction.11,12 Vibrational self-consistent field theory (VSCF)13−18 can also be used to create one-mode wave functions that see an averaged potential from the other modes, and these “modals” can then be correlated using perturbation theory (VMP2),3,19−21 configuration interaction (VCI),22−26 or coupled cluster (VCC)27,28 theories. Recently, unitarily transformed normal coordinates optimized for maximum spatial localization29−31 have been employed during VSCF calculations (L-VSCF) and found to reduce the number of terms needed to compute the anharmonic potential energy surface for a given level accuracy by increasing the spatial decay of terms that involve coupling between coordinates.32 Similarly, the anharmonic coupling between normal coordinates can be partially represented using Received: October 27, 2017 Published: January 31, 2018
(1) © XXXX American Chemical Society
A
DOI: 10.1021/acs.jctc.7b01075 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
methods were applied to the pyrene, C20H22, C10H12, and octatetraene (C8H10) molecules. The VSCF and L-VSCF methods were then applied using a nuclear basis set consisting of the 20 lowest harmonic oscillator wave functions along each of the modes, with a convergence threshold set to 1 × 10−6 cm −1 . Transition frequencies were obtained from the eigenvalues of the ground state modals. Anharmonic thirdand fourth-order derivatives of the energy with respect to nuclear displacement were used within a Taylor series expansion of the QFF PES, calculated using numerical differentiation of the Hessian matrices with a step size of 0.5291 Å, and the anharmonic frequencies were calculated using a 2MR representation of the QFF unless otherwise stated. The 6-31G, 3-21G, and STO-3G basis sets were also used when calculating anharmonic terms (vide inf ra), representing reductions in basis set quality. These compact Pople basis sets were chosen due to their computational efficiency compared with other larger basis sets available, as well as the increased scope for comparison with previous work in the field.32,47 In the case of C10H12, the reduced system size and mode-coupling allowed for comparison using the cc-pVTZ and cc-pVDZ basis sets, and these basis sets were used in place of 631G(d,p) and STO-3G, respectively, during that part of the analysis. Where different basis sets are used for harmonic and anharmonic terms, the equilibrium geometry will be different in the different basis sets. However, this has been observed to be insignificant for semirigid systems, with mixed method work having previously made use of numerical wavenumber anharmonicities calculated at one level of theory being applied to harmonic frequencies calculated at a different level,53 hybrid schemes combining more accurate methods for the harmonic force field with lower level anharmonic third- and fourthderivatives without compensation for this displacement,53,54 and basis set extrapolation schemes for anharmonic frequencies using different local mode formulations.55 Reduced basis set anharmonic terms have therefore been computed at the reference geometry of the full basis set. The localized modes used throughout this work, Q̃ , have been calculated using a unitary transformation, U, of normal modes, Q, which leads to the maximization of an objective function (ξ)29
only one-mode localized terms when including bilinear harmonic coupling terms.33 Localized modes have also been used as the basis for calculating two-dimensional infrared spectra (2D-IR) in the framework of exciton theory.34 Alternately, optimized coordinates that minimize the ground state or state averaged energy during a VSCF calculation have also been used in combination with VSCF (oc-VSCF),35−37 including variational optimization in internal coordinates.38 Correlations between the optimized or localized VSCF modals used to describe the nuclear wave function have also been accounted for using a number of schemes, including VCI expressed in localized coordinates (L-VCI),39 quasi-degenerate perturbation theory (oc-VQDPT),40 and coupled cluster theory (oc-VCC).41,42 Anharmonic mode−mode correlations have been shown in these cases to be smaller in localized or VSCF energy optimized coordinates, with faster convergence seen with respect to expansion of the excitation space describing the wave function.39,41,42 Methods have also been developed which alter the extent to which high order excitation space expansions of the wave function are needed in cases where harmonic coupling terms dominate, either by accounting for differences between pairs of related vibrational states in normal and optimized coordinates, as with the L-VSCF(HC) method,43 or by applying a variety of constrained coordinate optimizations including intermediate-constrained localized coordinates (ILVSCF),44−46 windowed localized coordinates (WL-VSCF),47 or hybrid optimized-localized coordinates (HOLC-VSCF).48 These techniques have focused on reducing the first two bottlenecks in anharmonic frequency calculations. However, localized vibrational modes and their associated coordinates often involve significantly displacing relatively few atoms compared to their normal coordinate counterparts. Cheng and Steele noted that the motion of a small number of atoms within a sufficiently large system may be expected to cause a relatively insignificant perturbation to the overall electronic structure.32 These small changes should also make computational efficiency improvements possible in the electronic structure calculations that define the nuclear PES. One way in which this is possible is by reducing the computational load during ab initio electronic structure calculations by reducing the basis set dependence of anharmonic calculations through lowering the number of basis functions needed to describe the electronic structure for regions or atoms that are not expected to significantly change the overall energy. In this work three different schemes for reducing the basis set dependence during anharmonic calculations through the use of localized vibrational coordinates have been investigated within the LVSCF framework, ranging from combining localized coordinate schemes with general mixed basis sets, using coordinate specific basis set combinations tailored to individual coordinates, and analyzing distance based force-field reductions in the basis set used to describe coupling between spatially remote modes. This final scheme shares some conceptual features with work being done on electronic many-body systems where it is known that the n-body expansion of molecular clusters converges more rapidly if point charges from atoms considered outside the region of interest are used in the energy calculations.49,50
Q̃ = QU
(2)
Two possible objective functions were originally proposed.29 The first function follows from the molecular orbital localization scheme of Pipek and Mezey by maximizing the sum of the squares of the atomic contributions to the localized modes31 ⎞2 ⎛ 2⎟ ⎜ ̃ ξat = ∑ ∑ ⎜ ∑ (Q iα , μ) ⎟ ⎠ μ i ⎝ α=x ,y,z N
(3)
where μ and i are localized modes and nuclei, respectively, and N is the total number of nuclei. The second function follows the orbital localization scheme of Boys by maximizing the distance between the centers of the modes30
■
COMPUTATIONAL DETAILS Equilibrium structures, normal coordinates, and harmonic frequencies were calculated using the 6-31G(d,p) basis set the B97-1 exchange-correlation functional,51 with a developmental version of the Q-Chem software package.52 These
ξdist
⎞2 ⎛N 2 = ∑ ⎜⎜∑ ∑ (Q̃ iα , μ) R i⎟⎟ ⎠ μ ⎝ i α=x ,y,z
(4)
where Ri is the position vector of each nucleus with respect to the molecular origin. B
DOI: 10.1021/acs.jctc.7b01075 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation The resulting localized modes can then be used to transform the mass-weighted Hessian matrix from atomic Cartesian displacement coordinates into localized mode coordinates H̃ = U†Q†HQU
(5)
The unitary transformations used in forming the localized coordinates were determined by repeated Jacobi sweeps until iteration no longer increased the objective function by more than 1 × 10−6, and the frequencies associated with the diagonal elements of the Hessian matrix expressed in localized mode coordinates are referred to here as pseudoharmonic frequencies. Convergence problems can arise during the VSCF procedure in cases where the coordinate system being used is insufficient to properly describe a particular type of motion, such as torsional or bending modes in the case of rectilinear coordinates,56 or when the PES expansion order is insufficiently high, either of which can happen when using localized coordinates. In order to minimize this risk, convergence was aided by slowing down the change in the effective potential during the VSCF procedure by computing the effective potential at each new cycle from the weighted sum of the unaltered new potential and the effective potential obtained during the previous iteration56
Figure 1. Pyrene with atoms using reduced basis functions marked in red/orange. Shown for the full basis set and three mixed basis set combinations denoted mix1, mix2, and mix3 in order of reducing size.
respectively, and have been used when calculating the thirdand fourth-order terms during the following anharmonic calculations. The results of the mixed basis set approaches are shown in terms of mean absolute deviation (MAD) from the full basis set results in Table 1 and Table 2, with signed error distribution Table 1. Mean Absolute Differences between the Mixed Basis Set and Full Basis Set VSCF Frequencies of Pyrene (in cm−1) basis STO-3G mix1 mix2 mix3 3-21G mix1 mix2 mix3 6-31G mix1 mix2 mix3
Um,s(Q m) = b[Um,s(Q m)]old + (1 − b)[Um,s(Q m)]new (6)
where Um,s(Qm) is the effective mean-field potential of mode m in state s, the subscripts new and old denote the current and previous VSCF cycles, respectively, and b is a real number, set to 0.8, that defines the extent of the mixing between the two potentials.
■
RESULTS AND DISCUSSION A. Mixed Basis Sets. The normal mode vibrations in the hydrogen stretching region of polycyclic aromatic hydrocarbons (PAHs) are delocalized across the C−H bound hydrogen atoms, yet they remain localized with respect to the remaining carbon atoms. This separation between the atom types leads to a chemically intuitive way of separating the basis set. The C−H stretching vibrational modes of PAHs have recently become of increased astrophysical interest since their neutral and cationic forms were proposed as candidates for the previously unidentified infrared bands observed in the interstellar medium, and these modes are therefore a potential target for mixed basis set and partial Hessian approaches.57−60 Previous work localizing C−H modes of small molecules has also indicated that localized C−H stretching modes tend to have less mode coupling than the normal modes and perform particularly well when treated within the L-VSCF and post L-VSCF framework. Pyrene has therefore been chosen as a test case to see if coordinate localization methods can be used to improve the accuracy of mixed basis set calculations. In addition to the complete basis set calculations, nine mixed basis set combinations have been used to test the sensitivity of pyrene anharmonic frequencies in the different coordinates. These mixed calculations employed a lower quality basis set of either 6-31G, 3-21G, or STO-3G quality to describe (i) the carbon atoms at the center of the molecule, (ii) all carbon atoms not directly bound to hydrogen atoms, and (iii) all of the carbon atoms in the system. These mixed basis sets, shown in Figure 1, have been designated “mix1″, “mix2″, and “mix3″,
all modes
>3000 cm−1
1k−3k cm−1
3000 cm−1
1k−3k cm−1
