Hybrid QM/MM Methods

1. Hybrid QM/MM Methods Paul Sherwood CLRC Daresbury Laboratory p.sherwood@dl.ac.uk

2. Overview • Overview of QM/MM Methodology • Historical Overview • QM/MM Coupling and link atoms • Choice of QM and MM models • Conformational complexity • Implementational issues • Experience with different approaches • GAMESS-UK/CHARMM • ChemShell • Exploitation of HPC • The QUASI project • Objectives, partners, timelines, funding, and workshop • Some recent QM/MM applications • Concluding remarks

3. The QM/MM Modelling Approach • Couple quantum mechanics and molecular mechanics approaches • QM treatment of the active site • reacting centre • excited state processes (e.g. spectroscopy) • problem structures (e.g. complex transition metal centre) • Classical MM treatment of environment • enzyme structure • zeolite framework • explicit solvent molecules • bulky organometallic ligands

4. Historical Overview • 1976 Warshel, Levitt MM+Semi-empirical • study of Lysozyme • 1986 Kollmann, Singh QUEST • ab-initio (Gaussian-80/UCSF) + AMBER • 1990 Field, Bash, Karplus CHARMM/AM1 • full semi-empirical dynamics implementation • 1992 Bernardi, Olivucci, Robb MMVB • simulation on MC-SCF results using VB • 1995 van Duijnen, de Vries HONDO/DRF • direct reaction field model of polarisation • 1995 Morokuma IMOMM • mechanical embedding with “hidden variable” optimisation • 1996 Bakowies and Thiel MNDO/MM • mechanical, electrostatic and polarised semi-empirical models • 1997 Eichler, Kölmel, and Sauer QM/Pot • subtractive coupling of GULP/TURBOMOLE

5. QM/MM - One idea, Many schemes! Termination Atoms Chemical type hydrogen atoms,… MM Charge perturbations charge deletion, charge shift, double link atoms… Constraint • QM/MM Couplings • Unpolarised • QM polarisation • choice of charges? • MM polarisation • shell model • dipole polarisabilities M2 Q2 M1 L Q1 M2 Q2 Total Energy Expression “Additive” E(M,MM) + E(QL,QM) Link atom corrected E(M,MM) + E(QL,QM) - E(L,MM) “Subtractive” E(MQ,MM) + E(QL,QM) - E(QL,MM)

6. Termination of the QM region • Boundary region approaches • Boundary atoms are present in both QM and MM calculations • Range of representations within QM code • Modified ab-initio atom with model potential • Semi-empirical parameterisation • Frozen orbitals • Re-parameterised MM potentials • Link atom schemes • Terminating (link) atoms are invisible to MM calculations • Hydrogen, pseudo-halogen, etc.

7. Energy Expression E(OI,MM) + E(IL,QM) - E(IL,MM) includes link atom correction can treat polarisation of both the MM and QM regions at the force-field level Termination Any (provided a force field model forIL is available) Advantages Potentially highly accurate and free from artefacts Can also be used for QM/QM schemes (e.g. IMOMO, Morokuma et al) Disadvantages Need for accurate forcefields (mismatch of QM and MM models can generate catastrophes on potential energy surface) No electrostatic influence on QM wavefunction included Applications Zeolites (Sauer et al) Organometallics (Morokuma et al) Subtractive QM/MM Coupling

8. Additive Energy Expressions • Energy Expressions • Without link atom correction E(O,MM) + E(I,QM) + E(IO,QM/MM) • Link atom correction E(O,MM) + E(I,QM) + E(IO,QM/MM) - E(L,MM) • Boundary methods E(OB,MM) + E(IB,QM) + E(IBO,QM/MM) • Highly variable in implementation • QM/MM couplings, • QM termination etc • Advantages • No requirement for forcefield for reacting centre • Can naturally build in electrostatic polarisation of QM region - effects of environment of excitations etc • Disadvantages • Electrostatic coupling of the two regions, E(IO,QM/MM) is problematic with link atoms • Need for boundary atom parameterisation • Applications • Solvation • Enzymes • Zeolites

9. QM/MM Non-bonded Interactions • Short-range forces (van der Waals) • Typically will follow MM conventions (pair potentials etc), sometimes reparameterisation is performed to reflect replacement of point charges interactions with QM/MM electrostatic terms. • Electrostatic interactions: • Mechanical Embedding in vacuo QM calculation coupled classically to MM via point charges at QM nuclear sites • Electrostatic Embedding MM atoms appear as centres generating electrostatic contribution to QM Hamiltonian • Polarised Embedding MM polarisability is coupled to QM charge density

