Analytical Gradients for Density Functional Calculations with

Article pubs.acs.org/JPCA

Analytical Gradients for Density Functional Calculations with Approximate Spin Projection Toru Saito and Walter Thiel* Max-Planck-Institut für Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Mülheim an der Ruhr, Germany S Supporting Information *

ABSTRACT: We have derived and implemented analytical gradients for broken-symmetry unrestricted density functional calculations (BS-UDFT) with removal of spin contamination by Yamaguchi’s approximate spin projection method. Geometry optimizations with these analytical gradients (AGAP-opt) yield results consistent with those obtained with the previously available numerical gradients (NAP-opt). The AGAP-opt approach is found to be more precise, efficient, and robust than NAP-opt. It allows full geometry optimizations for large openshell systems. We report results for three types of organic diradicals and for a binuclear vanadium(II) complex to demonstrate the merits of removing the spin contamination effects during geometry optimization (AGAP-opt vs BS-UDFT) and to illustrate the superior performance of the analytical gradients (AGAP-opt vs NAP-opt). The results for the vanadium(II) complex indicate that the AGAP-opt method is capable of handling pronounced spin contamination effects in large binuclear transition metal complexes with two magnetic centers.

1. INTRODUCTION Spin-unrestricted density functional theory (UDFT) has been widely applied to open-shell systems because of its simplicity.1−3 The broken-symmetry (BS) spin-polarized solution, which represents a mixture of pure low-spin and high-spin (HS) states, accounts for a certain amount of nondynamical correlation in the framework of a single Slater determinant. However, the BS solution is not an eigenfunction of the S2 operator because of spin contamination from the HS states.4 The expectation value ⟨S2̂ ⟩BS can be used to measure the degree of spin contamination. When it is large, UDFT calculations for the low-spin state can give spurious results for energies and geometries. It is well-known that the inherent shortcomings of any chosen exchange-correlation (XC) functional may cause errors in the calculated relative energies and barrier heights.5,6 To evaluate the accuracy of an XC functional, it is desirable to separate the errors caused by spin contamination from those of the XC functional itself. Several projection schemes for the removal of spin contamination from exchange coupling constants7−9 and the potential energy surface (PES) of a BS solution10−13 have been proposed. Yamaguchi and co-workers devised an approximate spin projection (AP) method7,11,13 and a geometry optimization technique based on the AP energy.14,15 They use a numerical finite difference procedure to evaluate the first derivatives of the ⟨Ŝ2⟩BS value (see below), so hereafter we call this approach NAP-opt. The available NAP-opt results indicate that the spin contamination effect should be eliminated from the energy and the energy derivatives when investigating © 2012 American Chemical Society

chemical reactions that involve open-shell systems. The NAPopt method becomes impractical for large systems, since 6N (N: number of atoms) energy evaluations are required for calculating the gradient in each cycle. Also, geometry optimizations with numerical gradients often encounter convergence problems for polyatomic molecules because of numerical inaccuracies. Therefore, analytical gradients are essential for the efficient investigation of large open-shell systems (such as transition metal complexes and active sites of metalloproteins). In this article, we present a geometry optimization technique based on the analytical gradients for the AP method (called AGAP-opt). We have implemented the AGAP-opt method in a development version of the Gaussian16 interface of the ChemShell program package.17,18 To illustrate the merits of the AGAP-opt method, we report four applications (see Figure 1), involving m-xylylene (1), the diradical intermediate in the [2 + 2] cycloaddition of singlet oxygen with a long-chain alkene (2), the 2,6-isomer of didehydropyridine (2.6-pyridyne) (3), and a binuclear vanadium(II) complex (4). Previous studies on these systems have shown that the geometries obtained from BS-UDFT calculations differ significantly from the experimental data or multireference (MR) ab initio results.13,19−23 We compare the AGAP-opt geometries with results from experiment and previous computations. Received: September 7, 2012 Revised: October 5, 2012 Published: October 22, 2012 10864

dx.doi.org/10.1021/jp308916s | J. Phys. Chem. A 2012, 116, 10864−10869

The Journal of Physical Chemistry A

Article

Figure 1. Calculated molecules: m-xylylene (1), 1,4-diradical intermediate (2), 2,6-isomer of didehydropyridine (2,6-pyridyne) (3), and [V(C5H5)]2Pn complex (4).

