We investigate the use of Hartree-Fock and density functional perturbative corrections for estimating the counterpoise correction (CPC) for interaction energies at the self-consistent field level. We test our approach using several popular basis sets on the S22 set of weakly bound systems, which can exhibit large basis set superposition errors. Our results show that the perturbative approaches typically recover over 95% of the CPC and can be up to twelve times faster to compute than the conventional methods and therefore provide an attractive alternative to calculating CPCs in the conventional way.P.M.W.G. thanks the Australian Research Council for
funding (Grant Nos. DP0984806 and DP1094170) and APAC
for a generous allocation of supercomputer resources. J.D.
thanks the ANU/RSC for a PhD scholarship

Deng, Jia

Gilbert, Andrew T. B.

Gill, Peter M. W.

The Journal of Chemical Physics

The Australian National University

THE JOURNAL OF CHEMICAL PHYSICS 135, 081105 (2011)Communication: Efficient counterpoise correctionsby a perturbative approachJia Deng, Andrew T. B. Gilbert, and Peter M. W. Gilla)Research School of Chemistry, The Australian National University, ACT 0200, Australia(Received 9 August 2011; accepted 11 August 2011; published online 26 August 2011)We investigate the use of Hartree-Fock and density functional perturbative corrections for estimatingthe counterpoise correction (CPC) for interaction energies at the self-consistent field level. We testour approach using several popular basis sets on the S22 set of weakly bound systems, which can ex-hibit large basis set superposition errors. Our results show that the perturbative approaches typicallyrecover over 95% of the CPC and can be up to twelve times faster to compute than the conventionalmethods and therefore provide an attractive alternative to calculating CPCs in the conventional way.© 2011 American Institute of Physics. [doi:10.1063/1.3632054]I. INTRODUCTIONThe interaction energy of a weakly bonded system is of-ten found by the so-called supermolecular approach, whichuses the energy differenceEint = EAB(A ∪ B) − EA(A) − EB(B), (1)between the entire system and its subsystems, where sub-scripts denote the (sub)system, parentheses denote the basisset, and the fragment structures are taken from the complex.1However, this straightforward approach is complicated by thebasis set imbalance between the system and its subsystems;subsystem A is modelled by the basis functions A ∪ B in thesupermolecule, but only those in A in the isolated subsystem.This inconsistency systematically overestimates the interac-tion energy and gives rise to the so-called “basis set super-position error” (BSSE).2, 3 The BSSE is largest when smallbasis sets are used, and can become large enough to prevent aqualitative characterization of the interaction energy.Because BSSE is a consequence of basis set incomplete-ness, it vanishes in the complete basis set limit. Therefore, astraightforward solution to reducing this error is to performthe calculations in sufficiently large bases. However, such abrute force approach is impractical for large systems and al-ternative methods for removing the BSSE must be used. Themost important of these are the a priori Chemical Hamil-tonian Approach (CHA),4 Symmetry-Adapted PerturbationTheory (SAPT),5 and the a posteriori counterpoise correc-tion (CPC) of Boys and Bernardi.6 Because of its conceptualsimplicity, and ease of implementation, the CPC is the mostwidely used method and is also the focus of this present work.For a system, AB, consisting of two subsystems A and B,the CPC is defined asECPC = EA(A) − EA(A ∪ B) + EB(B) − EB(A ∪ B),(2)where, for example, EA(A ∪ B) is the energy of A calculatedusing the basis functions for both A and B. Because the basisa)Electronic mail: peter.gill@anu.edu.au.functions associated with B have no nuclei in that calculation,they are referred to as ghost functions.The cost of calculating the CPC is typically greater thanthat of calculating the original interaction energy. Given thatthe CPC provides only an estimate of the BSSE, it seemswise to consider approximate methods which mitigate thisadditional cost. Recently, we developed perturbative meth-ods that “bootstrap” a small-basis calculation toward a cor-responding large-basis calculation7–10 and these approaches(and the related dual-basis methods11–16) have proven usefulin approximating large-basis results at greatly reduced cost.Errors due to this formalism have been shown to be small,while reducing the computational cost by factors between3 and 10. Originally proposed for Hartree-Fock theory, thetechniques have been extended to density functional theory,9second-order Møller-Plesset perturbation theory (MP2)(Refs. 10 and 14) and, recently, to complete active space self-consistent field (CAS-SCF) theory.17Although the ghost functions are a large fraction of thebasis in a CPC calculation, their energetic contribution issmall and such calculations are therefore ideal candidatesfor our perturbative methods or, indeed, for other dual-basisschemes.11–16 In this Communication, we use our Hartree-Fock perturbative correction (HFPC) (Refs. 7 and 8) anddensity functional perturbative correction (DFPC) (Ref. 9)schemes to evaluate the effect of the ghost functions in CPCcalculations.Because of their inability to capture long-range disper-sion effects, many SCF methods are poorly suited to describ-ing weakly bound systems. However, this deficiency is ame-liorated in the empirical DFT-based methods such as DFT-D(Refs. 18–20) and, because of their cost effectiveness, suchmethods are likely to continue to be popular for the studyof systems with several dozen atoms. Accordingly, we focuson HF and DFT levels of theory and apply these to the S22dataset21 of weakly bounded complexes.II. THEORYThe theory and performance of HFPC and DFPC calcula-tions are described in detail in Refs. 7–9. Here we outline the0021-9606/2011/135(8)/081105/4/$30.00 © 2011 American Institute of Physics135, 081105-1Downloaded 26 Sep 2011 to 150.203.35.38. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions081105-2 Deng, Gilbert, and Gill J. Chem. Phys. 135, 081105 (2011)HFPC scheme and refer the interested reader to our earlierpapers for details.An HFPC calculation begins with a converged SCF cal-culation in a primary basis of n functions, yielding molecularorbitals and an associated density matrix P. P is then used tobuild the new Fock matrixF[1]ab = hab +n∑λσPλσ[(ab|λσ ) − 12(aλ|bσ )], (3)where λ, σ denote primary basis functions and a, b denotefunctions in a larger (N > n) secondary basis. Diagonaliza-tion of F[1] yields improved molecular orbitals and an im-proved density matrix P[1]. The HFPC energyEHFPC =N∑abP[1]ab hab +12N∑abcdP[1]ab P[1]cd [(ab|cd) −12(ac|bd)],(4)is more accurate7, 8 than that from other dual-basis methodsbecause of its use of P[1]P[1], rather than PP[1].Although HFPC (and DFPC) are able to jump betweenany two basis sets, the selection of an optimal primary basisfor a given secondary basis and accuracy threshold normallyrequires careful calibration studies. However, in the case ofthe CPC calculation of EA(A ∪ B), it is it is natural to adoptA as the primary basis and A ∪ B as the secondary basis, be-cause the ghost functions, B, perturb the energy only slightly.Moreover, such a basis pairing leads to a large N/n ratiowhich yields significant savings in computational time.III. RESULTSA. Computational detailsWe considered three standard basis sets: 6-31G(d), 6-31+G(d), and 6-31+G(2df, p). For each of these, we cal-culated the exact CPC using fully converged SCF calcula-tions and used these results as our reference values. The CPCswere also estimated using our perturbative schemes, which re-quired evaluating each of the energies in Eq. (2) using eitherthe HFPC or DFPC method. The accuracy was measured bythe percentage of the CPC that was recovered by the pertur-bative method.For the DFT calculations we used the B3LYP hybridfunctional and the SG-1 quadrature grid.22 In addition to theBSSE, Kohn-Sham CPC calculations are also affected by ananalogous grid superposition error (GSE). Both the GSEs forthe conventional and perturbative calculations were correctedin the same way using a counterpoise correction.CPC values were calculated for the weakly bound sys-tems in the S22 dataset,21 which features a diverse selectionof non-bonded interactions and a range of system sizes. Allcalculations were performed using a locally modified copy ofthe Q-Chem program23 and the raw energies of all the systemsare provided in the supplementary material.24B. AccuracyColumns 2–7 of Table I list the CPCs and the percent-age recovered by the HFPC method for all the systems inthe S22 set and for the three basis sets considered. BSSE re-flects basis set incompleteness and our results confirm thatCPCs typically decrease from 6-31G(d) to 6-31+G(d) to 6-31+G(2df, p). This effect led Antony and Grimme to arguethat polarized triple-ζ basis sets are sufficiently large that onecan dispense with the CPC altogether25 but, unfortunately,such basis sets are impractical for large systems.Increasing the basis set size also systematically improvesthe relative performance of the HFPC method. Larger basissets reduce the usefulness of the ghost functions and therebyallow HFPC to capture greater proportions of the BSSE. Thisimprovement in relative performance is also a feature of theoriginal HFPC method.7 We see that, even for the small 6-31G(d) basis, where the CPCs are the largest, the HFPCmethod is still able to capture around 95%. For the larger basissets, the percentage recovery is even better.In absolute terms, the error introduced by using HFPCto calculate the CPC can be as large as 0.75 kJ/mol forthe 6-31G(d) basis. The errors for formamide dimer, indole-benzene and adenine-thymine are all large. For 6-31+G(d)the maximum error reduces to 0.2 kJ/mol and, given that thecounterpoise procedure is only an estimate of the BSSE, thisdiscrepancy will be tolerable in many chemical applications.Columns 8–13 of Table I show analogous results usingDFPC with the B3LYP functional. The DFPC and HFPC re-sults show similar trends but there are two significant differ-ences. First, the B3LYP CPCs are significantly larger thanthose from HF. This observation was also made by Garzaet al.26 who attributed it to the tendency of DFT orbitals tobe more diffuse than HF orbitals. Second, our perturbativeapproach recovers a greater percentage of the B3LYP CPCsthan of the HF CPCs. Even for 6-31G(d), DFPC recovers atleast 98% of the CPC and, for the larger basis sets, the re-covery is usually >99.9%. This is consistent with our earlierresults9 and can be attributed to the fact that DFT energiestend to converge slightly faster than HF energies with respectto basis set size.C. EfficiencyCalculating the exact CPC by Eq. (2) requires two addi-tional calculations that use all of the basis functions in the sys-tem. This is expensive, typically taking almost twice as longas the calculation of the uncorrected interaction energy (1).Our perturbative schemes are much cheaper because the pri-mary density matrix is generated as part of the interactionenergy calculation, and therefore the computational overheadassociated with our perturbative schemes amounts to a singleFock build and energy evaluation in the secondary basis.The final columns of Table I show timing comparisonsbetween CPC calculations via the conventional and pertur-bative approaches. The time taken to calculate the interac-tion energy, T int, is simply the time required to compute thefour energy terms in Eq. (1) but, because we can reuse someof these terms, the overhead of calculating Eq. (2) is onlyT CPC = tEA(A∪B) + tEB(A∪B). Columns 14–19 show the timingratios T CPC/T int where, for example, a ratio of 1 implies thatDownloaded 26 Sep 2011 to 150.203.35.38. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions081105-3 Efficient CPC calculations J. Chem. Phys. 135, 081105 (2011)TABLE I. For the systems in the S22 dataset: (Columns 2–7) counterpoise corrections (kJ/mol) at the HF level for various basis sets and the percentagerecovered by HFPC; (Columns 8–13) counterpoise corrections (kJ/mol) at the B3LYP level for various basis sets and the percentage recovered by DFPC; and(Columns 14–19) timing ratios for conventional and perturbative counterpoise calculations using the 6-31+G(2df, p) basis.HF B3LYP T CPC/T int6-31G(d) 6-31+G(d) +G(2df, p) 6-31G(d) 6-31+G(d) +G(2df, p) Exact Perturbative Pert./exactComplex CPC % CPC % CPC % CPC % CPC % CPC % HF B3LYP HF B3LYP HF B3LYPAmmonia dimer (C2h) 4.29 92.1 1.15 95.7 0.58 96.6 6.82 98.4 1.09 99.1 0.65 100.0 2.24 2.88 0.59 0.37 0.26 0.13Water dimer (Cs) 3.90 93.6 4.48 95.5 1.72 96.5 6.80 98.4 4.57 99.3 1.63 99.4 1.34 1.49 0.38 0.21 0.28 0.14Formic acid dimer (C2h) 11.06 94.4 7.30 97.7 2.92 97.9 17.53 98.5 5.45 99.8 2.07 100.0 2.18 1.97 0.38 0.29 0.17 0.15Formamide dimer (C2h) 10.68 93.0 4.09 97.8 1.86 98.4 17.36 98.0 3.30 100.0 1.83 100.0 2.63 2.24 0.33 0.25 0.13 0.11Uracil dimer (C2h) 9.16 94.1 5.18 98.1 2.40 98.3 14.82 98.3 3.95 100.0 2.39 100.0 1.59 2.19 0.25 0.24 0.16 0.112-Pyridoxine 9.48 93.5 4.94 98.0 1.99 98.0 15.15 98.2 3.68 100.0 2.34 100.0 1.22 1.08 0.15 0.12 0.13 0.112-aminopyridine (C1)Adenine thymine WC (C1) 11.15 94.0 5.59 98.2 2.11 100.0 17.97 98.4 4.25 100.0 2.63 100.0 1.23 1.03 0.17 0.09 0.14 0.09Methane dimer (D3d) 0.44 95.5 0.14 100.0 0.11 100.0 0.95 98.9 0.12 100.0 0.09 100.0 2.94 2.65 0.64 0.36 0.22 0.14Ethene dimer (D2d) 3.33 93.7 0.68 97.1 0.50 96.0 3.84 99.0 0.70 100.0 0.67 100.0 2.93 2.62 0.72 0.42 0.24 0.16Benzene methane (C3) 2.05 95.6 0.57 98.2 0.48 97.9 2.64 99.2 0.68 100.0 0.64 100.0 1.89 0.82 0.34 0.14 0.18 0.17Benzene dimer (C2h) 9.34 93.5 1.67 98.2 1.52 96.7 9.57 99.1 2.49 100.0 2.50 100.0 2.35 3.74 0.56 0.41 0.24 0.11Pyrazine dimer (Cs) 8.21 94.4 3.24 97.8 1.91 97.9 8.38 99.2 2.25 100.0 1.87 100.0 1.74 1.60 0.25 0.22 0.15 0.14Uracil dimer stack (C2) 11.92 95.6 6.81 97.9 4.54 98.0 14.75 99.3 5.52 100.0 4.39 100.0 2.00 2.21 0.37 0.24 0.19 0.11Indole benzene stack (C1) 11.82 93.7 3.02 98.0 2.56 98.0 11.52 98.9 2.80 100.0 2.86 100.0 1.13 1.03 0.21 0.12 0.18 0.12Adenine thymine stack (C1) 15.36 95.1 8.58 98.0 4.73 98.1 18.06 99.1 5.37 100.0 4.08 100.0 1.43 1.80 0.17 0.14 0.12 0.08Ethene ethyne (C2v) 2.67 92.1 0.52 98.1 0.38 97.4 2.80 98.6 0.81 100.0 0.64 100.0 1.57 1.27 0.36 0.18 0.23 0.14Benzene water (Cs) 3.62 96.1 2.49 96.8 1.27 96.9 4.61 99.3 2.56 99.6 1.49 99.3 0.81 0.70 0.14 0.11 0.17 0.15Benzene ammonia (Cs) 2.71 95.9 1.66 97.0 0.85 96.5 3.57 99.4 1.78 99.4 1.16 100.0 0.96 1.09 0.22 0.17 0.23 0.15Benzene HCN (Cs) 2.93 96.2 1.16 98.3 0.83 97.6 3.29 99.4 1.27 100.0 1.19 100.0 1.17 1.51 0.23 0.22 0.20 0.15Benzene dimer (C2v) 3.68 95.4 1.18 98.3 1.17 98.3 3.69 99.2 1.43 100.0 1.59 100.0 1.72 2.11 0.35 0.32 0.20 0.15Indole benzene T-shape (C1) 6.08 95.2 2.33 97.9 1.77 98.3 6.66 99.2 2.55 100.0 2.36 100.0 1.02 1.05 0.17 0.14 0.17 0.13Phenol dimer (C1) 7.52 95.3 5.31 97.4 2.52 97.2 10.86 99.0 5.28 99.8 3.02 100.0 0.89 1.19 0.14 0.17 0.16 0.14Minimum 92.1 95.5 96.0 98.0 99.1 99.3 0.81 0.70 0.14 0.09 0.12 0.08Mean 94.5 97.7 97.8 98.9 99.9 99.9 1.68 1.74 0.32 0.22 0.19 0.13Maximum 96.2 100.0 100.0 99.4 100.0 100.0 2.94 3.74 0.72 0.42 0.28 0.17the CPC is as expensive to obtain as all the energy terms inEq. (1).Interestingly, not only does the DFPC approach capturea slightly higher fraction of the CPC, it does so at a lowerrelative cost compared to HFPC. This is primarily becausethe DFT calculations employing ghost functions take longerto converge making the SCF calculation of EA(A ∪ B) andEB(A ∪ B) more expensive. Our perturbative scheme effec-tively includes only one of these SCF cycles and is thereforenot plagued by this slower convergence.IV. CONCLUDING REMARKSComputational studies of weakly bound systems are im-portant, especially in biological contexts and, in order tocalculate accurate interaction energies in these systems, theBSSE must be corrected. This is usually done via the CPCand, in the past, this has necessitated calculations that are sig-nificantly more expensive than calculating the original inter-action energy. In this work, we have applied the HFPC andDFPC perturbative schemes to estimate the CPC at between1/4 and 1/12 of the cost of a conventional CPC calculation.Tests on the S22 set of weakly bounded complexes indicatethat these schemes recover almost all of the CPC and oftenintroduce no significant error at all.Recently, other forms of CPCs, mainly dealing with in-tramolecular BSSEs, have been proposed.27, 28 These proce-dures involve dividing the molecule into many fragments andcalculate each using a large fraction of the full basis. Our ap-proach applies immediately to such calculations and offersmajor computational savings.We conclude that the HFPC and DFPC schemes are ef-fective tools for evaluating counterpoise corrections for bothHF and DFT calculations.ACKNOWLEDGMENTSP.M.W.G. thanks the Australian Research Council forfunding (Grant Nos. DP0984806 and DP1094170) and APACfor a generous allocation of supercomputer resources. J.D.thanks the ANU/RSC for a PhD scholarship.1We note that this definition of the interaction energy does not include therelaxation energy associated with structural changes to the fragments uponcomplexation.2H. B. Jansen and P. Ros, Chem. Phys. Lett. 3, 140 (1969).3B. Liu and A. D. McLean, J. Chem. Phys. 59, 4557 (1973).4I. Mayer, Int. J. Quantum. Chem. 23, 341 (1983).5B. Jeziorski, B. Moszynski, and K. Szalewicz, Chem. Rev. 94, 1887 (1994).6S. F. Boys and F. Bernardi, Mol. Phys. 19, 553 (1970).7J. Deng, A. T. B. Gilbert, and P. M. W. Gill, J. Chem. Phys. 130, 231101(2009).Downloaded 26 Sep 2011 to 150.203.35.38. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions081105-4 Deng, Gilbert, and Gill J. Chem. Phys. 135, 081105 (2011)8J. Deng, A. T. B. Gilbert, and P. M. W. Gill, J. Chem. Phys. 133, 044116(2010).9J. Deng, A. T. B. Gilbert, and P. M. W. Gill, Phys. Chem. Chem. Phys. 12,10759 (2010).10J. Deng and P. M. W. Gill, J. Chem. Phys. 134, 081103 (2011).11K. Wolinski and P. Pulay, J. Chem. Phys. 118, 9497 (2003).12W. Liang and M. Head-Gordon, J. Phys. Chem. A. 108, 3206 (2004).13R. P. Steele, R. A. DiStasio, Jr., Y. Shao, J. Kong, and M. Head-Gordon, J.Chem. Phys. 125, 074108 (2006).14R. P. Steele, Y. Shao, R. A. DiStasio, Jr., and M. Head-Gordon, J. Phys.Chem. A. 110, 13915 (2006).15R. P. Steele and M. Head-Gordon, J. Chem. Phys. 124, 184108 (2006).16R. P. Steele, R. A. DiStasio, Jr., and M. Head-Gordon, J. Chem. TheoryComput. 5, 1560 (2009).17L. Kong and E. F. Valeev, J. Chem. Phys. 133, 174126 (2010).18S. Grimme, J. Comput. Chem. 25, 1463 (2004).19S. Grimme, J. Comput. Chem. 27, 1787 (2006).20S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132,154104 (2010).21P. Jurecka, J. Sponer, J. Cerny, and P. Hobza, Phys. Chem. Chem. Phys. 8,1985 (2006).22P. M. W. Gill, B. G. Johnson, and J. A. Pople, Chem. Phys. Lett. 209, 506(1993).23Y. Shao, L. Fusti-Molnar, Y. Jung, J. Kussmann, C. Ochsenfeld,S. T. Brown, A. T. B. Gilbert, L. V. Slipchenko, S. V. Levchenko,D. P. O’Neill, R. A. DiStasio, Jr., R. C. Lochan, T. Wang, G. J. O.Beran, N. A. Besley, J. M. Herbert, C. Y. Lin, T. V. Voorhis, S. H. Chien,A. Sodt, R. P. Steele, V. A. Rassolov, P. E. Maslen, P. P. Korambath,R. D. Adamson, B. Austin, J. Baker, E. F. C. Byrd, H. Dachsel, R. J.Doerksen, A. Dreuw, B. D. Dunietz, A. D. Dutoi, T. R. Furlani, S. R.Gwaltney, A. Heyden, S. Hirata, C. P. Hsu, G. Kedziora, R. Z. Khalliulin,P. Klunzinger, A. M. Lee, M. S. Lee, W. Z. Liang, I. Lotan, N. Nair, B.Peters, E. I. Proynov, P. A. Pieniazek, Y. M. Rhee, J. Ritchie, E. Rosta,C. D. Sherrill, A. C. Simmonett, J. E. Subotnik, H. L. Woodcock III,W. Zhang, A. T. Bell, A. K. Chakraborty, D. M. Chipman, F. J. Keil,A. Warshel, W. J. Hehre, H. F. Schaefer III, J. Kong, A. I. Krylov,P. M.W. Gill, and M. Head-Gordon, Phys. Chem. Chem. Phys. 8, 3172(2006).24See supplementary material at http://dx.doi.org/10.1063/1.3632054 for rawenergies.25J. Antony and S. Grimme, Phys. Chem. Chem. Phys. 8, 5287 (2006).26J. Garza, J. Ramirez, and R. Vargas, J. Phys. Chem. A 109, 643(2005).27F. Jensen, J. Chem. Theory Comput. 6, 100 (2010).28A. Galano and J. R. Alvarez-Idaboy, J. Comput. Chem. 27, 1203(2006).Downloaded 26 Sep 2011 to 150.203.35.38. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions

Communication: Efficient counterpoise corrections by a perturbative approach

Abstract

Similar works

Full text

Available Versions

The Australian National University