10. Mechanical Embedding Advantages • MM and QM energies are separable • separate MM relaxation, annealing etc possible • QM/MM terms can be integrated directly into the forcefield • No interactions between link atoms and MM centres • QM energies, gradient, Hessian are the same cost as gas phase Drawbacks • No model for polarisation of QM region • QM/MM electrostatic coupling requires atomic charges for QM atoms • generally these will be dependent on reaction coordinate Examples • IMOMM (Morokuma) • MNDO/MM (Bakowies and Thiel)

11. Electrostatic Embedding (i) Assign MM Charges for pure MM system • Derived from empirical schemes (e.g. as part of forcefield) • Fitted to electrostatic potentials • Formal charges (e.g. shell model potentials) • Electronegativity equalisation (e.g. QEq) (ii) Delete MM charges on atoms in inner region • Attempt to ensure that MM “defect” + terminated QM region has • correct total charge • approximately correct dipole moment (iii) Insert charges on MM centres into QM Hamiltonian • Explicit point charges • Smeared point charges • Semi-empirical core interaction terms • Make adjustments to closest charges (deletion, shift etc)

12. Creation of neutral embedding site (i) Neutral charge groups • Deletion according to force-field neutral charge-group definitions R O H C C N N C C H R O

13. Creation of neutral embedding site (i) Neutral charge groups • Total charge conserved, poor dipole moments R O H C C N H H N C C H R O

14. Creation of neutral embedding site (ii) Polar forcefields bond dipole models, e.g. for zeolites (Si +0.5x, O -0.5x) O-x O-x Si+2x O-x O-x

15. Creation of neutral embedding site (ii) Polar forcefields O-0.5x H O-0.5x H Si H O-0.5x H O-0.5x

16. Creation of neutral embedding site (iii) Double link atoms • Suggestion from Brooks (NIH) for general deletion (not on a force-field neutral charge-group boundary) R O H C C N N C C H R O

17. Creation of neutral embedding site (iii) Double link atoms • All fragments are common chemical entities, automatic charge assignment is possible. O R H C H C H C N H H N C R H O

18. H H H H C C H H Creation of neutral embedding siteDouble Link Atoms • Conventional QM/MM schemes • Break into H3Cand CH3 • Neutral fragments (in this simple case) • Non-zero individual dipoles • Replace QM CH3 with some form of terminated group (L-CH3) • Finite total dipole moment • Often further adjustment to MM charges is required (may create additional charge and dipole errors). • Double link atom • Dipole generally well approximated by superposition of bonds H H H H H C C H H H H H H H H C C H H QM MM

19. Boundary adjustments • Some of the classical centres will lie close to link atom (L) • Artefacts can result if charge at the M1 centre is included in Hamiltonian, many adjustment schemes have been suggested • Adjustments to polarising field can be made independently from specification of MM…MM interactions • Similar adjustments may are needed if M1 is classified as a boundary atom, depending on M1 treatment. Q2 M2 M1 L Q1 M2 Q2 M3 Q3

20. Boundary Adjustments (i)Selective deletion of 1e integrals • L1: Delete integrals for which basis functions i or j are sited on the link atom L • found to be effective for semi-empirical wavefunctions • difference in potential acting on nearby basis functions causes unphysical polarisation for ab-initio QM models • L3: Delete integrals for which basis functions i and j are cited on the link atom and qAis the neighbouring MM atom (M1) • less consistent results observed in practice † † Classification from Antes and Thiel, in Combined Quantum Mechanical and Molecular Mechanical Methods, J. Gao and M. Thompson, eds. ACS Symp. Ser., Washington DC, 1998.

21. Boundary Adjustments (ii)Deletion of first neutral charge group • L2 - Exclude charges on all atoms in the neutral group containing M1 • Maintains correct MM charge • leading error is the missing dipole moment of the first charge group • Generally reliable • free from artefacts arising from close contacts • Limitations • only applicable in neutral group case (e.g. AMBER, CHARMM) • neutral groups are highly forcefield dependent • problematic if a charge group needs to be split • Application • biomolecular systems

22. Boundary adjustments (iii)Charge shift • Delete charge on M1 • Add an equal fraction of q(M1) to all atoms M2 • Add correcting dipole to M2 sites (implemented as a pair of charges) • charge and dipole of classical system preserved • Leading sources of residual error is that Q---L dipole moment is not equivalent to Q------M Q2 M2 M1 L Q1 M2 Q2 M3 Q3

23. Adjusted Connection Atoms

24. Boundary adjustments (iv)Gaussian Blur • Delocalise point charge using Gaussian shape function • Large Gaussian width : electrostatic coupling disappears • Narrow Gaussian width : recover point charge behaviour • Intermediate values • short range interactions are attenuated • long range electrostatics are preserved • Importance of balance - apply to entire MM system or to first neutral group • Particularly valuable for double-link atom scheme where MM link atom charge lies within QM molecular envelope