integration over the spin density.25,26 Thus, we rely on the conventional formula, in which the ⟨Ŝ2⟩BS value is defined as

2. METHOD In the framework of the Heisenberg Hamiltonian, the AP energy (EAP)7,11 is given by E AP = αEBS − (α − 1)EHS

(1)

2

⟨Ŝ ⟩HS − ⟨Ŝ ⟩BS

2 occ occ ∂Sij ∂⟨Ŝ ⟩BS = −2 ∑ ∑ Sij ∂x ∂x i∈α j∈β

(2)

where EX and ⟨Ŝ2⟩X denote the total energy and the expectation value of total spin angular momentum for the spin state X (X = BS,HS), respectively, and SBS z is the z-component of total spin for the BS state defined as

∂Sij

N α − Nβ = (3) 2 α β where N and N represent the number of α and β electrons, respectively. In the AP method, the ⟨S2̂ ⟩HS value is regarded as constant (geometry-independent) because the reference HS state is assumed to be free from spin contamination. The SBS z value is constant since Nα and Nβ are fixed. Thus, the first derivative of EAP with respect to a geometrical variable can be written as BS

∂x

∂cμq ∂x

(

)

⎛ ∂c*

∑ ⎜⎜ μν

μi

⎝ ∂x

Sμνcνj + cμ*i

∂Sμν ∂x

cνj + cμ*iSμν

∂cνj ⎞ ⎟ ∂x ⎟⎠

(8)

=

∑ cμpU pqx p

(9)

where the indices p and q refer to orbitals which are either occupied or virtual. For both the α and β orbitals, the occupiedvirtual blocks of Uxpq are obtained by solving the coupledperturbed Hartree−Fock (CPHF) equations,27 while the occupied-occupied and virtual−virtual blocks are obtained from the orbital orthogonality condition.28 The energy and gradient for the reference HS state (EHS and ∂EHS/∂x), the uncorrected energy and gradient for the target BS state (EBS and ∂EBS/∂x), and the CPHF solutions are computed using Gaussian 09.16 The ∂EAP/∂x term is then evaluated from eqs 7−9. The AP energy and its analytical gradient are used to perform geometry optimizations with the DL-FIND optimizer,29 which is part of the ChemShell package.17,18

(4)

with 2 2 ⟨Ŝ ⟩HS − SzBS(SzBS + 1) ∂⟨Ŝ ⟩BS ∂α = 2 2 2 ∂x ∂x ⟨Ŝ ⟩HS − ⟨Ŝ ⟩BS

AO

=

where Sμν represents an atomic orbital (AO) overlap matrix element and ∂Sμν/∂x its derivative; cμq denotes an MO coefficient. The derivative ∂cμq/∂x term can be expressed as

HS

∂E ∂E ∂E ∂α BS =α − (α − 1) + (E − EHS) ∂x ∂x ∂x ∂x

(7)

where Sij denotes the overlap between α and β molecular orbitals (MOs). The derivative of Sij is written as

SzBS

AP

(6)

in analogy with wave function theory.4 Its derivative is given by

2 ⟨Ŝ ⟩HS − SzBS(SzBS + 1) 2

∑ ∑ |Sij|2 i∈α j∈β

with α=

occ occ

2

⟨Ŝ ⟩BS = SzBS(SzBS + 1) + N β −

(5)

The evaluation of the ⟨Ŝ2⟩ value of UDFT solutions is complicated because of the two-electron nature of the operator. Recently, Nishihara et al.24 provided justification for the conventional scheme based on the single Kohn−Sham determinant, through comparisons with approaches that involve 10865



Article

3. RESULTS AND DISCUSSION 3. 1. m-Xylylene (1). First, we examined m-xylylene (1), which has a triplet ground state. To confirm the validity of the AGAP-opt method, we applied it to 1 at the B3LYP/6-311G** level of theory30 and compared the results with those obtained from the numerical NAP-opt procedure. Table 1 shows the

