QM/MM Modelling

1. QM/MM Modelling Paul SherwoodCLRC Daresbury Lab.P.sherwood@dl.ac.uk

2. Plan for the course • Tuesday p.m • Lecture 1 • Lecture 2 • Wednesday a.m. • Lecture 3 • Wednesday p.m. • Practical Session

3. QM/MM ModellingLecture 1Concepts and Theory

4. Overview • Rationale • Notation • inner, outer boundary, link etc • Historical overview • Energy expression • Additive vs Subtractive • QM/MM Coupling • Mechanical, electrostatic, polarised • Link Atom treatments • Choice of QM and MM models • Conformational complexity • Geometry optimisation algorithms

5. 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

6. Terminology • Usually we can view QM/MM calculations as a combination of • a QM calculation on a small region of the system • A surrounding environment modelled by a MM calculation • Optionally, some treatment, e.g. dielectric response of the region outside the outer region is included.

7. 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.

8. QM/MM - One idea, Many schemes! Termination Atoms Chemical type hydrogen atoms,… MM Charge perturbations charge deletion, charge shift, double link atoms… Constraints • QM/MM Couplings • Unpolarised • QM polarisation • choice of charges? • MM polarisation • shell model • dipole polarisabilities O2 I2 O1 L I1 O2 I2 Total Energy Expression “Additive” E(O,MM) + E(IL,QM) Link atom corrected E(O,MM) + E(IL,QM) - E(L,MM) “Subtractive” E(OI,MM) + E(IL,QM) - E(IL,MM)

9. 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

10. 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(IL,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

11. 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

12. 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

13. 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 • Shell Model

14. Valence Force Fields (i) • Used for covalently bound molecules and networks • Terms associated with bonded groups • bonds • e.g. harmonic, quartic • angles • e.g. harmonic, quartic • dihedral (torsion) angles • sin,cos (rotational barriers) • harmonic (e..g planarity constraints) • sometimes other cross-terms • bond-bond coupling • bond-angle coupling

15. Valence forcefields (ii) • Non-bonded terms • Summed over all non-covalently-bound pairs • always exclude bonded pairs • exclude 1,3 interactions for angles • sometimes scale 1,4 interactions for dihedrals • van der Waals • Buckingham, Lennard-Jones • electrostatics • simple coulomb (qiqj/r) • Need to decide the atomic charges … • distance-dependent dielectric • Approximate correction for solvation • Examples • MM2, AMBER, CHARMM, UFF, CFF, CVFF

16. Shell Model Force fields • Typically used for ionic solids • Leading terms are non-bonded • Electrostatics • often based on formal charges • polarisability of ions included by splitting total ion charge in • Core (often +ve) and Shell (-ve), modelling the valence electrons • Shell can shift in response to electrostatic forces, restoring forces from harmonic “spring” • van der Waals • sometimes compute using shell position • Can also incorporate 3-body terms • some bond angles are preferred over others, introducing some covalent character Core position Shell position

17. Choice of MM Model • 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

18. 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

19. 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)

20. 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)

21. 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

22. 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

23. 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

24. 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

25. 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

26. 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

27. 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

28. 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.

29. 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

30. 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

31. 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

32. 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)

33. 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

34. 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.

35. Computed Proton Affinities

36. 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

37. 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

38. 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

39. 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

40. 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

41. 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

42. Adjusted Connection Atoms • Semi-empirical parameterisation of boundary atom • Implemented in the MNDO package (Thiel el a) • No link atoms needed, boundary atom sited at MM centre • Typically boundary atom is C, parameterised to mimic CH3

43. Localised Orbital Approaches (i) • LSCF (Rivail et al) • Semi-empirical • Single orbital on QM boundary atom (pointing outwards) is frozen, based on calculation on a fragment (case-by-case set up) • GHO (Generalised Hybrid Orbital, Gao et al. • Semi-empirical (being extended to ab-initio) • Single orbital (sp3 hybrid) on MM boundary centre (pointing inwards) • Remaining 3 hybrid (“auxiliary orbitals”) are populated with fixed density matrix elements to produce correct MM charge • Semi-empirical parameters of the MM centre are adjusted based on model compounds (expected to be transferable)

44. Localised orbital Approaches (ii) • QSite implementation (Friesner et al) • Ab-initio implementation, in Jaguar package • Based on calculations on model fragments, using a particular basis set • Local orbitals include contribution from connected atoms (not just the QM and MM centres • Adjustment of MM parameters performed on a case-by-case basis, currently being used for protein system

45. 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)

46. 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

47. 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

48. 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

49. HDLCopt optimiser (I) hdlcopt • Hybrid Delocalised Internal Coordinate scheme (Alex Turner, Walter Thiel, Salomon Billeter) • Developed within QUASI project • O(N) overall scaling per step • Key elements: • Order N algorithm: • use delocalised internal coordinates, defined separately for each reside, O(N3 ) with residue, O(N) overall. • mix in cartesian coordinates to the internal coordinate generation (increases redundance but introduces dependent on translations and rotations • Iterate optimisation of each reside

50. HDLC Optimiser (ii) • Residue specification, often taken from a pdb file (pdb_to_res) allows separate delocalised coordinates to be generated for each residue • Can perform P-RFO TS search in the first residue with relaxation of the others • increased stability for TS searching • Much smaller Hessians • Further information on algorithm • S.R. Billeter, A.J. Turner and W. Thiel, PCCP 2000, 2, p 2177