25. Gaussian Blur: QM/MM Ethane H H Aggregate dipole moment, with MM atoms broadened • Large Gaussian width dipole converges to that of the QM methane ( 0 ) • Small Gaussian width point charges polarise the QM region, away from the C-H bond H H H H C C H H MM QM Compare with dipole of MM methyl group (0.18)

26. Localised Orbital Approaches

27. Electrostatic EmbeddingSummary Advantages • Capable of treating changes in charge density of QM • important for solvation energies etc • No need for a charge model of QM region • can readily model reactions that involve charge separation Drawbacks • Charges must provide a reliable model of electrostatics • reparameterisation may be needed for some forcefields • Danger of spurious interactions between link atoms and charges • QM evaluation needed to obtain accurate MM forces • QM energy, gradient, Hessian are more costly than gas phase QM

28. Polarised embedding schemes • Incorporate polarisation of classical region • most appropriate used when the forcefield itself is based on explicit polarisability • back-coupling of polarised charge density to QM calculation can sometimes be omitted • Approaches • Iterative solution of dipole polarisabilities • Direct Reaction Field Hamiltonian (van Duijnen, de Vries) • solution of coupled polarisabilities using relay matrix • possibility of including 2-electron dispersion terms • implemented in HONDO and GAMESS-UK • Shell model-based schemes • atomic charge is split into core and valence electron shell, connected by a harmonic spring • QUASI solid-state embedding scheme

29. Classical cluster termination Base model on finite MM cluster QM region sees fitted correction charges at outer boundary QM region termination Ionic pseudopotentials (e.g. Zn2+, O2-) associated with atoms in the boundary region Forcefield Shell model polarisation Classical estimate of long-range dielectric effects (Mott/Littleton) Energy Expression Uncorrected Advantages suitable for ionic materials Disadvantages require specialised pseudopotentials Applications metal oxide surfaces Solid-state Embedding Scheme MM QM

30. Implementation of solid-state embedding • Under development by Royal Institution and Daresbury • Based on shell model code GULP, from Julian Gale (Imperial College) • Both shell and core positions appear as point charges in QM code (GAMESS-UK) • Self-consistent coupling of shell relaxation • Import electrostatic forces on shells from GAMESS-UK • relax shell positions GAMESS-UK SCF & shell forces GULPshell relaxation GAMESS-UK atomic forces GULP forces

31. Polarised Embedding SchemesSummary • Advantages • more accurate treatment of solvation effects • allows coupling to systems where the best forcefields are based on polarisation (e.g. shell model potentials for metal oxide systems) • Drawbacks • Additional cost • solution of coupled polarisabilities • some schemes will require additional SCF iterations • requirement for polarised force-field • danger of electrostatic instabilities close to boundaries • difficult to apply reliably when using link atoms

32. QM Termination Schemes • Link atom schemes • Hydrogen atoms • Adjusted electronegativity • Hamiltonian shift operator • pseudohalogen (Hyperchem) • Methyl groups (Cummins, Gready) • Boundary schemes • Frozen Orbitals • Local SCF scheme (Rivail) • Generalised hybrid orbital (Gao) • ab-initio implementation (Friesner) • Pseudopotentials • Gaussian basis (Yang), Plane-wave (Rothlisberger) • Adjusted connection atoms (Thiel) • semi-empirical mimic for attached methyl group

33. Positioning of link atoms • Initial placement • Usually on terminated bond • Unconstrained • Additional degrees of freedom present in geometry optimisation and MD • e.g. CHARMM, QUEST • Constrained • Need to take into account forces on link atoms, • shared internal coordinate definitions (IMOMM) • chain-rule differentiation (QM/Pot, ChemShell)

34. Deletion of MM terms • General rule: Keep all MM terms involving one or more M atoms • When using link atoms constrained to bond direction, delete (or adjust) to avoid double counting • angle Q2-Q1-M1 (duplicates Q2-Q1-L) • torsion Q3-Q2-Q1-M1 • When L is free to optimise sometimes additional restraint terms (e.g. Q2-Q1-L) are added to ensure stability. Q2 M2 M1 L Q1 M2 Q2 M3 Q3

35. Choice of QM Model Applicability • Most QM methods are suitable • semi-empirical • Empirical valence bond (Warshel, MOLARIS) • MM-VB (Robb, fitted to CASSCF) • ab-initio • DFT • Gaussian basis • Plane wave (CP) - Zeigler, Parrinello, Rothlisburger Attributes • For electrostatic embedding need to insert extra nuclei in Hamiltonian • Cost implications, eg derivatives, Hessians

36. Choice of MM model Choice of parameterisation based on chemical applications • Current implementations based on • Macromolecular force fields • Charmm, Quest (AMBER) • MM2/MM3 • MNDO/MM, IMOMM • Commercial generation eg CFF (Discover) • Electronegativity equalisation • Practical considerations • for link atoms schemes must be able to remove selected forcefield terms from topology • need vdW parameters for interaction with QM (always) • QM charges and/or forcefield terms (sometimes) • numerical noise (e.g. cutoffs) important for transition states etc. • Future prospects • DMA, polarisabilities