C2−C1−O1 angle, and the dihedral angle C2−C1−O1−O2 by 0.021 Å, 4.8°, and 4.7°, respectively. Furthermore, we studied the [2 + 2] cycloaddition of singlet oxygen to 1-undecene (C11H22) as an example of long-chain alkenes. The geometry of the corresponding 1,4-diradical intermediate (2) was optimized at the B3LYP/6-31G* level. The numerical NAP-opt approach failed to locate a local minimum for 2, while the analytical AGAP-opt calculation converged smoothly. Such superior performance of AGAP-opt vs NAP-opt is often found in large systems. Compared with the uncorrected BS-UB3LYP/6-31G* results, removal of spin contamination via AGAP-opt leads to changes in the C1−O1 bond length, the C2−C1−O1 angle, and the dihedral angle C2−C1−O1−O2 by 0.025 Å, 5.4°, and 5.5°, respectively. As expected, these changes are similar to those found in the ethylene case (see above). They seem too large to be ignored. The other geometrical parameters are less affected (bond lengths by up to 0.01 Å and bond angles by up to 0.5°, see the optimized structures in Supporting Information for details). 3. 3. 2,6-Pyridyne (3). The equilibrium geometry of 2,6pyridyne (3) may be monocyclic (diradicaloid) or bicyclic (closed-shell).23,31−33 Spin-restricted DFT calculations give the bicyclic form with a short C2−C6 distance (ca. 1.6 Å) and a small C2−N−C6 angle (ca. 70°), while ab initio MR methods yield both bicyclic and monocyclic structures.23,31 Geometry optimizations performed with a state-specific multireference coupled cluster singles and doubles approach (Mk-MRCCSD) and analytic gradients predicted the monocyclic structure to be more stable.31 Here, we assess the performance of UDFT calculations with regard to the equilibrium structure of 3, and investigate the influence of spin contamination on the geometric parameters through comparison with ab initio MR results. We examined three hybrid XC functionals: B3LYP (20% HF exchange), BHandHLYP (50% HF exchange), and LC-ωPBE (100% HF exchange for long-range interactions),30,34 in combination with the cc-pVTZ basis set.35 Calculations for an initial structure close to the bicyclic form always gave the closed-shell solution with ⟨Ŝ2⟩BS = 0 regardless of the chosen XC functional, and subsequent geometry optimizations yielded the bicyclic forms as in the case of previous spin-restricted DFT calculations.31,32 On the other hand, whenever using an initial geometry close to the monocyclic form, the spin-polarized BS solutions were obtained, which led to the monocyclic structures upon geometry optimizations (BS-opt) as seen from the C2− C6 bond length and C2−N−C6 angle (Table 2).

Table 1. Optimized C−C Bond Lengths (Å) in m-Xylylene (1) parametersa

a

geometry

a

b

c

d

BS-opt NAP-opt AGAP-opt extrapolated valuesb

1.389 1.388 1.389 1.389

1.432 1.408 1.409 1.406

1.418 1.406 1.404 1.403

1.400 1.437 1.438 1.440

See Figure 1. bFrom ref 21.

structural parameters with (NAP-opt, AGAP-opt) and without (BS-opt) corrections; the optimized Cartesian coordinates are presented in the Supporting Information. It is obvious that the NAP-opt and AGAP-opt results for the bond lengths are in good agreement, while being quite different from the uncorrected BS-opt values in the case of b-d. Very recently, Malrieu et al.21 proposed a simple method to predict an optimized geometry based on the AP energy, without explicitly evaluating AP gradients. The extrapolated parameters obtained with their procedure at the same level of theory are also listed in Table 1. Their results agree well with the present NAP-opt and AGAP-opt results. 3. 2. 1,4-Diradical Intermediate (2). For the [2 + 2] cycloaddition of singlet oxygen (1Δg O2) to olefins, previous UDFT calculations suggested a stepwise mechanism involving a diradical intermediate.13,19,20 One earlier study applied the NAP-opt method to the 1,4-diradical intermediate found in the reaction of ethylene (C2H4) with 1Δg O2 at the B3LYP/6311+G** level of theory.30 We reoptimized the diradical intermediate by means of the AGAP-opt method at the same level of theory. The differences in structural parameters between the NAP-opt and AGAP-opt results are at most 0.004 Å for bond lengths and 1.0° for angles, respectively. The deviation in the total energy is negligibly small (less than 0.01 kcal mol−1) confirming that the numerical approach is useful and practical for small molecules. Compared with the uncorrected BS-UB3LYP/6-311+G* results, the application of the AGAP-opt method changes the C1−O1 bond length, the Table 2. Optimized Structural Parameters in 2,6-Pyridyne (3)a B3LYPc parameters

