Report

A Density-Functional Theory for Covalent and Noncovalent Chemistry Non-empirical and fast • Review of the XDM (exchange-hole dipole moment) dispersion model of Becke and Johnson • Combining the model with non-empirical exchange-correlation GGAs (Kannemann and Becke) • Tests on standard bio-organic benchmark sets Dispersion Interactions from the Exchange-Hole Dipole Moment: The “XDM” model Axel D. Becke and Erin R. Johnson* Email: [email protected] Department of Chemistry, Dalhousie University, Halifax, Nova Scotia, Canada (*now at University of California, Merced) What is the source of the “instantaneous” multipole moments that generate the dispersion interaction? Suggestion: Becke and Johnson, J. Chem. Phys. 122, 154104 (2005) > Becke and Johnson, J. Chem. Phys. 127, 154108 (2007) < The dipole moment of the exchange hole! (Position, rather than time, dependent) In Hartree-Fock theory, the total energy of many-electron system is given by: EHF 1 1 (r1 ) (r2 ) 3 3 2 3 3 i i d r Vnuc d r d r1d r2 E X 2 i , 2 r12 EX EX EX E X N i2 i 1 i (r1 ) j (r1 ) i (r2 ) j (r2 ) 3 3 1 d r1d r2 2 ij r12 The “exchange” energy The Exchange Hole The ( -spin) exchange energy can be rewritten as follows: E X where hX (r1 , r2 ) 3 3 1 (r1 ) d r2 d r1 2 r12 hX (r1 , r2 ) 1 i (r1 ) j (r1 ) i (r2 ) j (r2 ) (r1 ) ij is called the “exchange hole”. Physical interpretation: each electron interacts with a “hole” whose shape (in terms of r2) depends on the electron’s position r1. When an electron is at r1, the hole measures the depletion of probability, with respect to the total electron density, of finding another electron of the same spin at r2. This arises from exchange antisymmetry. It is simple to prove that… • The hole is always negative. • The probability of finding another same-spin electron at r2=r1 (“on top” of the reference electron) is completely extinguished: hX (r1 , r1 ) (r1 ) Pauli or exchange “repulsion”! • The hole always contains exactly (minus) one electron: 3 h ( r , r ) d X 1 2 r2 1 Dispersion Model: the Basic Idea An electron plus its exchange hole always has zero total charge but in general a non-zero dipole moment! This r1-dependent dipole moment is easily obtained by integrating over r2: 1 d X (r1 ) rij i (r1 ) j (r1 ) r1 (r1 ) ij rij r2 i (r2 ) j (r2 )d 3 r2 Note that only occupied orbitals are involved. (This reduces, for the H atom, to the exact dipole moment of the H atom when the electron is at r1.) In a spherical atom, consider the following simplified “2-point” picture: e- (r , W) dX h (r - dX , W) Notice that this picture generates higher multipole moments as well (with respect to the nucleus as origin) given by Nucleus [r (r d X ) ] and that all these moments depend only on the magnitude of the exchange-hole dipole moment. This is significant because the magnitude dX(r) of the exchange-hole dipole moment can be approximated using local densities and the Becke-Roussel exchange-hole model [Phys. Rev. A 39, 3761 (1989)], a 2nd-order GGA: d X (r ) b where b is the displacement from the reference point of the mean position of the BR model hole [Becke and Johnson, J. Chem. Phys. 123, 154101 (2005)]. Therefore, the entire van der Waals theory that follows has two variants: Orbital (XX) based, or Density-functional (BR) based XX performs slightly better in rare-gas systems. BR performs better in intermolecular complexes. All our current work employs the BR (DFT) variant. The Dispersion Interaction: Spherical Atoms rA Vint(rA,rB) rB Vint(rA,rB) = multipole moments of electron+hole at rA interacting with multipole moments of electron+hole at rB 2nd-Order Ground-State Perturbation Theory in the Closure (Ünsold) Approximation If the first-order, ground-state energy correction arising from a perturbation Vpert is zero: E (1) V pert 0 Then the second-order correction is approximately given by V pert 2 E ( 2) E avg where the expectation values are in the ground state and Eavg is the average excitation energy. To evaluate the expectation value <Vint2> square the multipole-multipole interaction Vint(rA,rB): Vint2(rA,rB) = (dipole-dipole + dipole-quadrupole + dipole-octopole + quadrupole-quadrupole +…)2 Then integrate the squared interaction over all rA and rB. This is a “semiclassical” calculation of <Vint2>. The result is where E disp C6 C8 C10 6 8 10 R R R 2 12 A 12 B C6 3 Eavg 12 A 22 B 22 A 12 B C8 Eavg 4 12 A 32 B 32 A 12 B 14 22 A 22 B C10 3 Eavg 5 Eavg with atomic moment integrals given by 2 (r )[r (r d X ) ]2 d 3 r What about Eavg ? Assume that Eavg E A EB and that, for each atom, 2 12 E 3 where α is its dipole polarizability. This easily follows from the same “semiclassical” 2nd-order perturbation theory applied to the polarizability of each atom. Thus our C6, C8, C10’s depend on atomic polarizabilities and moment integrations from Hartree-Fock (or KS) calculations! No fitted parameters or explicitly correlated wavefunctions! How well does it work? On the 21 pairs of the atoms H, He, Ne, Ar, Kr, Xe, the mean absolute percent errors are: C6 3.4 % C8 21.5 % C10 21.5 % From Free Atoms to Atoms in Molecules • Partition a molecular system into “atoms” using Hirshfeld weight functions: iat (r ) wi (r ) nat (r ) n where at is a spherical free atomic density placed at the appropriate nucleus and the n summation is over all nuclei. wi (r ) has value close to 1 at points near nucleus i and close to 0 elsewhere. Also, w (r ) 1 i i • Assume that an intermolecular Cm (m=6,8,10) can be written as a sum of interatomic Cm,ij : A B i j Cm Cm,ij i in A j in B In the previous expressions for Cm replace A and B with i and j : C 6,ij C8,ij C10,ij 2 2 2 1 i 1 j 3 (E i E j ) 12 i 22 j 22 i 12 j (Ei E j ) 2 2 2 2 2 2 4 1 i 3 j 3 i 1 j 14 2 i 2 j 3 (Ei E j ) 5 (Ei E j ) with 2 i generalized to 2 i wi (r ) (r )[r (r d X ) ]2 d 3 r and, where i 2 12 i Ei 3 i is the effective polarizability of atom i in A. We propose that r 3 i i i , free 3 r i , free r 3 i r 3 wi ( r ) ( r )d 3 r r 3 i , free r 3 i , free (r )d 3 r These are effective volume integrations. This is motivated by the well known qualitative (if not quantitative) general relationship between polarizability and volume. See Kannemann and Becke, JCP 136, 034109 (2012). All radii r in the above integrals and in the integrals for 2 i are with respect to the position of nucleus i. Atom-Molecule and Molecule-Molecule Dispersion Coefficients Test set: H2 and N2 with He, Ne, Ar, Kr, and Xe. Cl2 with He, Ne, Ar, Kr, and Xe (except C10). H2-H2 H2-N2 N2-N2 (only H2-H2 for C10). Fully-numerical Hartree-Fock calculations on the monomers using the NUMOL program (Becke and Dickson, 1989). MAPEs with respect to dispersion coefficients from frequency dependent MBPT polarizabilities: 12.7% for C6 16.5% for C8 11.9% for C10 (On a much more extensive test set of 178 intermolecular C6’s, the model has a MAPE of 9.1%) Everything, so far, has been about the asymptotic dispersion series between atom pairs, E disp C6 C8 C10 6 8 10 R R R The asymptotic series needs to be damped in order to avoid divergences when R is small. i.e. need information about characteristic R values inside of which the asymptotic series is no longer valid. The usual approach is to use empirical vdW radii. However… “Critical” Interatomic Separation Since we can compute C6, C8, and C10 non-empirically, we can obtain non-empirical range information. There is a “critical” Rc,ij where the three dispersion terms are approximately equal: Take Rc,ij as the average of , , and The asymptotic dispersion series is obviously meaningless inside Rc,ij. Therefore … General Dispersion Energy Formula … we use Rc,ij to damp the dispersion energy at small internuclear separations as follows: Edisp C6,ij C8,ij C10,ij 6 8 10 6 8 RvdW ,ij Rij RvdW ,ij Rij10 i j RvdW ,ij Rij RvdW ,ij a1 Rc,ij a2 where Rvdw,ij is an effective van der Waals separation with only two universal fit parameters. Best-fit a1 and a2 values depend on the exchange-correlation theories with which the above is combined. … which brings us to the second part of the talk … Key XDM Publications A. D. Becke and E. R. Johnson, J. Chem. Phys. 127, 154108 (2007) Dispersion coefficients E. R. Johnson and A. D. Becke, J. Chem. Phys. 124, 174104 (2006) Damping functions A. D. Becke and E. R. Johnson, J. Chem. Phys. 123, 154101 (2005) Transform to a DFT Exchange-Correlation GGAs + XDM Dispersion Axel D. Becke and Felix O. Kannemann Email: [email protected] Department of Chemistry, Dalhousie University, Halifax, Nova Scotia, Canada It would be nice if we could use “standard” XC-GGAs (instead of, eg., Hartree-Fock) plus the XDM dispersion model to treat vdW interactions. Let’s look at the exchange part first. Exchange GGA Functionals B86, B86b, B88 E XB86 E XLDA 4 / 3 2 0.0036 1 0.0042 E XB 86 b E XLDA 0.0036 E B88 X E LDA X 4 / 3 2 1 0.006 2 4/5 4 / 3 2 0.0042 1 0.0252 arcsin h 1/ 3 where the local (spin) density part is E XLDA 3 3 2 4 and the “reduced” (spin) density gradient is 4/3 4 / 3 4 / 3 Exchange GGA Functionals PW86, PW91, PBE(96) E X e XLDA ( ) g X ( s) g g PW 91 X PW 86 X e where LDA X 3 3 1/ 3 4 / 3 ( ) ( ) 4 (s) 1 1.296s 14s 0.2s 2 4 6 1 / 15 1 0.19645s sinh 1 (7.7956s) [0.2743 0.1508exp(100s 2 )]s 2 ( s) 1 0.19645s sinh 1 (7.7956s) 0.004s 4 g PBE X 0.804 ( s ) 1 0.804 1 0.21951 s 2 / 0.804 The “reduced” (spin) density gradient is and, for spin polarized systems, E X ( , ) s 2(3 2 )1 / 3 4 / 3 1 E X (2 ) E X (2 ) 2 revPBE (Yang group) g XPBE ( s ) 1 0.804 0.804 1 0.21951 s 2 / 0.804 g XrevPBE ( s) 1 1.245 1.245 1 0.21951 s 2 / 1.245 Now plot the exchange-only interaction energy in Ne2 Note that Hartree-Fock (exact exchange) is repulsive! Plot is from Kannemann and Becke, JCTC 5, 719 (2009) Ne2 Can get anything! from artifactual binding to massive over-repulsion, depending on the choice of functional! Which exchange GGA best reproduces exact Hartree-Fock? Lacks and Gordon, PRA 47, 4681 (1993) Kannemann and Becke, JCTC 5, 719 (2009) Murray, Lee, and Langreth, JCTC 5, 2754 (2009) PW86, followed by B86b Exchange Enhancement Factor [Zhang, Pan, Yang, JCP 107, 7921 (1997)] PW86 is a completely non-empirical exchange functional! Its 4 parameters are fit to a theoretical exchange-hole model. Perdew and Wang, PRB 33, 8800 (1986) E X e XLDA ( ) g X ( s) g PW 86 X (s) 1 1.296s 14s 0.2s 2 s 4 6 1 / 15 2(3 2 )1 / 3 4 / 3 That PW86 accurately reproduces Hartree-Fock repulsion energies is remarkable. The underlying theoretical model (truncated “GEA” hole) knows nothing about closed-shell atomic or molecular interactions! Could be a fortuitous accident?! Nevertheless, no parameters need to be fit to data Now, what about (dynamical) correlation? Use the non-empirical “PBE” correlation functional: Perdew, Burke, and Ernzerhof, PRL 77, 3865 (1996) Therefore we have, EXC E E XDM disp PW 86 X E PBE C E XDM disp C6,ij C8,ij C10 ,ij 1 6 8 10 6 8 10 2 i j RvdW ,ij Rij RvdW ,ij Rij RvdW ,ij Rij RvdW ,ij a1 Rc,ij a2 with a1 and a2 to be determined. How to determine a1 and a2? Fit to the binding energies of the prototypical dispersion-bound rare-gas systems He2 HeNe HeAr Ne2 NeAr Ar2 (reference data from Tang and Toennies, JCP 118, 4976 (2003)) At the CBS limit, we find the best-fit values a1=0.65 a2=1.68 (Kannemann and Becke, to be published) Gaussian 09, aug-cc-pV5Z, counterpoise, ultrafine grid, BEs in microHartree a1 = 0.65 a2 = 1.68Å RMS%E = 4.2% Note that in all subsequent benchmarking there is no (re)fitting of parameters! Our functional is, from here on, essentially nonempirical in all its parts. GGA + XDM Publications F.O. Kannemann and A.D. Becke, JCTC 5, 719 (2009) Rare-gas diatomics (numerical post-LDA) F.O. Kannemann and A. D. Becke, JCTC 6, 1081 (2010) Intermolecular complexes (numerical post-LDA) A.D. Becke, A.A. Arabi and F.O.Kannemann, Can. J. Chem. 88, 1057 (2010) Dunning aDZ and aTZ basis-set calculations (post-G09) Tests on Standard BioOrganic Benchmark Sets Axel D. Becke and Felix O. Kannemann Email: [email protected] Department of Chemistry, Dalhousie University, Halifax, Nova Scotia, Canada The “S22” and “S66” vdW Benchmark Sets of Hobza et al S22 S66 hydrogen bonding dispersion other noncovalent interactions References S22: Jurecka, Sponer, Cerny, Hobza, PCCP 8, 1985 (2006) S66(x8): Rezac, Riley, Hobza JCTC 7, 2427 (2011) The next slides contain Mean Absolute Deviations (MADs) for the S22 and the S66 benchmark sets in comparison with other popular DFT methods. Data for all methods other than ours are from Goerigk, Kruse & Grimme, CPC 12, 3421 (2011) All computations employ the def2-QZVP basis set PW86PBE-XDM no CP (0.30) CP (0.28) S66 PW86PBE-XDM no CP (0.27) CP (0.23) Aside Can we do better by combining XDM with hybrid functionals? (rather than the pure GGA, PW86+PBE) Burns, Vazquez-Mayagoitia, Sumpter, Sherrill, JCP 134, 084107 (2011) S22 MAD (kcal/mol) for B3LYP-XDM and other DFT methods However • The pure GGA, PW86+PBE, is completely nonempirical (B3LYP, with 3 fitted parameters, is not) • Density-fit basis sets can speed up pure GGA calculations, with no loss of accuracy, by an order of magnitude! Tests on the S66x8 Set (Nonequilibrium Geometries) Axel D. Becke and Felix O. Kannemann Email: [email protected] Department of Chemistry, Dalhousie University, Halifax, Nova Scotia, Canada S66x8 A benchmark set of 66 vdW complexes of bio-organic interest, at 8 intermonomer separations: 0.90, 0.95, 1.00, 1.05, 1.10, 1.25, 1.50, 2.00 (relative to the equilibrium intermonomer separation) Rezac, Riley, Hobza JCTC 7, 2427 (2011) Important because complex systems and materials may contain many vdW interactions between groups at nonequilibrium (especially stretched) geometries! Mean Percent Errors (MPEs) versus geometry 10 0 MPE [%] -10 -20 PW86PBE+XDM -30 M06-2X M05-2X -40 M06L -50 0.90 1.00 1.10 1.25 distance multiplier 1.50 2.00 25 20 15 MPE [%] 10 5 0 PW86PBE+XDM -5 M06-2X-D3 -10 M05-2X-D3 -15 M06L-D3 -20 0.90 1.00 1.10 1.25 distance multiplier 1.50 2.00 25 20 15 MPE [%] 10 5 0 PW86PBE+XDM -5 PBE-D3 -10 BLYP-D3 -15 B97-D3 -20 0.90 1.00 1.10 1.25 distance multiplier 1.50 2.00 10 5 MPE [%] 0 -5 PW86PBE+XDM LC-ωPBE-D3 -10 B3LYP-D3 -15 B2-PLYP-D3 -20 0.90 1.00 1.10 1.25 distance multiplier 1.50 2.00 Mean Absolute Percent Errors (MAPEs) versus geometry 80 PW86PBE+XDM 70 M06-2X 60 M05-2X MAPE [%] 50 M06L 40 30 20 10 0 0.90 1.00 1.10 1.25 distance multiplier 1.50 2.00 35 PW86PBE+XDM 30 M06-2X-D3 MAPE [%] 25 M05-2X-D3 M06L-D3 20 15 10 5 0 0.90 1.00 1.10 1.25 distance multiplier 1.50 2.00 35 PW86PBE+XDM 30 PBE-D3 25 BLYP-D3 MAPE [%] B97-D3 20 15 10 5 0 0.90 1.00 1.10 1.25 distance multiplier 1.50 2.00 35 PW86PBE+XDM 30 LC-ωPBE-D3 MAPE [%] 25 B3LYP-D3 B2-PLYP-D3 20 15 10 5 0 0.90 1.00 1.10 1.25 distance multiplier 1.50 2.00 pdf file of all S66x8 curves (Notice the importance of the dispersion term!) Interaction types in S66x8: • H-bonding 23 • Dispersion 23 • Mixed (“other”) 20 Dispersion is further divided into: pi-pi (10), aliphatic-aliphatic (5), and pi-aliphatic (8) How good is PW86+PBE+XDM for ordinary thermochemistry? Consider the functional XDM EXC EXGGA ECPBE Edisp with a variety of standard exchange GGAs in the first term. On the “G3/99” benchmark set of 222 atomization energies of organic/inorganic molecules (Curtiss, Raghavachari, Pople) we obtain the following error statistics, in kcal/mol: For standard hybrid functionals (eg., B3LYP) the MAE is of order 5-6 kcal/mol. The very best DFTs have MAE as small as 2-3 kcal/mol, but with fitted params! Availability of XDM code • Has been implemented (B3LYP-XDM) in Q-Chem by Kong, Gan, Proynov, Freindorf, and Furlani, PRA 79, 042510 (2009) • A “post-Gaussian09” code will be available from us by the end of 2012. Uses G09 to perform the PW86+PBE part, then adds XDM perturbatively. Can do Berny geometry optimizations using the EXTERNAL keyword! Very fast with density-fit basis sets! Many thanks to: Natural Sciences and Engineering Research Council of Canada the Killam Trust of Dalhousie University (Killam Chair) ACEnet (the Atlantic Computational Excellence Network)