Modelling polarization effects in molecular simulations with ab initio parametrized approaches

1. CECAM Workshop Ab Initio Meets Classical Simulations: The Development of Empirical Potentials for Atomistic Systems Modelling polarization effects in molecular simulations with ab initio parametrized approaches Gloria Tabacchi Lyon Oct.17-19 2005

2. First principlesvs. Empirical Potentials  • low computational cost: • long time scale simulation • of large systems  • Electronic structure is directly accessible • Allow to accurately describe chemically complex systems Polarizable models • Pair additivity approximation • Limited accuracy • Transferability problems • Parametrization  • High computational cost: small systems (~102 atoms) 

3. Polarizable models: many have been proposed! A promising starting point could be Use a DFT- based framework • Provides a general scheme from which different models can be derived • CPE methodsD.M. York, W. Yang, JCP 104, 159 (1996) • 2nd order expansion of E[ρ] in dr and dn • a polarizable force field model • a polarizable Kim-Gordon method

4. The system is defined as a collection of N individual subsystems (atoms, ions or molecules) Electronic density: sum of a reference density r0 (unperturbed, isolated subsystems) and a response density dr(polarization) Additivity of subsystems’ rA0 Additivity of subsystems’ drA drA r0A Definition of the system

5. Derivation of the energy expression • Start from the energy functional of the system with: 2. expand up to 2nd order in dr and dn, using FSR(1)[r] FSR(2)[r] 3. Perform approximations and solve for dr

6. A DFT-based polarizable force field modelG.Tabacchi, J.Hutter, CJ Mundy, M. Parrinello JCP 117, 1416 (2002) • Expand the response density in a finite basis of localized functions ⊲Derive from this expression a computationally efficient scheme for molecular simulations satisfying the linear response property • Define the external potential n and the total density rt

7. collect electrostatic terms in EH and get: • introduce effective densities: • Choose approximated expressions for the short range terms FSR(n) Approximations SR terms Extended Hückel approximation (F: empirical parameters) Nab is related to the hardness matrix: t is a property of the subsystem and does not depend on atomic positions. Modelled by a pair potential Φ(R) (to be parametrized) FSR[r0]=T+EXC

8. Energy expression for classical DFT based polarizable model: • c’s: variables describing polarization. • For a given ionic configuration, E=Min E[c]→ MD implementation →extended Lagrangian • c treated as additional dynamical variables • m: associated fictitious mass • get equations of motion for {R}, {c} & propagate • Apply Nose-Hoover thermostat to c’s to prevent • energy transfer between ionic/fictitious DOF

9. To be defined (parametrized) • Response basis functions φi • Extended Hückel Parameters Fa • First order SR kernel t • Effective charge distributions • Empirical two-body force field Φ DFT/LR-DFT calculations on the isolated subsystem Force matching to a database of first-principles forces from FP-MD simulations on the full system • For each configuration:Minimize E(c)Calculate the contributions to the forces depending on c,FcReplace FREF withFREF ’ =FREF-Fc and minimize Ω w.r.t. p

10. Application: LiI A challenging test case: • Standard chemistry textbook: Li-I is considered a ionic compound. • However, LiI is not properly described by a “Rigid Ion” model • I-: the largest and most polarizable among the halides • Li+: largest charge/mass ratio among the alkali cations, i.e. the highest polarizing ability Phases/systems studied • LiI molecule in the gas phase • Crystal (B1) and liquid phase Properties: • Energy vs r. curve of LiI molecule • Structural and thermodynamic properties of condensed phases Choice of reference subsystems: Li+ and I-

11. Parametrization of LiI ⊲ LDA approximation ⊲ PW (30 Ry cutoff) ⊲ Norm conserving pseudopotentials with p nonlocality ⊲ Dt=0.12 fs; m=500 au Reference KS-DFT calculations performed on: • the isolated subsystemsI-, Li+(LR) • LiI gas phase molecule • Condensed phases (32 LiI per unit cell) liquid (T=800 K and 3000 K) solid (rocksalt), 300 K CP-MD simulations (20 ps) 1. Get the parameters of the subsystems I- response basis {fI} → 1 s-type CCGF (npgf=5) + px, py, pz (npgf=2) Li+ response basis {fLi} → none or 1 s-type CCGF (npgf=5) → 1 pgf integrating to -1 and +1 respectively tI-, tLi+ → TF approximation for T FI-, FLi+ →

12. 2. Get the parameters of the pairwise term The database of first principles forces contains: • 30 configs. for the LiI molecule (different rLi-I) • 8 configs. taken from FP-MD simulations of the crystal • 34 configs taken from FP-MD simulations of liquid LiI Φ parameters from force-matching 3. Test the full set of parameters 4. Run simulations and compare to KS-DFT results 5. Can PFF do better than a Rigid Ion model? A non polarizable two body force field was derived by force matching to the same force database

13. LiI molecule: results • First set of parameters (PFF0): OK for molecule, NOT for condensed phase • solution: • keep the same response basis, readjust F,t, Φ by force matching • New set PFF1: transferable to condensed phases.

14. r(Å) r(Å) PFF-MD simulations results • conservative dynamics (~ 10-5 au) (Dt=0.6 fs; m=10-100 au; NH thermostats on c) • PFF can satisfactorily reproduce structural properties predicted by KS-DFT • transferability of PFF1 among different LiI phases LiI (solid) LiI (liquid) PFF can well reproduce the different shape of the like-ions g(r) in the solid phase Short-distance tail in the FPMD Li-Li g(r) due to formation of multicenter bonds: CT interactions are by construction out of the scopes of PFF