b

r(N−C2) r(C2−C3) r(C3−C4) r(C2−C6) ∠C2−N−C6 ∠N−C2−C3 ∠C2−C3−C4 ∠C3−C4−C5

BHandHLYP

LC-ωPBE

BS-opt

BS-opt

AGAP-opt

BS-opt

AGAP-opt

Mk-MRCCSDd

1.302 1.381 1.400 2.217 116.8 125.7 115.9 119.9

1.296 1.370 1.386 2.183 114.8 127.0 115.8 119.5

1.325 1.366 1.383 2.059 102.0 134.8 116.6 115.2

1.296 1.372 1.387 2.189 115.3 126.8 115.7 119.8

1.327 1.368 1.384 2.025 99.5 136.5 116.4 114.5

1.336 1.373 1.391 2.017 98.0 137.5 116.8 113.5

The cc-pVTZ basis set was used. bIn Å for bond lengths and degree for angles. cThe AGAP-opt method did not find a local minimum corresponding to the monocyclic form. dFrom ref 31. a

10866



Article

N−C6 angles of less than 100°, the EAP value gradually increased again. This behavior occurred even when starting from the BS-opt minimum structure obtained by the relaxed PES scan. Hence the fully relaxed spin-contamination-free B3LYP PES has only one minimum (the bicyclic structure). Analogous results were found when using pure DFT methods and other hybrid XC functionals with a small amount of HF exchange. Apparently, inclusion of a large amount of HF exchange term is required to stabilize the monocyclic form of 3. Such hybrid functionals correctly reproduce the equilibrium geometry of 3 obtained from state-of-the-art Mk-MRCCSD calculations, provided that spin contamination is removed during the optimization by the AGAP-opt approach. 3. 4. [V(C5H5)]2Pn Complex (4). Finally, we applied the AGAP-opt method to the binuclear vanadium(II) complex [V(C5H5)]2Pn (Pn2‑ = dianionic pentalene ligand) (4). Numerical NAP-opt optimizations are prohibitively expensive for such a large system. The localized S = 3/2 spins at each vanadium(II) ion can be coupled to singlet, triplet, quintet, and septet states. According to experimental36 and theoretical22 studies, the ground state is a singlet and has a formal triple bond. However, the strength of the V−V interaction depends on the chosen XC functional. Gorelsky22 concluded that hybrid XC functionals are less suitable for this multinuclear 3d transition metal complex. In fact, the V−V distance obtained with the widely used BS-UB3LYP approach was 0.45 Å larger than the experimental value, while pure DFT calculations gave a bond length in accordance with the X-ray data (2.54 Å). This type of phenomenon is found in many calculations on transition metal complexes.37−39 Recent studies suggest that inclusion of an empirical dispersion correction (DFT-D) can improve the geometries in such complexes,40,41 but spin contamination may also be responsible for poor performance. Here, we determine how much the structure for the BS singlet (S = 0) state is affected by spin contamination from the septet (S = 3) state, which is the major cause for the overestimation of the V−V distance. Strictly speaking, spin contamination from the intermediate spin (S = 1, S = 2) states should also be considered, but the AP method is incapable of removing these effects. The BP86, B3LYP, and M06 functionals30,42,43 were used together with the SVP basis set.44 Figure 3 shows the natural orbitals corresponding to the bonding σ, π, and δ orbitals. The occupation numbers obtained at each optimized geometry are also shown in Figure 3. The ⟨S2̂ ⟩BS values obtained at the UBP86, UB3LYP, and UM06 level are 1.531, 2.820, and 2.797, respectively. In the case of UBP86, the σ and π orbitals remain intact during the geometry optimization, and the optimized V−V distance reproduces the X-ray data well. On the other hand, the BS solutions for UB3LYP and UM06 both show strong spin

