Ab Initio Molecular Dynamics with a Continuum Solvation Model

1. CSC 2002 Conference, Vancouver, BC, 2 June 2002 Ab Initio Molecular Dynamics with a Continuum Solvation Model Hans Martin Senn, Tom Ziegler University of Calgary

2. Outline • Background • Car–Parrinello ab initio molecular dynamics • The Projector-Augmented Wave (PAW) method • Continuum solvation / The Conductor-like Screening Model (COSMO) • Continuum solvation within ab initio MD • Surface with smooth analytic derivatives • Surface charges as dynamic variables • Tests and first applications • Hydration energies • Conformers of glycine • CO insertion into Ir–CH3

3. Nuclei Classical mechanics Electrons Density-functional theory Ab Initio Molecular Dynamics R. Car, M. Parrinello, Phys. Rev. Lett. 1985, 55, 2471 • Ab initio molecular dynamics Ab initio molecular dynamics Forces on nuclei derived frominstantaneous electronic structure • Car–Parrinello scheme for ab initio MD • Fictitious dynamics of the wavefunctions • Simultaneous treatment of atomic and electronic structure • Friction dynamics for minimizations

4. The Projector-Augmented Wave (PAW) Method P.E. Blöchl, Phys. Rev. B 1994, 50, 17953 • Basis set: Augmented plane waves • Spatial decomposition of the wavefunction according to its shape • Smooth plane-wave expansion far from nuclei • Rapidly oscillating augmentation functions around the nuclei • Integrations for total energy and matrix elements decompose accordingly • Smooth parts in Fourier representation • One-centre parts on radial grids in spherical harmonics expansions • All-electron method • Frozen-core approximation

5. Continuum Solvation Models • Implicit electrostatic solvation • “Solvent” = homogeneous, isotropic dielectric medium • Fully characterized by relative permittivity (r) • Electrostatic interactions • Solute polarizes the dielectric ( reaction field) Epol > 0 • Solute charge density interacts with polarized medium Eint < 0 • Solute electronic and atomic structure adapt (“back-polarization”) ∆Een > 0 • Repeat to self-consistency! • Electrostatic free energy • For a linear dielectric:

6. The Solvation Energy • Non-electrostatic contributions • Cavity formation • Dispersion and exchange-repulsion • Approximately proportional to cavity surface area • Free energy of solvation • Neglect changes in thermal motion upon solvation

7. The Conductor-Like Screening Model (COSMO) • Main approximations in COSMO 1. Dielectric with infinite permittivity (er ∞) • No volume polarization; polarization expressed as surface charge density only • Vanishing potential at the cavity boundary 2. Recover true dielectric behaviour by scaling • Derived for rigid multipoles in spherical cavity • f = f(er) = (er – 1) / (er + x) [x = 0, 1/2 for monopole, dipole] 3. Discretize cavity surface and surface charge • Surface segments carrying point charges qi • Energy expression is variational in qi

8. Surface Charges as Dynamic Variables • Extended Lagrangean including solvation • Equations of motion for charges • Additional forces acting on the nuclei • Surface segments and charges implicitly depend on atomic positions via the surface construction • Forces obtained as analytic derivatives • Energy-conserving dynamics requires smooth derivatives • Segments/charges must not vanish or be created during the simulation • Each segment/charge must be unambiguously assigned to an atom • Conventional schemes for building the cavity are not compatible with this!

9. Surface Construction • A surface construction compatible with AIMD • Surface = union of atomic spheres • Each sphere is triangulated • All segments and charges are retainedEach segment/charge moves rigidly along with its atom 4. A switching function effectively removes charges lying within another sphere from the energy expression • Qi = 1 if the charge is exposed on the molecular surface • Qi = 0 if the charge is buried • Smooth transition between “on” and “off”

10. Switching Functions • Building the switching function • Rectangular pulse • Centred at d0 = Rsolv – c • Vanishes smoothly at d0 n = 1 n = 10 n = 24 h d • Atomic switching function d0 Rsolv • Total switching function • Charge is smoothly switched off if it lies within any one of the other spheres RAsolv RBsolv

11. Rijc pot. |si – sj| Rounding the Edges… • Behaviour of hidden charges • Magnitude of switched-off charges physically undetermined • Problematic if charge becomes again exposed • Restoring potential • Added to the Lagrangean • Artificial self-interaction, acting only on hidden charges • Modified Coulomb potential at short range • Charges can collide at seams • Modified Coulomb potential • True 1/r behaviour replaced by linear continuationfor |si – sj| < Rijc • Rijc derived from average of self-interaction potentials

12. Representation of the Solute Density • Decoupling of periodic images • Plane-wave-based methods always create periodic images • Artificial in molecular calculations • Electrostatic decoupling scheme in PAW • Model density • Superposition of atom-centred Gaussians • Coefficients determined by fit to true density • Reproduces long-range behaviour • Surface charges couple to model density • Not instrumental • Computationally convenient • 3–5 Gaussians per atom • Fitting procedure relays forces acting on the wave functions due to the charges

13. 4.2 kJ/mol MUE: 3.2 kJ/mol 5.2 kJ/mol Parameterization and Validation • Parameterization • Free energies of hydration • Parameters: Solvation radii RIsolv, non-polar parameters b,g • Fit set: Neutral organic molecules containing C, H, O, N • Mean unsigned error (MUE) over whole set: 4.1 kJ/mol • On par with similar DFT/COSMO results

14. A C Z B Conformers of Glycine • H-bonding patterns • Relative stabilities (kJ/mol) • Hydration energies (kJ/mol) • Zwitterion in solution • Most stable conformer • Not present in the gas phase

15. Energy Conservation • Conservation of total energy • Glycine without and with continuum solvation • ∆t = 7 a.u., T ≈ 300 K (no thermostats, no friction) • Drift in total energy: solvation off +2.78  5Eh/ps per atom solvation on 1.43  5Eh/ps per atom • Energy is conserved

16. RC A First Application • Carbonylation of methanol • “Monsanto” / ”Cativa” acetic-acid process • M = Rh, Ir • Via [MI3(CH3)(CO)2]– [MI3(COCH3)(CO)]– • AIMD simulation of CO insertion step • Solvent: iodomethane; T = 300 K • Slow-growth thermodynamic integration along d(CMe–CCO) • ∆A‡ = 111  (128) kJ/mol∆E‡ = 126  (155) kJ/mol; ∆S‡ = 60  (91) J/(K mol) • Dissociation of I– trans to incipient acyl! • Only in solution at finite temperature • Enthalpy of Ir–I bond traded for entropy of free (solvated) I–

17. Summary • COSMO continuum solvation in AIMD • Surface charges as dynamic variables • Surface construction • Switching function smoothly disables unexposed charges • Smooth analytic derivatives wrt. atomic positions • Energy-conserving • Solvation energies • On par with other DFT/COSMO implementations • Modelling of finite-temperature and solvation effects

18. Acknowledgments • Peter E. Blöchl,Institute of Theoretical Physics, Clausthal University of Technology, Germany • Peter M. Margl,Corporate R&D Computing, Modeling and Information Sciences, Dow • Rochus Schmid,Institute of Inorganic Chemistry, Munich University of Technology, Germany • Swiss National Science Foundation • NSERC • Petroleum Research Fund