An introduction to QM and QM/MM models Prof. Dr. Ville R. I. Kaila - - PowerPoint PPT Presentation
An introduction to QM and QM/MM models Prof. Dr. Ville R. I. Kaila Department of Chemistry Prof. Ville R. I. Kaila Technical University of Munich (TU Mnchen) CSC Spring School 16.3.2018 Outline of the lecture Introduction to QM and QM/MM
Department of Chemistry Technical University of Munich (TU München)
CSC Spring School 16.3.2018
Quantum mechanical models of (bio)chemical systems (QM, QM/MM, QM/QM methods) QM/MM models QM cluster models How to construct quantum (bio)chemical models? Some words about what "QM" can refer to in the QM/MM simulations Introduction to QM and QM/MM calculations
PES scans, QM/MM-FEP, LRA, QM/MM-US/WHAM and string simulations 3) Reaction pathways and free energies at QM/MM level 4) Characterizing the QM models by molecular property calculations Characterizing protein conformations by calculation of optical properties 1 & 2) Structures, energetics, and dynamics from QM and QM/MM What do I get out from QM/MM simulations?
Ubiquitin in water Setup of a protein QM/MM MD simulation
Modeling a SN2 reaction in water Exploring chemical reactions with QM/MM free energies
Understand how QM and QM/MM models work Become familiar with what type of information is obtained from biomolecular QM or QM/MM calculations
Know strategies to obtain potential energy surface scans from QM/MM Know methods to compute free energy from QM/MM simulations Become familiar with ways to characterize spectroscopic properties from QM/MM Become familiar with how to setup a QM/MM simulations for (bio)chemical systems How to compute free energies using QM/MM methods
The Kaila team at TU Munich
Enzyme Catalysis Energy converting enzymes Structure/Function/Dynamics Photobiology Bio-inspired catalysts Designing new proteins
Di Luca et al. PNAS (2017) Gamiz et al. JACS (2017) Supekar et al. Angew Chemie (2016) Sharma et al. PNAS (2015) Kaila et al. PNAS (2014) Mader et al. Nature Comm. (2018) Gamiz-Hernandez et al. JCPB (2014) Kaila & Hummer JACS 133 (2011) Kaila & Hummer PCCP 13 (2011) Suomivuori et al. PNAS (2017) Zhou et al. JACS (2014) Kaila et al. Nature Chemistry (2014) Gamiz et al. Angew. Chemie (2015)
How is light-energy captured? How to mimic nature? Differences in chem. vs. bio.? What is special about the biological environment?
To accurately describe (bio)chemical reactions one needs to rely on quantum chemical models
Example of applications: enzyme catalysis; photobiological systems; systems where biomolecular force fields fail (e.g. accurate differences in conformation of ligands, protein side chains)
Cut out a central part of the protein active site, Fix terminal atoms to account for protein framework model the surroundings as a polarizable medium Types of quantum chemical models used for (bio)chemical systems
Cut out a central part of the protein active site, Fix terminal atoms to account for protein framework model the surroundings as a polarizable medium Cut out a central part of the protein active site, model the surroundings explicitly using force fields, couple QM and MM regions together Types of quantum chemical models used for (bio)chemical systems
Cut out a central part of the protein active site, Fix terminal atoms to account for protein framework model the surroundings as a polarizable medium Cut out a central part of the protein active site, model the surroundings explicitly using force fields, couple QM and MM regions together
Cut out a central part of the protein active site, model the surroundings explicitly using approximate QM Account for the interaction between QM and QM regions Types of quantum chemical models used for (bio)chemical systems
Cut out a central part of the protein active site, Fix terminal atoms to account for protein framework model the surroundings as a polarizable medium Cut out a central part of the protein active site, model the surroundings explicitly using force fields, couple QM and MM regions together
Cut out a central part of the protein active site, model the surroundings explicitly using approximate QM Account for the interaction between QM and QM regions
Ab initio or density functional theory (DFT) treatment: N3-N6 (see lectures by Dr. Johansson)
Semi-empirical QM methods (e.g. AM1, MNDO, PM6/PM7, (SCC)-DFTB); linear scaling – N2 Reactive force fields (e.g. EVB methods); linear scaling It is challenging to apply any method with a formal computational scaling of higher than N4 in molecular simulations of chemical problems
DFT MP2 CCSD CCSD(T) FCI
Basis sets N5
HF
N4 N6 N7 N3 QM Method N2 N N!
Semi- empirical
Not-that-applicable for large biochemical systems Applicable for large biochemical systems
QM/MM discussion in these lectures will focus on DFT/MM but.....
Self Consistent Charge-DFTB model: (SCC-DFTB) Use point charges to describe the density Minimizing the energy of DFT eqn., but only with respect to the shape of the KS orbitals, not by changing the density sampling of 100-1000 ps range is possible with SCC-DFTB + accuracy well benchmarked for standard parametrizations
for all elements (Fe, S, etc.) SCC-DFTB developers:
Classical force fields have a defined topology – energy steeply rises if bonds are stretched
Empirical Valence Bonds (EVB) (Warshel and Weiss, 1981) X H Y 2 b1 b2 r3
Empirical Valence Bonds (EVB) methods have been to many biological proton transfer reactions
Correlated quantum chemistry Quantum Chemical DFT models QM/MM DFT/MM
QM MM
Classical Molecular Dynamics Accelerated sampling Continuum and coarse- grained method
Explore local chemistries Explore conformational phase-space
QM/QM
System size
Applicable timescale
Time scale Quantum chemical techniques Classical molecular simulation techniques
Gamiz et al. JACS (2017) Supekar et al. Angew. Chemie (2017) Di Luca et al. PNAS (2017) Suomivuori et al. PNAS (2017) Supekar et al. Angew. Chemie (2016) Sharma et al. PNAS (2015) Gamiz et al. Angew. Chemie (2015) Kaila et al. PNAS (2014) Zhou et al. JACS (2014) Kaila et al. Nature Chem. (2014)
Classical molecular dynamics simulations
input: experimental structure
different conformational/catalytic/non- equilibrium states
molecular structures
Quantum Classical
Classical molecular dynamics simulations
Hybrid quantum/classical simulations (QM/MM)
Free-energies/barriers for catalysis, mechanisms
Quantum Classical
coupling between local & non-local effects
site-directed mutagenesis biophysical experiments structural studies
Classical molecular dynamics simulations
molecular structures
Where should I cut the residues? Proteins are polymers
Where should I cut the residues? Where should I cut the residues? Cut at Cb atom minimal model where protein strain can easily be included Model residue size: 33 è 16 atoms
Where should I cut the residues? Where should I cut the residues? Cut at Cb atom minimal model where protein strain can easily be included Model residue size: 33 è 16 atoms
System of interest: enzyme active site/biochromophore Bioinformatical sequence comparisons can be highly informative in deciding, which residues to include in the QM model Cluster models: Include first protein solvation sphere, charged, hydrogen- bonding, stacking residues
Simplest environmental effect: Model the protein environment as a homogenous low-dielectric polarizable medium (e=4) Conductor-like Screening MOdel (COSMO) Polarizable Continuum Model (PCM) * * * * * * * Restrain or fix terminal carbons atoms in model to simulate the rigidity of the protein framework
Increasing the QM system size leads to saturation of the models
Schotte, Cho, Kaila et al. PNAS 109 (2012), Kaila et al. Nature Chemistry (2014) 124 atom model 142 atom model 176 atom model 250 atom model
DFT
X-ray density map
Explicit QM modeling of surroundings Extend the QM model by (conserved) second sphere ligands Saturation of cluster models?
Himo & Siegbahn, J. Biol. Inorg. Chem. 14, 636 (2009) but cf. also Liao & Thiel JCTC 8, (2012)
Frozen Density Embedding (FDET) Method is an alternative method to extend the size of the QM system, by dividing the system into an active system embedded in a frozen surrounding density
Wesolowski & Washel, JPC 97, (1993)
Ground state polarization of the surroundings can be considered by iterative freeze-and-thaw cycles
Wesolowski & Weber, Chem. Phys. Lett. 248, 71 (1996)
Although it’s possible to add many subsystems together, the protein framework has, nevertheless, to be cut implementations available in ADF and QChem
Frozen Density Embedding (FDET) methodology (Wesolowski & Warshel JPC, 97, 1993, Wesolowski et al. Chem Rev 2015)
Nonadditive exchange correlation component Embedded system Electron density
embedding system Embedding system (surrounding) Embedding potential Electron density
embedding system Kinetic energy component
there is also a quantum mechanical exchange between the active and surrounding systems
Hybrid quantum mechanics/molecular mechanics (QM/MM) – combine QM clusters with a classical environment description
Overview: Senn&Thiel Angew. Chemie 48 (2009)
MM treatment of the surroundings QM treatment of the active site in which covalent bonds can break/formation QM .... MM interaction
Forces
QM system from Hellmann-Feymann on the fly from the wavefunction
additive QM/MM this was MNDO2 considered "classically" standard biomolecular force field
Esurr The additive QM/MM energy expression:
Eactive = EDFT
non-bonded + link (covalent terms)
electrons-nuclei of QM, electrons-point charges of MM LJ-parameters assigned for QM atoms à this is called and electrostatic embedding and requires changes in electronic structure code
MM atoms L - Link atom QM atoms Link atoms are introduced to cap the QM system so that electrons do not flow from the unsaturated bond I – inner “active” QM system
C-C bond 1.54 Å C-H bond 1.09 Å a = 0.71 make boundary error smaller: charge shifting scheme NOTE: cut simple bonds Ca-Cb is ok, peptide bond is not recommended
In addition to link atoms frozen localized orbitals can also be employed to prevent electrons from flowing
Senn&Thiel
Accuracy of frozen localized orbitals is however similar as in the link atom approach NOTE: boundary region always introduces errors à move boundary as far as computationally feasible
The subtractive QM/MM energy expression: QM MM MM
In the original form this mechanical embedding scheme is called ONIOM (Our-own N-layer Integrated molecular Orbital and molecular Mechanics) by Morokuma and coworkers Benchmarking indicate mechanical embedding is not very accurate – but subtractive schemes with electrostatic embedding produce similar accuracies as additive schemes
Additive QM/MM Subtractive QM/MM
Mechanical embedding Electrostatic embedding Polarizable embedding (with polarizable force fields) +
The QM/MM methodology can also be extended to multiple QM regions
Röpke, Bärwinkel, Kaila MM QM
Michael Röpke
The QM/MM methodology can also be extended to multiple QM regions
Röpke, Bärwinkel, Kaila
QM1 QM2
MM QM
QM Flaig, Beer, Ochsenfeld JCTC 8 (2012)
Errors introduced in QM/MM partitioning? MM surroundings may help converging QM system
fixed QM Full optimization fixed QM Full
fixed QM Full optimization Fixed surr. Fixed QM
fixed surroundings
QM Cluster Models QM/MM Models
sufficient for convergence
convergence? Convergence is usually faster in QM/MM
to large surroundings. Requires fixing, micro- iterations etc.
comparison of total energies
Total energies can be difficult to interpret
terminal atoms to account for protein strain; capping by hydrogen atoms
capping with e.g. link atoms (over-polarization); hard to introduce additional continuum models
converge/study the effect of the surroundings; difficult to account for distant charged groups
explicitly included (protonation states?). ß perform in parallel! à
Saura et al. 2018
QM/MM methods are normally limited to order of magnitude slower sampling times than MD
How to overcome?
à Study the system first at MD level followed by QM/MM
Hydration dynamics on 300-1000 ns timescales
Di Luca, Gamiz-Hernandez, Kaila PNAS (2017)
NuoH
Di Luca, Gamiz, Kaila PNAS (2017) Kaila et al. PNAS (2014)
NuoL NuoM NuoN NuoH
The water wires support Grotthuss-type proton transfer
1) Structures and energies of intermediate states 2) Dynamics of a intermediates - specific events 3) Potential energy and/or free energy surfaces 4) Molecular property calculations: calculate optical or spectroscopic properties
Example: enzymatic reaction mechanisms of reaction 1 reaction 2
Example: How is quinone reduced in the active site of complex I? ToDo:
Can the process take place, if it is not observed in the QM/MM MD?
Gamiz-Hernandez, Jussupow, Johansson, Kaila, JACS (2017) Sharma, Belevich, Gamiz-Hernandez, Rog, Vattulainen, Verkhovskaya, Hummer, Kaila, PNAS (2015)
ToDo:
pathway optimizations = restrained optimizations, NEB, etc.
Tyr87 - Q His38 - Q Example: How is quinone reduced in the active site of complex I? QM/MM potential energy surface (PES) scans
∆EAV ‡ = −RTln 1 n # exp − ∆Ei ‡ RT n i=1
Saura-Martinez et al.
Example: C-H bond activation in mammalian 15-lipoxygenase-1 Sampling: 10 ns classical MD QM/MM PES on 20 structures picked from MD
Challenge: How to obtain representative PES in a highly fluctuating environment? Barrier computed from average Boltzmann weighted sum: variation based
but this converges rather slowly...
(e.g. Ryde JCTC 2017)
Back-and-forth sampling of the PES can also be used to help convergence... starts here scan this way scan that way until converged
In QM calculations the computational sampling time is normally limited to 1-50 ps Possible solution: do sampling using semi-empirical QM methods (SCC-DFTB, PM7, etc.) or reactive force fields (e.g. Empirical Valence Bond (EVB)), and compute the free energy of transferring the system into a first-principle QM/MM system QM/MM Free Energy Perturbation (QM/MM-FEP) (e.g. Warshel and coworker, Yang and coworker, Ryde and coworker) Perform a free energy perturbation by linear response: Perform a free energy perturbation calculation between reference and target surfaces:
Mlyńsky,́ et al.
QM(SCS-MP2)/MM QM(AM1)/MM QM(SCC-DFTB)/MM
Developed by Torrie and Valleau (J. Comput Phys 1977)
Introduces a bias potential to enhance sampling of the phase space The bias potential can be of any form – but in practice very often harmonic potentials are used
The restrains flatten the free energy surface
The biased runs introduce on-Boltzmann statistics - what would be unbiased probability to obtain taken the known biasing potential
Unbiased probability Biased probability The unbiased probability in terms of the biased probability
Using A = - 1/b ln p we get
Biasing potential Biased probability “constant” defining absolute height of PMF Unbiased free energy
The absolute height of the unbiased free energy profiles are not known – must be stitched together into a complete PMF profile
Developed by Torrie and Valleau (JCP 1954)
Method to minimize the statistical error of pu solved self-consistently
Phosphate bond break/formation r2-r1 (Å) Multi-dimensional reaction coordinate for all “bond breaking” - “bond forming” interatomic distances:
R = (r4-r3) + (r2-r1)
Associative path Dissociative path Reaction coordinate on bond formation and bond breaking reactions Proton transfer r4-r3 (Å) Water attacks gP Followed by water deprotonation Deprotonation
Followed by OH- attack on gP E33 example:
Proton transfer r4-r3 (Å) Phosphate bond break/formation r2-r1 (Å) Concerted path Concerted reaction mechanism E33 R380 R32 Mg2+ N37 ATP
Umbrella sampling applied to QM/MM followed by WHAM
R-H ... Y
r1 r2 place a bias on R = r1-r2 Ub = ½ k(R-R0)2 Example of QM/MM US
for an enzymatic process
Mader et al. (in preparation)
Free SVP TZVP Free energy (kcal mol-1) microsolvation calculated the quantum chemical way in implicit solvent H + TdS + ZPE calculated the QM/MM US in explicit solvent microsolvation the overlap needs to be good to obtain accurate free energies
QM/MM MD
excited state calculations from dipole moment autocorrelation function: ⎰ dt < µ(t) µ(0) > Compute IR spectra
Time-Dependent Density Functional Theory, TDDFT (N4 scaling) Solve: To obtain: transition vectors, XY excitation energies, w Coupled-cluster- based theories (e.g. second-order approximate coupled cluster theory, CC2)
Multireference calculations (CASSCF, CASPT2)
VEEs RVS-ADC(2)/def2-TZVP/CHARMM For each state: 5 ps QM/MM MD (B3LYP-D3/MM) with chromphore + nearby protein residues in QM; remaining system in MM 500 snapshots/state
Suomivuori, Gamiz-Hernandez, Sundholm, Kaila PNAS 2017
Three hierarchies of quantum chemical models
Additive & subtractive QM/MM schemes Mechanical, electrostatic, and polarizable embedding Link atoms or frozen localized orbitals are used to cap & saturate the QM region Difficult to properly minimize large QM/MM systems Micro-iterative schemes can help Dynamics, PES scans, Free energies from QM/MM Employ QM & QM/MM models in combination! Perform also classical dynamics to study long-timescale behavior of your QM/MM system