However, the BS-opt values for the C2−C6 bond length and C2−N−C6 angle differ significantly from the reference MkMRCCSD results, being too large by more than 0.15 Å and 15°, respectively. This discrepancy may be due to the insufficient accuracy of the XC functionals and/or spin contamination from the triplet state, which reduces the interaction between two local spins in 3. We note that all of the present XC functionals predict the triplet state to be the ground state at the BS-opt structures. To clarify the origin of these discrepancies, we performed NAP-opt and AGAP-opt geometry optimizations as well as relaxed PES scans along the C2−N−C6 angle. These latter scans were done with the Gaussian 09 program package,16 optimizing all other geometric parameters with the BS-opt method. Single-point HS and AP energies were then calculated at the BS-opt structures. Figure 2 shows the resulting energy

Figure 2. Potential energy profiles along the C2−N−C6 angle; all other geometric parameters were optimized for the BS (S = 0) state at the BHandHLYP/cc-pVTZ level. The HS (S = 1) and AP (S = 0) energies were calculated at the same level and the same BS-opt structures.

profiles computed at the BHandHLYP/cc-pVTZ level. Similar plots are obtained for B3LYP and LC-ωPBE (see Figures S1 and S2 in the Supporting Information). The usual BS-opt procedure converges to a spurious local minimum with a large C2−N−C6 angle because of the excessive stabilization caused by the triplet ground state. The NAP-opt method failed to locate a local minimum in each case, indicating that the ∂⟨Ŝ2⟩BS/∂x term should be properly evaluated even for such small molecules. The AGAP-opt method converged to local minima with the BHandHLYP and LC-ωPBE functional: the C2−C6 bond length (the C2−N−C6 angle) was dramatically reduced by 0.124 Å (12.8°) for BHandHLYP and 0.164 Å (15.8°) for LC-ωPBE, respectively (see Table 2). The AGAP-opt structures are quite similar to the optimized Mk-MRCCSD geometry. Application of the AGAPopt method with BHandHLYP and LC-ωPBE thus corrected the potential energy surface by removing the spin contamination effect, providing local minima that are in good agreement with the Mk-MRCCSD results. However, in the case of B3LYP, the AGAP-opt method did not find a stationary point corresponding to the monocyclic form: both EAP and ∂EAP/∂x were lowered at first, but for C2−

Figure 3. Natural orbitals for σ, π, and δ bonding orbitals and their occupation numbers obtained by UBP86, UB3LYP, and UM06 calculations at each optimized structure. 10867



Article

Tests on four open-shell singlet molecules demonstrate the superior performance of the AGAP-opt method. The NAP-opt and AGAP-opt results are consistent with each other for mxylylene (1) and the investigated 1,4-diradical (2). The spin contamination effect on the geometry is found to be significant in the case of 2,6-pyridyne (3) and the metal−metal bond in the binuclear vanadium(II) complex (4). The magnitude of this effect depends on the chosen XC functional, which may be relevant for the validation of XC functionals (especially in transition-metal chemistry). The AGAP-opt method has been implemented in the ChemShell package and can thus be applied in combined quantum mechanics/molecular mechanics (QM/MM) studies,45 for example on metalloproteins. This work is in progress.

polarization in these orbitals, indicating significant spin contamination caused by the septet state. Thus, we carried out AGAP-opt geometry optimizations using the septet state as the HS state. Dispersion-corrected40 BP86-D and B3LYP-D calculations were performed for comparison. Table 3 Table 3. Optimized V−V Bond Length (Å) in the [V(C5H5)]2Pn Complex functional

geometry

r(V−V)

BP86

BS-opt AGAP-opt BS-opt AGAP-opt BS-opt AGAP-opt BS-opt AGAP-opt BS-opt AGAP-opt

2.55 2.50 2.48 2.44 2.96 2.76 2.84 2.61 2.87 2.72 2.54

BP86-D B3LYP B3LYP-D M06 Expt.a a

■

ASSOCIATED CONTENT

S Supporting Information *

Optimized Cartesian coordinates for 1−4 and energy profiles from relaxed PES scans along the C2−N−C6 angle for B3LYP and LC-ωPBE (Figures S1 and S2). This material is available free of charge via the Internet at http://pubs.acs.org.

From ref 36.

■

