1 / 82

Bartolomeo Civalleri Theoretical Chemistry Group Department of Chemistry IFM

B3LYP augmented with an empirical dispersion term (B3LYP-D*) as applied to solids. Bartolomeo Civalleri Theoretical Chemistry Group Department of Chemistry IFM & NIS Centre of Excellence University of Torino bartolomeo.civalleri@unito.it. Weak interactions in crystalline solids.

lfloyd
Download Presentation

Bartolomeo Civalleri Theoretical Chemistry Group Department of Chemistry IFM

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. B3LYP augmented with an empirical dispersion term (B3LYP-D*) as applied to solids Bartolomeo Civalleri Theoretical Chemistry Group Department of Chemistry IFM & NIS Centre of Excellence University of Torino bartolomeo.civalleri@unito.it Vallico Sotto July 2009

  2. Weak interactions in crystalline solids • Cohesive forces • long-range: electrostatic, induction, dispersion • short-range: exchange repulsion, charge transfer • Weak interactions play an important role in the solid state (see T. Steiner, Angew. Chem. Int. Ed.41 (2002) 48) • Molecular recognition  crystal packing • Supramolecular chemistry and crystalline engineering • Molecular crystals (polymorphism) • Layered and composite/intercalated materials • Adsorption and reactivity on surfaces • Very important for many properties of interest: structure, interaction energies, vibrational frequencies and thermodynamics, elastic constants, relative stability, … Vallico Sotto July 2009

  3. State of the art in ab initio calculations of MCs The “standard” ingredients (with many variants) : a) HF or DFT (Kohn-Sham) Hamiltonians b) Plane-Wave + Pseudopotentials(no BSSE) or Localized functions (Gaussian) + All-Electron (BSSE) [CRYSTAL06] c) Analytic derivatives of energy and other observables (e.g. phonons) DFT is the most common choice to include electron correlation Only LDA, GGA and hybrid-GGA (e.g. B3LYP) methods are routinely available in solid state codes Present DFT functionals do not account for dispersion energy BSSE can give binding where there is none. Vallico Sotto July 2009

  4. Hydrogen bonded molecular crystals: B3LYP results 3D, 2D and 1D HB molecular crystals Volume in Å3. Energy in kJ/mol. Deviation from experimental volumes in parentheses. • BSSE correctedcohesiveenergies are independent from the adopted basis set but they are markedly underestimated • With the large TZP basis set, BSSE is very small but predicted unit cell sizes are largely overestimated • With small basis sets, BSSE artificially compensates the missing dispersion energy • When HB is dominating, B3LYP gives good lattice parameters (not shown) Vallico Sotto July 2009

  5. DFT vs vdW forces: new hopes… • How to deal with dispersion interactions in DFT? • New functionals: • vdW-DFT (Langreth, Lundqvist and co-workers) • beyond m-GGA (Perdew’s “Jacob’s Ladder” – fifth rung) • screened Coulomb (or CAM) functionals (Scuseria, Handy, Savin, Hirao, …)- Truhlar’s family (M05, M05-X, M06, …) • Perturbational electron-interaction corrections - on top of range-separated hybrids (Savin and coworkers: e.g. Goll et al. PCCP (2005)) • - double-hybrids (Grimme JCP 124 (2006) 034108, Head-Gordon JPC-A (2008) ) • A pragmatic approach: Wilson-Levy (WL) correlation functional (a-posteriori HF) (T.A. Walsh, PCCP 7 (2005) 403; B. Civalleri et al. CPL 451 (2008) 287) • Empirical corrections by adding a -C6/R6 term (Grimme, Neumann, Yang, Zimmerli, Scoles, …) • A DF model of the dispersion interaction: C6 in terms of exchange-hole dipole moment (Becke-Johnson, JCP 123 (2005) 024101, JCP 124 (2006) 014104, JCP 127 (2007) 154108)or C6 from MLWFs (Silvestrelli, PRL 100 (2008) 053002) • Dispersion-corrected atom centered pseudo-potentials (U. Rothlisberger and co-workers: e.g. Tapavicza et al. JCTC (2007), G. DiLabio CPL (2008)) Vallico Sotto July 2009

  6. Empirical –C6/R6 correction: Grimme’s model S. Grimme, J. Comput. Chem., 2004, 25, 1463 and J. Comput. Chem., 2006, 27, 1787 Atom-atom additive damped empirical potential of the form -f(R)C6/R6 where • s6: scaling factor for each DFT method (s6=1.05 for B3LYP) • C6ij are computed from atomic dispersion coefficients: C6ij = C6i·C6j • Rvdw is the sum of atomic van der Waals radii: Rvdw=Rivdw+Rjvdw • d determines the steepness of the damping function (d=20) • summation over g truncated at 25 Å (estimated error < 0.02 kJ/mol on DE) • Grimme proposed a set of parameters (i.e. C6i and Rivdw) from H to Xe Total energy is then computed as: Implemented in CRYSTAL06 for energy and gradients (atoms and cell): B. Civalleri, C.M. Zicovich-Wilson, et al., CrystEngComm, 2008, DOI: 10.1039/b715018k (see supplementary material for erratum) Vallico Sotto July 2009

  7. GRIMME input block E.g.: Urea – B3LYP-D Urea CRYSTAL 0 0 0 113 5.565 4.684 5 6 0.0000 0.5000 0.3260 8 0.0000 0.5000 0.5953 7 0.1459 0.6459 0.1766 1 0.2575 0.7575 0.2827 1 0.1441 0.6441-0.0380 Optional keywords END (ENDG) Basis set END … GRIMME 1.05 20. 25. 4 1 0.14 1.001 6 1.75 1.452 7 1.23 1.397 8 0.70 1.342 … END Grimme empirical dispersion keyword s6 (scaling factor) d (steepness) Rcut (cut-off radius, Å) Nr. of atomic species Atomic number C6 (Jnm6 mol−1)Rvdw (Å) Atomic number C6 Rvdw Atomic number C6 Rvdw Atomic number C6 Rvdw End of SCF&method input section Vallico Sotto July 2009

  8. GRIMMEdataset Parameters available from H to Xe Rvdwvalues arederived from the radius of the 0.01a0−3 electron density contour fromROHF/TZV computations of the atoms in the ground state C6 coefficients derivedfrom the London formula for dispersion. DFT/PBE0calculations of atomic ionization potentials Ip and static dipole polarizabilitiesα. The C6 coefficient for atom i (in Jnm6 mol−1) is thengiven as (Ip and α in atomic units) C6i = 0.05NIpi αi where N has values 2, 10, 18, 36, and 54 for atoms from rows1–5 of the periodic table Suitable for solids? S. Grimme, J. Comput. Chem., 2006, 27, 1787 Vallico Sotto July 2009

  9. Propane Boric acid C2H2 CO2 NH3 Formic acid 1,4-dichloro-benzene Urea Naphthalene Formamide C6H6 Succinic anhydride 1,4-dicyano-benzene Urotropine Tests on a set of selected molecular crystals 14 molecular crystals both dispersion bonded and hydrogen bonded • Experimental sublimation energies at 298K available from published data (estimated error bar: ±4 kJ/mol) • For some of them accurate low temperature structural data from NPD Vallico Sotto July 2009

  10. Cohesive energies: B3LYP vs B3LYP-D Grimme BSSE corrected cohesive energies vs Experimental data Exp.: -DE=DH0sub(T)+2RT from data at 298K Cell fixed geometry optimization of the atomic positions at B3LYP/6-31G(d,p) • B3LYP: MD=54.4 • Empirical correction definitely improves cohesive energies • Tendency of B3LYP-D Grimme to overestimate cohesive energy(MD=-6.0 & MAD=8.9) especially for HB molecular crystals • Small basis sets suffer from large BSSE • BSSE corrected data are less basis set dependent -110 < DE < -25 kJ/mol Vallico Sotto July 2009

  11. Grimme’s model: the role of the damping function • The damping function is needed: • to avoid near singularities for small interatomic distances • some short-range correlation effects are already contained in the density functional • However: • crystal packing leads to larger overlap between molecular charge densities • damping function must act to longer-range where the B3LYP method does not contribute to the intermolecular interactions • atomic vdW radii define where the –f(R)C6/R6 contribution becomes dominant • atomic vdW radius for H very important • Strategy: scaling the atomic RvdW CarbonRvdw(C)=161 pm RvdW From: S. Grimme, J. Comput. Chem. 25 (2004) 1463 See also: P. Jurecka et al. J. Comput. Chem. 28 (2007) 555 Vallico Sotto July 2009

  12. Determination of the atomic vdW radii scaling factor • s6=1.00 • Atomic vdW radii (RvdW) were progressively increased to find the best agreement between computed and experimental data • larger scaling for the vdW radius of H (RH) • better balance between dispersion bonded and hydrogen bonded molecular crystals • SRvdW=1.05; SRH=1.30 MD:Mean Deviation; MAD:Mean Absolute Deviation; RMS:Root-Mean-Square Deviation from experiment (kJ/mol)J. S. Chickos and W. E. Acree, J. Phys. Chem. Ref. Data, 2002, 31, 537 B3LYP-D* Vallico Sotto July 2009

  13. Cohesive energies with B3LYP-D* BSSE corrected cohesive energies vs Experimental data Exp.: -DE=DH0sub(T)+2RT from data at 298K Cell fixed geometry optimization of the atomic positions at B3LYP/6-31G(d,p) • B3LYP-D* gives cohesive energies in excellent agreement with experimental data • MD=2.2 & MAD=6.3 • Better balance between hydrogen bonded and dispersion bonded molecular crystals Vallico Sotto July 2009

  14. Urotropine C6H6 C2H2 CO2 Urea NH3 Geometry optimization: B3LYP-D Grimme vs B3LYP-D* Lattice parameters • TZP basis set suffers from a remarkably small BSSE • B3LYP/TZP largely overestimates lattice parameters • B3LYP-D* lattice parameters are in excellent agreement with experimental data • B3LYP-D (Grimme) gives too short lattice constants TZP ______ 6-31G(d,p) - - - - Vallico Sotto July 2009

  15. Interlayer interaction in graphite: B3LYP-D* Exp.: a = 2.46 (fixed) c = 6.71 Å HF BS: 6-31G(d) hybrids BSSE corrected GGA __ LDA B3LYP-D* B3LYP-D* gives results in excellent agreement wrt experiment At long-range empirical correction correctly decays as -1/R4 Vallico Sotto July 2009

  16. Interlayer interaction in graphite: B3LYP-D* BS: 6-31G(d) Exp.: in Å a = 2.46 c = 6.70 Opt.: a = 2.453 c = 6.640 B3LYP-D* (CPC) Exp. B3LYP-D* __ B3LYP-D* gives results in excellent agreement wrt experiment Vallico Sotto July 2009

  17. Interlayer interaction in graphite: B3LYP-D* long-range 1/R4 At long-range empirical London-type formula correctly decays as -1/R4 Vallico Sotto July 2009

  18. CO adsorption on MgO(001) Distances in Å, interaction energies in kJ/mol, vibrational frequency shifts in cm-1 MgO basis set: Alhrichs’ TVZ; *MgO top-most layer: Alhrichs’ QZVP B3LYP-MP2 (slab): DE(CPC)=-12.2 kJ/molM06-HF (cluster): DE(CPC)=-24.6 kJ/mol; Dwh=22 cm-1CI (cluster): DE(CPC)=-10.5 kJ/mol; Dwh= 19 cm-1 Exp: R. Wichtendahl et al.Surf. Sci. 423, 90 (1999); G. Spoto et al. Prog. Surf. Sci. 76, 71 (2004) Vallico Sotto July 2009

  19. For molecular crystals: • Grimme’s scheme  Recalibration needed • Useful tool to correct the PES. Electron density is indirectly influenced • It gives results in excellent agreement wrt experiment for cohesive energies and structures. • Lattice modes also well reproduced • A large basis set should be adopted (e.g. TZP) to reduce the BSSE In perspective: • Work is in progress to test the transferability of B3LYP-D* to alkali halides (e.g. which C6 for Li+, Na+, …?) • C6 from non-empirical models (e.g. Becke-Johnson, Silvestrelli, …) Conclusions Dispersion interactions are crucial and must be taken into account Vallico Sotto July 2009

  20. Calculation of vibrational frequencies and tools for their analysis with CRYSTAL06R. Dovesi (Torino )L. Valenzano (Torino) C. Zicovich (Cuernavaca)Y. Noël (Paris)F. Pascale (Nancy) Vallico Sotto July 2009

  21. The CRYSTAL code: Quantum-Mechanical, ab-initio, periodic, using a local basis set (“Atomic Orbitals”)

  22. CRYSTA06 web site www.crystal.unito.it Vallico Sotto July 2009

  23. A few historical references Formulation and implementation (graphite) • C. Pisani and R. DovesiExact exchange Hartree-Fock calculations for periodic systems.I. Illustration of the method.Int. J. Quantum Chem. 17, 501-516 (1980). • R. Dovesi, C. Pisani and C. RoettiExact exchange Hartree-Fock calculations for periodic systems.II. Results for graphite and hexagonal boron nitrideInt. J. Quantum Chem 17, 517-529 (1980). The Coulomb problem: multipolar expansion+Ewald • R. Dovesi, C. Pisani, C. Roetti and V.R. SaundersTreatment of Coulomb interactions in Hartree-Fock calculations of periodic systems.Phys. Rev. B 28, 5781-5792 (1983). • V.R. Saunders, C. Freyria-Fava, R. Dovesi, L. Salasco and C. RoettiOn the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions. Mol. Phys. 77, 629-665 (1992) • V.R. Saunders, C. Freyria-Fava, R. Dovesi and C. RoettiOn the electrostaticpotential in linear periodic polymers. Comp. Phys. Comm. 84, 156-172 (1994) Towards the hybrids....... M. Causà, R. Dovesi, C. Pisani, R. Colle and A. FortunelliCorrelation correction to the Hartree-Fock total energy of solids.Phys. Rev. B 36, 891-897 (1987). Vallico Sotto July 2009

  24. Consistent treatment of Periodicity 3D - Crystalline solids (230 space groups) 2D - Films and surfaces (80 layer groups) 1D – Polymers(75 rod groups) 0D – Molecules(32 point groups) Infinite sums of particle interactions Ewald's method Specific formulæ for 1D, 2D, 3D Full exploitation of symmetry in direct space in reciprocal space The periodic model Vallico Sotto July 2009

  25. Exchange functionals Slater [L] von Barth-Hedin [L] Becke '88 [G] Perdew-Wang '91 [G] Perdew-Burke-Ernzerhof [G] Correlation functionals Vosko-Willk-Nusair (VWN5 parameterization) [L] Perdew-Wang [L] Perdew-Zunger '81 [L] von Barth-Hedin [L] Lee-Yang-Parr [G] Perdew '86 [G] Perdew-Wang '91 [G] Perdew-Burke-Ernzerhof [G] Hamiltonians • Restricted and Unrestricted Hartree-Fock Theory • Total and Spin Density Functional Theory • Local functionals [L] and gradient-corrected [G]exchange-correlation functionals • Hybrid DFT-HF exchange functionals • B3PW, B3LYP (using the VWN5 functional) • User-definable hybrid functionals Vallico Sotto July 2009

  26. The basis set Crystalline orbitals as linear combinations of Bloch Functions as linear combinations of Atomic Orbitals as contractions of Hermite Gaussian functions Vallico Sotto July 2009

  27. Software performance Memory management: dynamic allocation Efficient storage of integrals or Direct SCF Full parallelization (MPI) Replicated data version Massive parallel version up 2048 processors (soon available) Supported platforms Pentium and Athlon based systems with Linux IBM workstations and clusters with AIX 4.2 or 4.3 SGI workstations and servers DEC Alpha workstations HP-UX systems Sun Solaris Linux Alpha Running CRYSTAL2006 Vallico Sotto July 2009

  28. The problem of H It is well known that the stretching modes involving hydrogen atoms are strongly anharmonic: typically for the O-H stretching anharmonicity can be as large as 180 cm-1. However this difficulty is compensated by the full separability of this mode. Vallico Sotto July 2009

  29. Anharmonic correction for hydroxyls OH stretching is considered as decoupled from any other normal modes exe=(2 01- 02) / 2 E2 A wide range (0.5 Å) of OH distances must be explored to properly evaluate E1 and E2 02 E1 01 E0 Direct comparison with experiment for fundamental frequency, first overtone and anharmonicity constant This procedure is automatically implemented in the code Vallico Sotto July 2009

  30. H O M Edingtonite surface Chabazite M=Mg Brucite M=Ca Portlandite Isolated OH groups in crystals: model structures/1 All calculations with 6-31G(d,p) basis set Vallico Sotto July 2009

  31. B3LYP vs experimental OH frequencies Vallico Sotto July 2009

  32. Is the choice of the Hamiltonian critical? BRUCITE, Mg(OH)2 Fundamental OH stretching frequencies, cm-1 No hydrogen bond Vallico Sotto July 2009

  33. Is the choice of the Hamiltonian critical? Be(OH)2 Fundamental OD stretching frequencies. All data in cm-1 Hydrogen bonded OH groups • Only B3LYP is in good agreement with experimental free OH frequency • All Hamiltonians are unable to predict shifts due to strong hydrogen bond • The 1D approximation is not appropriated to describe the OH stretching properties in the case of strong interaction of the H atom (as HB). The anharmonic constant is overestimated. Vallico Sotto July 2009

  34. numerical second derivative analytical first derivative Harmonic frequencies at the central zone are obtained by diagonalising the mass weighted Hessian matrix, W Isotopic shift can be calculated at no cost! Harmonic frequency in solids with CRYSTAL Building the Hessian matrix Vallico Sotto July 2009

  35. The dynamical matrix The behavior of the phonons of a wave vector k close to the Γ point can be described as follows: Center-zone phonons: ANALYTICAL Dependence on the direction of k:limiting cases k→0 NON ANALYTICAL *greek indices: atoms in the primitive cell **latin indices: cartesian coordinates Vallico Sotto July 2009

  36. The analytical part of the dynamical matrix *Mx= mass of the x atom **H=Hessian matrix Vallico Sotto July 2009

  37. The Born charges • The atomic Born tensors are the key quantities for : • calculation of the IR intensities • calculation of the static dielectric tensor • calculation of the Longitudinal Optical (LO) modes They are defined, in the cartesian basis, as: *Ei=component of an applied external field **μ=cell dipole moment Vallico Sotto July 2009

  38. μdepends on the choice of the cell BUT the dipole moment difference between two geometries of the same periodic system (polarization per unit cell) is a defined observable. The partial second derivatives appearing in the the Born tensors are estimated numerically from the polarizations generated by small atomic displacements (the same as for the second energy derivative) LOCALIZED WANNIER FUNCTIONS (WF) to compute polarization Vallico Sotto July 2009

  39. Procedure for the polarization derivative calculation • full localization scheme for the equilibrium point → centroids of the resulting WFs • WFs of the central point are projected onto the corresponding occupied manifolds of the distorted structures → centroids of the resulting WFs • difference between the sum of the reference WF centroids at the two geometries • C.M. Zicovich-Wilson, R. Dovesi, V.R. SaundersA general method to obtain well localized Wannier functions for composite energy bands in linear combination of atomic orbital periodic calculationsJ. Chem. Phys., 115, 9708-9719 (2001) • Alternative scheme; through Berry phase • S. Dall’Olio, R. Dovesi, R. Resta • Spontaneous polarization as a Berry phase of the HF wavefunction. Phys, Rev B56, 10105 (1997) Vallico Sotto July 2009

  40. The non-analytical contribution and the LO modes Dynamical matrix: Analytical contribution: Non-analytical contribution: Vallico Sotto July 2009

  41. Trasverse Optical (TO) modes: the non-analytic part vanishes K and Zp,m are perpendicular Longitudinal Optical (TO) modes: the non-analytic part is ≠0 K and Zp,m are parallel Vallico Sotto July 2009

  42. Frequencies, symmetry analysys, IR intensities, IR and Raman activities CRYSTAL frequency calculation output Vallico Sotto July 2009

  43. AIMS • Document the numerical stability of the computational process • Document the accuracy (with respect to experiment, when experiment is accurate) • Interpret the spectrum and attribute the modes Vallico Sotto July 2009

  44. Garnets: X3Y2(SiO4)3 Space Group: Ia-3d 80 atoms in the primitive cell (240 modes) Γrid = 3A1g + 5A2g + 8Eg + 14 F1g + 14 F2g +5A1u + 5 A2u+ 10Eu + 18F1u + 16F2u 17 IR (F1u) and 25 RAMAN (A1g, Eg, F2g) active modes Vallico Sotto July 2009

  45. Silicate garnet grossular structure: Ca3Al2(SiO4)3 O Si O Ca tetrahedra O distorted dodecahedra Al • Cubic Ia-3d • 160 atoms in the UC (80 in the primitive) • O general position (48 equivalent) • Ca (24e) Al (16a) Si (24d) site positions octahedra Vallico Sotto July 2009

  46. The interest for garnets+TM compounds • M.D. Towler, N.L. Allan, N.M. Harrison, V.R. Saunders, W.C. Mackrodt and E. ApràAn ab initio Hartree-Fock study of MnO and NiO.Phys. Rev. B 50, 5041-5054 (1994) • R. Dovesi, J.M. Ricart, V.R. Saunders and R. OrlandoSuperexchange interaction in K2NiF4 . An ab initio Hartree-Fock study J. Phys. Cond. Matter 7, 7997-8007 (1995) • Ph. D'Arco, F. Freyria Fava, R. Dovesi and V. R. SaundersStructural and electronic properties of Mg3Al2Si3O12 pyrope garnet: an ab initio studyJ. Phys.: Cond. Matter 8, 8815-8828 (1996) Vallico Sotto July 2009

  47. Symmetry is crucial for solids R. Dovesi On the role of symmetry in the ab initio Hartree-Fock linear combination of atomic orbitals treatment of periodic systems. Int. J. Quantum Chem. 29, 1755-1774 (1986).INTEGRALS C. Zicovich-Wilson and R. Dovesi,On the use of symmetry adapted crystalline orbitals in SCF-LCAO periodic calculations. I. The construction of the symmetrized orbitals.Int. J. Quantum Chem. 67, 299-309 (1998). K SPACE DIAG-IRREPS C. Zicovich-Wilson and R. Dovesi,On the use of symmetry adapted crystalline orbitals in SCF-LCAO periodic calculations.II. Implementation of the Self-Consistent-Field Scheme and examplesInt. J. Quantum Chem. 67, 311-320 (1998). SYM LABELS TO STATES R. Dovesi, F. Pascale, C. M. Zicovich-Wilson The ab inizio calculation of the vibrational spectrum of cristalline compounds; the role of symmetry and related computational aspects. Beyond standard quantum chemistry: applications from gas to condensed phases ISBN: 978-81-7895-293-2 Editor: Ramon Hernandez-Lamoneda (2007) HESSIAN Vallico Sotto July 2009

  48. Hessian construction and Symmetry(Garnet example) Each SCF+Gradient calculation provides one line of Hik 80 atoms = 240+1 SCF+G calculations with low (null) symmetry • Point symmetry is used to generate lines of atoms symmetry related • Other symmetries (among x, y, z lines; translational invariance) further reduce the number of required lines At the end only 9 out of 241SCF+G calculations are required Vallico Sotto July 2009

  49. Cost of the calculations The 9 SCF+GRAD calculations: Spessartine (open shell). Elapsed time, in seconds, per point and per processor. Load balancing Total CPU time NODE 0 CPU TIME = 286054.450 NODE 1 CPU TIME = 284862.080 NODE 2 CPU TIME = 285570.430 NODE 3 CPU TIME = 285803.140 ... ... ... ... ... ... 79 h = 3.3 days on 16 processors (Dual Core AMD Opteron 875, 2210 Mhz, 64 bit, shared memory) Vallico Sotto July 2009

  50. Numerical stability of the computational process is the mean absolute deviation of frequencies between 2 values of the indicated option. (in cm-1) • DFT integration grid Large =0.7 XLarge Standard =0.6 • Standard grid is enough(For the pyrope case) • SCF convergence (total energy, in hartree) =0.2 Tol∆E=10-10 Tol∆E=10-11 Vallico Sotto July 2009

More Related