Article pubs.acs.org/JPCA
Approximate Force Constants from Uncoupled Self-Consistent Field Perturbation Theory Using Nonhybrid Density Functional Theory Published as part of The Journal of Physical Chemistry virtual special issue “Mark S. Gordon Festschrift”. Zhigang Ni,†,§ Krzysztof Wolinski,‡ and Peter Pulay*,† †
Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701, United States Department of Theoretical Chemistry, Faculty of Chemistry, Maria Curie-Sklodowska University, 20-031 Lublin, Poland
‡
S Supporting Information *
ABSTRACT: Second derivatives of the molecular energy with respect to the nuclear coordinates (the nuclear Hessian or force constant matrix) are important for predicting infrared and Raman spectra, for calculating thermodynamic properties, for characterizing stationary states, and for guiding geometry optimization. However, their calculation for larger systems scales with molecular size one power higher than the calculation of the energy and the forces. The step responsible for the steep scaling of the nuclear Hessian is the coupled-perturbed self-consistent field (CP-SCF) iteration. This is omitted in the uncoupled SCF (UC-SCF) approximation. We have found that, though UC-SCF performs rather poorly at the Hartree−Fock and hybrid DFT levels, its performance for “pure” (non-hybrid) DFT is remarkably good. This is valid also for imaginary frequencies that characterize transition states. UC-SCF vibrational frequencies and normal modes are compared with coupled calculations for various exchange−correlation functionals including Hartree−Fock, and with basis sets ranging from simple to large for a variety of organic and some organometallic molecules. Their unexpectedly good performance makes them good candidates for calculating thermodynamic properties and for guiding difficult geometry optimizations, including the determination of transition states.
1. INTRODUCTION The second derivatives of the molecular energy with respect to external perturbations are important in computational quantum chemistry.1 A number of molecular properties, among others polarizabilities, NMR chemical shifts, and spin−spin coupling constants, are second-order quantities. Most important are the second derivatives with respect to nuclear coordinates, the (quadratic) force constants, also called collectively the nuclear Hessian. They determine the harmonic vibrational frequencies and associated thermodynamic quantities of molecules. They are the key quantities for characterizing stationary points on molecular potential energy surfaces as minima, transition states, or higher-order saddle points. One of their most important roles is to guide geometry optimization. For this purpose, and also for most thermodynamic quantities, approximate values suffice. For the equilibrium structures of standard organic molecules, Hessians from molecular mechanics work generally well, and the gain in the optimization step seldom justifies the cost of the quantum mechanical evaluation of the Hessian.1 However, for more challenging optimization problems, in particular for the location of transition states, and for the characterization of stationary points on molecular potential surfaces, accurate Hessians are necessary. Unfortunately, they are also the computationally most expensive second-order © XXXX American Chemical Society
properties, particularly for large systems. In principle, they can be calculated by numerically differentiating the molecular energy twice. However, this is such an expensive and numerically sensitive procedure that the calculation of force constants for polyatomic molecules became practicable only after the introduction of gradient techniques, as the numerical first derivatives of the forces on the nuclei.2,3 Fully analytical second derivatives with respect to the nuclear coordinates offer further advantages in computational efficiency and numerical precision. However, the computational scaling with molecular size remains the same as for the numerical differentiation of the gradients, a fact not widely acknowledged. In this paper we will concentrate on SCF wave functions, as these are the techniques applicable for large systems. Their theory was developed fully by Bratož4 in 1958 but this paper was so much ahead of the computational realities of that time that it remained virtually unknown. This was also the fate of the first implementation by Swanstrøm et al.5 Gerratt and Mills6 advocated using the analytical derivatives of the Hellmann− Feynman7,8 force. Although this is not a good approximation Received: October 31, 2016 Revised: November 23, 2016 Published: December 2, 2016 A
DOI: 10.1021/acs.jpca.6b10959 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
various common SCF methods (Hartree−Fock, hybrid density functional theory (DFT), and pure DFT that does not contain Hartree−Fock exchange) by comparing them to coupledperturbed SCF (CP-SCF) results. The CP-SCF results yield the true second derivatives of the molecular energy. We will also investigate whether the error of the UC-SCF scheme depends significantly on the quality of the one-electron basis set. Although the UC-SCF results are quite accurate for pure density functional methods, they are not exact and cannot be applied for all purposes. Section 4 discusses the potential applications of inexpensive UC-SCF Hessians, and methods to eliminate their main deficiencies.
for force constants, their formulation of the key computational step, the coupled-perturbed SCF for nuclear perturbations, became the standard in the field. The first practical analytical second derivative program for geometry perturbations was developed by Pople and co-workers.9 This paper also introduced a practical method for solving the coupledperturbed Hartree−Fock equation by a combination of iterative and finite methods. This technique is applicable also for density functional theory (DFT); we will call it the coupled-perturbed self-consistent field (CP-SCF) method. For large molecules, the computational cost of SCF analytical second derivatives is determined by the CP-SCF step, which has to be carried out for each nuclear coordinate perturbation. The formal computational cost of CP-SCF in the usual atomic basis set SCF method scales as O(N4) for a single perturbation, where N is the number of basis functions, and is proportional to the molecular size at constant basis set quality. For the full Hessian, the cost scales as O(MN4), where M is the number of nuclei. This is one power higher than the energy or gradient. In practice, natural sparsity limits the scaling to ∼O(N3) or asymptotically O(N2 ) in the SCF step but the ratio T(Hessian)/T(SCF) is still proportional to the number of nuclei for larger molecules. For smaller molecules, the O(N4) overhead masks the formal fifth-order scaling of the CP-SCF procedure. In addition, the CP-SCF procedure for large molecules requires a large amount of fast memory. If this is not available, the calculation must be carried out in batches with the corresponding reduction of computational efficiency. The efficiency of Hessian calculations can be much increased, particularly for large molecules, if the coupled-perturbed SCF step, i.e., the back-coupling of the first-order change in the CPSCF equations, is omitted. This procedure is also called the sum-over-states (SOS) method, as it can be formulated by summing over all single excitations in the SCF procedure (conversely, the CP-SCF method is known in the physics literature as the random-phase approximation). The terms coupled and uncoupled Hartree−Fock were introduced by Dalgarno10 for perturbation independent basis sets. This method included additional approximations that were revised and analyzed by Langhoff et al.11 However, omitting the coupling between orbitals introduces an error. This error is large for polarizabilities (it may be as large as factor of 2),12 and this discouraged the use of the uncoupled SCF (UC-SCF) procedure. The error is also significant for nuclear Hessians at the Hartree−Fock and hybrid DFT levels, although smaller. We have found, to our surprise, that uncoupled Hessians are quite accurate for pure (nonhybrid) density functionals. Though we have no clear theoretical explanation, this is not totally unexpected. For magnetic (i.e., imaginary) perturbations, the coupling between orbitals is caused by the Hartree−Fock exchange term, and it vanishes for pure DFT. In this case, SOS is equivalent to the full CP-SCF procedure. Pure (nonhybrid) density functional theory offers another big potential advantage. The most expensive part of the calculations, the evaluation of the Coulomb energy and its derivatives, can be speeded up very much (up to 2 orders of magnitude) by expanding the electron density in an intermediate plane wave basis.13−15 To our knowledge, no analytical Hessian program exists that makes use of a plane wave expansion, although several gradient programs have been implemented.16,17 Section 2 summarizes the theory, and in section 3 we investigate the performance of the UC-SCF procedure for
2. THEORY In this section, we briefly summarize our formalism for analytical second derivatives of the SCF energy with perturbation dependent basis sets. We prefer the density matrix formalism18 we used for, e.g., NMR chemical shifts,19 as it has no problems with degeneracy as some orbital-based formulations have. The second derivative of the closed-shell SCF energy with respect to two perturbations, denoted by the superscripts a and b is given20 by the equation below. Note that a term is missing in ref 20. 1 G(D,g ab)D 2 1 + Tr (ha + G(D, g a ))Db − Sa (DbFD + DFbD + DFDb) 2
ab Eab = V nuc + Tr habD − 2SabW +
Here Tr⟨A⟩ denotes a matrix trace, V is the nucleus−nucleus repulsion, h is the one-electron part of the Hamiltonian, and G(D,g) is the two-electron part that depends both on the unperturbed first-order reduced density matrix D = 2CCT and on the basis functions {g}. W = 2CεCT is the weighted density matrix, C is the occupied part of the orbital coefficient matrix, F = h + G is the Fock matrix, and the superscript T denotes transposition. The derivative Fb = hb + G(Db,g) + G(D,gb) includes the dependence of the Fock matrix on the nuclear coordinates through both the basis functions and the density matrix. For DFT, the exchange part of G is replaced by the exchange−correlation matrix. There are equivalent formulas that are symmetrical in the perturbation.20 The first trace in the above equation contains only the zeroth-order density together with the derivatives of the Hamiltonian and the basis functions. The second trace contains derivatives of the density matrix. They have to be calculated for each of the 3N nuclear perturbations and this makes the evaluation of the Hessian expensive for large systems. The density matrix derivatives, e.g., Db, are determined by solving the CP-SCF equations: 1 Db = − DSbD + 2 (CiCTv
+
∑ CTi (Fb − εiSb)Cv (εi − εv)−1 i,v
Cv CTi )
This part requires the introduction of the occupied orbital vectors Ci and the virtual vectors Cv, which are columns of SCF orbital coefficient matrix. Fb depends on Db, and therefore the equation has to be solved iteratively. In the UC-SCF procedure, we neglect the dependence of Fb on Db, and determine the density matrix derivative as B
DOI: 10.1021/acs.jpca.6b10959 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
Table 1. Comparison of Coupled (CP-SCF) SCF and Uncoupled SCF (UC-SCF) Harmonic Vibrational Frequencies for transButadiene, gauche-Dimethyl Disulfide, and Naphthalene Using Different Density Functionals and Polarized Triple-ζ Basis Setsa fingerprint
X−H method
RMS
MSD
MAD
RMS
MSD
low MAD
RMS
MSD
MAD
94 37 13 12 13
238 100 13 12 9
227 91 11 9 7
227 91 11 9 7
66 25 10 11 12
244 92 10 11 13
232 88 −4 −9 −11
232 88 9 8 11
102 41 18 18 18
201 87 9 9 7
201 87 9 9 7
201 87 9 9 7
C4H6 H−F B3LYP BLYP BP86 PBE
142 91 83 88 87
142 90 80 85 84
142 90 80 85 84
H−F B3LYP BLYP BP86 PBE
147 92 82 87 83
147 90 78 83 79
147 90 78 83 79
H−F B3LYP BLYP BP86 PBE
147 98 88 94 92
147 98 88 94 91
147 98 88 94 91
120 94 47 37 22 13 21 12 22 12 CH3SSCH3 75 66 30 25 18 8 19 8 19 7 C10H8 108 102 45 41 26 17 24 16 25 17
a The butadiene and naphthalene calculations used the cc-pVTZ basis; the dimethyl disulfide, the def2-tzvp basis. Root-mean-square (RMS), mean signed (MSD), and mean absolute (MAD) deviations (UC minus CP) are given in cm−1 separately for the C−H stretching vibrations, for the fingerprint region (350−2000 cm−1), and for low vibrations (