summarizes the optimized V−V distances with (AGAP-opt) and without (BS-opt) corrections. The spin contamination effects on the V−V distance are 0.05, 0.20, and 0.15 Å for the BP86, B3LYP, and M06 functionals, respectively. The results correlate with the magnitude of spin polarization of the σ and π orbitals. Removing the spin contamination effect of the septet state (AGAP-opt vs BS-opt) shortens the V−V bond length by 0.04−0.05 Å for the pure functionals (BP86, BP86-D) and by 0.15−0.23 Å for the hybrid XC functionals (B3LYP, B3LYP-D, M06). Inclusion of dispersion corrections further shortens the V−V distance by 0.06−0.07 Å (BP86 vs BP86-D) and 0.12− 0.15 Å (B3LYP vs B3LYP-D). Comparing with the experimental value (2.54 Å), an uncorrected BS-UBP86 optimization yields almost perfect agreement (2.55 Å), which deteriorates after correcting for spin contamination and dispersion (2.44 Å). On the other hand, the strongly overestimated BS-UB3LYP value (2.96 Å) moves much closer to the experimental value after including both corrections (2.61 Å). The BS-UM06 results (2.87 Å) also benefits from removing spin contamination, but to a lesser extent (2.72 Å). In a pragmatic sense, one may argue that uncorrected BS-UDFT optimizations in such systems are appropriate for pure DFT functionals such as BP86, for which spin contamination effects are small. The significant changes found for the hybrid XC functionals (B3LYP and M06) strongly indicate that the spin contamination effect should be removed in these cases. It is reassuring that B3LYP gives a reasonable V−V distance (2.61 Å), only slightly larger than the experimental value (2.54 Å), after accounting for the large effects of spin contamination and dispersion.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS T.S. is grateful to the Japan Society for the Promotion of Science (JSPS) for a Postdoctoral Fellowship for Research Abroad.

■

REFERENCES

(1) Neese, F. Coord. Chem. Rev. 2009, 253, 526. (2) Solomon, E. I.; Scott, R. A.; King, R. B., Eds.; Computational Inorganic and Bioinorganic Chemistry; Wiley: New York, 2009. (3) Siegbahn, P. E. M. J. Biol. Inorg. Chem. 2006, 11, 695. (4) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry; Dover Publications, Inc.: New York, 1996. (5) Zhao, Y.; González-Garcia, N.; Truhlar, D. G. J. Phys. Chem. A 2005, 109, 2012. (6) Zheng, J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2007, 3, 569. (7) Yamaguchi, K.; Takahara, Y.; Fueno, T. In Applied Quantum Chemistry; Smith, V. H., Schaefer, H. F., III, Morokuma, K., Eds.; D. Reidel: Boston, MA, 1986; p 15. (8) Noodleman, L.; Davidson, E. R. Chem. Phys. 1986, 109, 131. (9) Ruiz, E.; Cano, J.; Alvalez, S.; Alemany, P. J. Comput. Chem. 1999, 20, 1391. (10) Schlegel, H. B. J. Chem. Phys. 1986, 84, 4530. (11) Yamaguchi, K.; Jensen, F.; Dorigo, A.; Houk, K. N. Chem. Phys. Lett. 1988, 149, 537. (12) Li, J.; Noodleman, L. In Spectroscopic Methods in Bioinorganic Chemistry, ACS Symposium Series; Solomon, E. I., Hodgson, K. O., Eds.; American Chemical Society: Washington, DC, 1998; Vol. 692, p 179. (13) Saito, T.; Nishihara, S.; Kataoka, Y.; Nakanishi, Y.; Kitagawa, Y.; Kawakami, T.; Yamanaka, S.; Okumura, M.; Yamaguchi, K. J. Phys. Chem. A 2010, 114, 7967. (14) Kitagawa, Y.; Saito, T.; Ito, M.; Nakanishi, Y.; Shoji, M.; Koizumi, K.; Yamanaka, S.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Chem. Phys. Lett. 2007, 442, 445.

4. CONCLUSIONS In this work, we have derived and implemented analytical gradients for BS-UDFT calculations with approximate spin projection. The analytical gradient evaluation includes a computationally demanding CPHF step, but it is more precise and less costly than the numerical finite-difference evaluation, especially for large systems (with 6N energy calculations). Geometry optimizations have been carried out both with analytical (AGAP-opt) and numerical (NAP-opt) gradients. 10868



Article