37. QM/MM coupling - Proton Affinity Tests • Simplified alcohol test set based on [1] • AMBER charge model • QM region is small (HOCH3 in all cases) • Include systems with 1, 2 and 3 link atoms, Ethanol, i-Propanol, t-Butanol • Fixed geometries (from 3-21G QM optimisations) • Compare • Pure QM • L2 (delete first charge group) • Shift (move charge from first atom to neighbours, add dipole. • Double link atom + Gaussian Blur CH3 H H C O H CH3 H CH3 C O H CH3 H CH3 C O CH3 • [1] I. Antes and W. Thiel, in “Hybrid Quantum Mechanical and Molecular Mechanical Methods” J. Gao (ed.) ACS Symp.Ser. 712, ACS, Washington, DC, 1998.

38. Computed Proton Affinities

39. Conformational Complexity The Problem! • QM/MM schemes combine • the cost of QM • the computational complexity of large MM systems • Special treatments are needed for flexible molecules • For mechanical embedding, hidden variables can reduce number of coordinates by exploring QM PES in presence of adiabatically relaxed MM environment (eg IMOMM) • For electrostatic embedding, a charge fit to the QM charge density can be used to accelerate MM relaxation around a frozen QM core • A.J. Turner et al, PCCP 1, p1323, 1999 • Zhang et al, J. Chem. Phys., 112, p3483, 2000 • Molecular dynamics and free energy techniques will be increasingly important • computationally demanding for ab-initio QM, requirement for HPC implementations

40. QM/MM Geometry Optimisation Small Molecules • Internal coordinates (delocalised, redundant etc) • Full Hessian • O(N3) cost per step • BFGS, P-RFO • Macromolecules • Cartesian coordinates • Partial Hessian (e.g. diagonal) • O(N) cost per step • Conjugate gradient • L-BFGS • Coupled QM/MM schemes • Combine cartesian and internal coordinates • Reduce cost of manipulating B, G matrices • Define subspace (core region) and relax environment at each step • reduce size of Hessian • exploit greater stability of minimisation vs. TS algorithms • Use approximate scheme for environmental relaxation

41. HDLC Optimiser

42. Hybrid Computational Schemes Termination Atoms Chemical type hydrogen atoms, pseudopotentials adjusted connection atoms, localised orbitals Charge perturbations none, charge deletion, charge shift, selection of 1e integrals, double link atoms • QM/MM Couplings • Unpolarised • QM polarisation • choice of charges? • MM polarisation • shell model • dipole polarisabilities M2 Q2 M1 L Q1 M2 Q2 Total Energy Expression Uncorrected E(M,MM) + E(QL,QM) Boundary corrected E(M,MM) + E(QL,QM) - E(L,MM) Subtractive E(MQ,MM) + E(QL,QM) - E(QL,MM)

43. Termination of the QM region • Boundary region approaches • Boundary atoms are present in both QM and MM calculations • Range of representations within QM code • Modified ab-initio atom with model potential • Semi-empirical parameterisation • Frozen orbitals • Re-parameterised MM potentials • Link atom schemes • Terminating (link) atoms are invisible to MM calculations • Hydrogen, pseudo-halogen, etc.

44. Electrostatic EmbeddingSummary Advantages • Capable of treating changes in charge density of QM • important for solvation energies etc • No need for a charge model of QM region • can readily model reactions that involve charge separation Drawbacks • Charges must provide a reliable model of electrostatics • reparameterisation may be needed for some forcefields • Danger of spurious interactions between link atoms and charges • QM evaluation needed to obtain accurate MM forces • QM energy, gradient, Hessian are more costly than gas phase QM

45. Classical cluster termination Base model on finite MM cluster QM region sees fitted correction charges at outer boundary QM region termination Ionic pseudopotentials (e.g. Zn2+, O2-) associated with atoms in the boundary region Forcefield Shell model polarisation Classical estimate of long-range dielectric effects (Mott/Littleton) Energy Expression Uncorrected Advantages suitable for ionic materials Disadvantages require specialised pseudopotentials Applications metal oxide surfaces Solid-state Embedding Scheme MM QM

46. GULP - GAMESS-UK coupling • Features • From Julian Gale (Imperial College), classical energies including shell model • Self-consistent coupling of shell relaxation to GAMESS-UK • Import electrostatic forces on shells from GAMESS-UK • relax shell positions • Defect energies (Mott Littleton scheme) a posteriori • MPI parallel version available. • Use coupling=shift (without any bonds, and with some shells in MM region, include_qm_force=yes for GULP ChemShell - hybrid module GAMESS-UK SCF & shell forces GULPshell relaxation GAMESS-UK atomic forces GULP forces