15. Calculation of thermodynamic properties and comparison with experiment • Aim: build a model able to reproduce results of first-principles MD simulations on condensed phase systems • Model: validated through comparison to FPMD simulations • How such model can reproduce experimental results? ⊲Calculation of the relative stability of B1, B2, B4 structures the model predicts the correct crystal structure! ⊲ Calculation of representative thermodynamic properties: Rather satisfactory agreement NB: no experimental data was used in the parametrization!

16. DFT-derived scheme Basis functions for dr consistent with LR theory Conservative, fast MD A single parametrization describes polarization effects in different phases of LiI Satisfactory comparison with KS-DFT/ experiment Use of a pair potential for non electrostatic energy contributions due to r0 Force matching needs an extended database of FP forces (full system) Parametrization is not trivial and application to much complex systems not straightforward DFT-based PFF: +’s and –’s  

17. How to improve the model Wishlist: • More accuracy and transferability • Simpler parametrization • avoid the pairwise additivity approximation (no force matching to FP forces of full system) • Parameters of subsystems only Explicit treatment of SR (T+EXC) terms r0: frozen electron gas approach (Kim-Gordon)

18. Kim-Gordon method (KG)R.G. Gordon, Y.S. Kim, JCP 56, 3122 (1972) Y.S. Kim, R.G. Gordon JCP 60, 1842 (1974) • additivity of frozen electronic densities of individual atoms • local density approximations for T, EXC; rA from HF calculations • works well for noble gas atoms Empirical Kim-Gordon method (EKG)D.Barker, M. Sprik Mol. Phys. 101, 1183 (2003) • Treats subsystem density (atomic or molecular) • as an adjustable parameter. • Application: liquid water. Reference density was • adjusted to reproduce the dipole moment of H2O • in liquid water (change dipole from 1.8 to 2.9 D)

19. drA r0A Polarizable Kim-Gordon method (PKG)G.Tabacchi, J.Hutter, C.J. Mundy JCP 123, 074108 (2005) Idea: Response density: defined as in PFF Frozen reference density “a-la KG” • r0 models the valence electron density: not an adjustable parameter • Core electrons contribution: → define a local pseudopotential n=nc+npsSR • nc due to the core charge rc • Define rt and collect electrostatic terms in Ees • start from

20. PKG requires: • Basis set for response density • definition of subsystems’ reference density • Local pseudopotentials Derived from calculations on the individual subsystems: no KS-DFT on the specific system under investigation Subsystems’ r0 • Calculate the KS-DFT valence density of each subsystem • Parabolic confining potential for negatively charged subsystems • Project onto a finite basis of cartesian Gaussians (ng=5-6)

21. Local atomic pseudopotentials • next should also account for core electrons • Atomic PP of conventional KS-DFT are non local • Introduce a novel type of local atomic PP nc npSR • Use a generalized form of the local part of GTH PP • S. Goedecker, M. Teter, J. Hutter PRB 54, 1703 (1996) • Follow the GTH prescription: • the potential parameters should minimize the difference between all-electron and pseudo eigenvalues/charges for all valence states

22. Test PKG model: alkali halidesLiI & LiCl molecules • Reference KS-DFT calculations: LDA, PW (60 Ry), GTH PP (Li semicore) • PKG calculations: TF functional for KE. Reference subsystems: M+, X- • same Li parameters (PP, r0, response basis) adopted for LiCl, LiI • Satisfactory agreement with KS-DFT (with TF!) • Local PPs (obtained for neutral atoms) transferable to ions • Accuracy ≳ PFF, but no force matching!

23. PolKG-MD simulations: LiI 32 formula units as in reference calculations. Dt=0.6 fs, m=100 au. LiI (solid) LiI (liquid) -- PolKG ― KS-DFT -- PolKG ― KS-DFT PKG-MD: good agreement with FPMD structural properties

24. PolKG-MD simulations: LiCl LiCl (solid) LiCl (liquid) -- PolKG ― KS-DFT -- PolKG ― KS-DFT PKG-MD: good transferability of subsystems’ parameters

25. Thermodynamic properties and comparison with experiment • Correct crystal structure is predicted • satisfactory agreement with KS-DFT & EXP

26. PKG: an improvement over PFF • no pair potentials • Better accuracy • Much simpler parametrization (only on isolated subsystems) • Transferability of subsystems’ parameters among different compounds • still competitive in terms of computational convenience

27. CONCLUSION • computationally cheap models derived from approximations to DFT can reproduce properties of polarizable ionic systems with an accuracy comparable to KS-DFT in perspective: • Modelling more complex systems: should be feasible even though non-trivial • Interface to KS-DFT

28. Acknowledgments • J. Hutter, C.J. Mundy, M. Parrinello • CP2K developers’ group implementation available through the open source project CP2K: http://cp2k.berlios.de/ • CECAM, COST-MolSimu Thanks for your attention

29. Polarization effects in large scale systems • Speed up first principles (order-N) • Build better empirical force fields • Develop novel techniques aimed at combining the advantages of both approaches Goal: achieve an explicit (and accurate) description of the distortions in the electronic valence density of a system due to changes in its chemical environment “approximated DFT methods”

30. Force matchingF.Ercolessi, J.B. Adams Europhys. Lett. 26, 583 (1994) Given an empirical potential characterized by a set of parameters {p}, determine the optimal set {p*} that match the forces FREF generated from reference first principles calculations on a set of ionic configurations, i.e: Nc:number of configurations Nj: atoms in configuration j In our scheme polarization effects are explicitly accounted for: Parameters of the pair potential F have a different physical meaning w.r.t. a traditional rigid ion model Fc • For each configuration in the database:Minimize E(c)Calculate the contributions to the forces depending on c,FcReplace FREF withFREF ’ =FREF-Fc and minimize Ω w.r.t. p

32. PKG method: energy and forces on coefficents

33. PKG method: Forces on particles