(15) Saito, T.; Nishihara, S.; Kataoka, Y.; Nakanishi, Y.; Matsui, T.; Kitagawa, Y.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Chem. Phys. Lett. 2009, 483, 168. (16) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision C.01; Gaussian, Inc.: Wallingford, CT, 2009. (17) Sherwood, P.; de Vries, A. H.; Guest, M. F.; Schreckenbach, G.; Catlow, C. R. C.; French, S. A.; Sokol, A. A.; Bromley, S. T.; Thiel, W.; Turner, A. J.; Billeter, S.; Terstegen, F.; Thiel, S.; Kendrck, J.; Rogers, S. C.; Casci, J.; Watson, M.; King, F.; Karlsen, E.; Sjøvoll, M.; Fahmi, A.; Schäfer, A.; Lennartz, C. J. Mol. Struct. (THEOCHEM) 2003, 632, 1. (18) Chemshell, a Computational Chemistry Shell; Science &Technology Facilities Council: Swindon, U.K.; www.chemshell.org. (19) Yoshioka, Y.; Tsunesada, T.; Yamaguchi, K.; Saito, I. Int. J. Quantum Chem. 1997, 65, 787. (20) Maranzana, A.; Ghigo, G.; Tonaccini, G. J. Am. Chem. Soc. 2000, 122, 1414. (21) Malrieu, J. -P.; Trinquier, G. J. Phys. Chem. A 2012, 116, 8226. (22) Gorelsky, S. I. J. Chem. Theory Comput. 2012, 8, 908. (23) Debbert, S. L.; Cramer, C. J. Int. J. Mass Spectrom. 2000, 201, 1. (24) Nishihara, S.; Saito, T.; Yamanaka, S.; Kitagawa, Y.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Mol. Phys. 2010, 108, 2559. (25) Wang, J.; Becke, A. D.; Smith, V. H. J. Chem. Phys. 1995, 102, 3477. (26) Cohen, A. J.; Tozer, D. J.; Handy, N. C. J. Chem. Phys. 2007, 126, 214104. (27) Gerratt, J.; Mills, J. M. J. Chem. Phys. 1968, 49, 1719. (28) Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Int. J. Quantum Chem. 1979, 13, 225. (29) Kästner, J.; Carr, J. M.; Keal, T. W.; Thiel, W.; Wander, A.; Sherwood, P. J. Phys. Chem. A 2009, 113, 11856. (30) (a) Becke, A. Phys. Rev. A 1988, 38, 3098. (b) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785. (c) Becke, A. J. Chem. Phys. 1993, 93, 5648. (31) Prochnow, E.; Evangelista, F. A.; Schaefer, H. F., III; Allen, W. D.; Gauss, J. J. Chem. Phys. 2009, 131, 064109. (32) Chattopadhyay, S.; Chaudhuri, R. K.; Mahapatra, U. S. Chem. Phys. Lett. 2010, 491, 102. (33) Jagau, T. -C.; Gauss, J. Chem. Phys. 2012, 401, 73. (34) Vydrov, O. A.; Scuseria, G. E. J. Chem. Phys. 2006, 125, 234109. (35) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (36) Jones, S. C.; O’Hare, D. Chem. Commun. 2003, 2208. (37) Bühl, M.; Reimann, C.; Pantazis, D. A.; Bredow, T.; Neese, F. J. Chem. Theory Comput. 2008, 4, 1449. (38) Sandala, G. M.; Hopmann, K. H.; Ghosh, A.; Noodleman, L. J. Chem. Theory Comput. 2011, 7, 3232. (39) Minenkov, Y.; Singstad, Å.; Occhipinti, G.; Jensen, V. R. Dalton Trans. 2012, 41, 5526. (40) Grimme, S. J. Comput. Chem. 2006, 27, 1787. (41) Lonsdale, R.; Harvey, J. N.; Mulholland, A. J. J. Phys. Chem. Lett. 2010, 1, 3232. (42) Perdew, J. P. Phys. Rev. B 1986, 33, 8822. (43) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215. (44) Schäfer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571. (45) Senn, H. M.; Thiel, W. Angew. Chem., Int. Ed. 2009, 48, 1198. 10869


Analytical Gradients for Density Functional Calculations with

Recommend Documents