C7790 Introduction to Molecular Modelling -1C7790 Introduction to Molecular Modelling TSM Modelling Molecular Structures Petr Kulhánek kulhanek@chemi.muni.cz National Centre for Biomolecular Research, Faculty of Science Masaryk University, Kamenice 5, CZ-62500 Brno PS/2020 Distant Form of Teaching: Rev1 Lesson 27 Various (Molecular Dynamics II) C7790 Introduction to Molecular Modelling -2- Contents ➢ Free Energy Calculations ➢ Sampling Problem ➢ Potential of Mean Force Methods ➢ Metadynamics ➢ Constrained Dynamics ➢ Adaptive Biasing Force ➢ Multiple Walkers Approach ➢ Molecular Dynamics and Reactions ➢ Empirical Valence Bond ➢ CPMD versus BOMD -2- C7790 Introduction to Molecular Modelling -3Free Energy Calculations C7790 Introduction to Molecular Modelling -4Free Energy KRTAr ln−= RT A B e h Tk k   − = 1 rA  A Free energy is related to equilibrium and rate constants. Knowledge of free energy allows to quantify: • chemical reactivity (e.g. enzymatic activities) • thermodynamics (e.g. binding affinities) C7790 Introduction to Molecular Modelling -5Free Energy Calculations 9th November 2010, Bratislava -5- 1 2 21 ln   RTA −= → 0)(ln)( ARTA +−=  Density of state (probability) It can be calculated from molecular dynamics or Monte Carlo simulations, but … reaction coordinate collective variables C7790 Introduction to Molecular Modelling -6Sampling Problem d1 10 ns long simulation is able to discover free energy landscape with depth only about 3 kcal/mol. C7790 Introduction to Molecular Modelling -7Free Energy Calculations Available methods: ➢constrained dynamics system is biased by constraining reaction coordinate ➢adaptive biasing force system is biased by force which is opposite to potential of mean force ➢umbrella sampling system is biased by restraining reaction coordinate ➢metadynamics system is biased by Gaussian hills, which fill free energy landscape A system has to be biased achieving efficient sampling in the region of interest. We need to know how to obtain the unbiased free energy from such biased simulation. C7790 Introduction to Molecular Modelling -8Free Energy Calculations ➢Alchemical Transformation one system is slowly changed to another one (changes are very often unrealistic, atoms are created and/or annihilated) what: mostly changes in binding free energies: how: thermodynamic integration (TI), free energy perturbation (FEP) ➢ Potential of Mean Force system is changed along reaction coordinate what: free energy of conformation changes, reaction free energies how: constrained dynamics, adaptive biasing force, umbrella sampling, metadynamics, steered dynamics ➢ End-points Methods free energy of every state is calculated independently what: mostly binding free energies how: MM/XXSA; XX=PB, GB, LRA C7790 Introduction to Molecular Modelling -9- Metadynamics C7790 Introduction to Molecular Modelling -10Metadynamics, theory ( ) ( ) i i i x xV = dt txd m   −2 2 Free energy landscape is filled by Gaussian hills. ( ) ( ) ( ) ix,V+xV x = dt txd m h i i i   −2 2 Equations of motion MTD history potential Equations of MTD motion (direct approach) ( ) ( )         − − 2 σ ss H=is,V t i =t th 2 exp 2 1 History-dependent term converges to FES ( ) ( )si,V=sA h i − → lim C7790 Introduction to Molecular Modelling -11Metadynamics, example DIS (distance) hill height 0.01 kcal/mol, width 0.5 x 0.5 Å MTD frequency 500 fs 2 ns long simulation 300 K, vacuum, GAFF force field, time step 0.5 fs C7790 Introduction to Molecular Modelling -12Constrained Dynamics C7790 Introduction to Molecular Modelling -13Constrained Dynamics, theory ( ) ( ) i i i x xV = dt txd m   −2 2 Reaction coordinate is fixed (constrained) at the value of interest. ( ) ( ) ( ) ( )         −  xσtλ+xV x = dt txd m k k k i i i 2 2 ( ) ( ) 00 =ξxξ=xσ − Equations of motion Constraint condition method of Lagrange multipliers Equations of constrained motion = Lagrange multipliers holonomic constraint kλ C7790 Introduction to Molecular Modelling -14Constrained Dynamics, theory ( )dξ ξ dξ ξdA =ΔG ξ uc  2 1 Derivative of unbiased free energy is also given by (concise formulation): ( ) ( ) ( ) 0 2/1 0 ln ξξ ucccuc Z dξ d RTλ= dξ ξdA + dξ ξdA = dξ ξdA − → −− Final free energy is obtained by numerical integration: second derivatives are not required C7790 Introduction to Molecular Modelling -15Constrained Dynamics, example Small numbers are calculated by averaging big numbers. 70 points d177 points,  0.1 Å Method B: 5 ps shift, 5 ps equilibration, 20 ps production 300 K, vacuum, GAFF force field, time step 1 fs / 0.5 fs DIS (distance) C7790 Introduction to Molecular Modelling -16Constrained Dynamics, example d177 points,  0.1 Å Method B: 5 ps shift, 5 ps equilibration, 20 ps production 300 K, vacuum, GAFF force field, time step 1 fs / 0.5 fs DIS (distance) C7790 Introduction to Molecular Modelling -17Adaptive Biasing Force C7790 Introduction to Molecular Modelling -18ABF, theory ( ) ( ) i i i x xV = dt txd m   −2 2 Movement along reaction coordinate is the subject of diffusion process. Equations of motion Free energy and force along RC Equations of ABF motion ( ) ( ) ( )ξF+ x xV = dt txd m ABF i i i   −2 2 force along reaction coordinate is subtracted from the system ( ) ( ) dx dξ dξ ξdA =xFABF − C7790 Introduction to Molecular Modelling -19ABF, theory ξ ξ dt dξ Zdt d = ξ A         −   1 Free energy is given by: it contains the second derivatives of reaction coordinate if treated analytically equation is solved numerically Procedure: • range of reaction coordinate is divided into bins • a value of reaction coordinate determine a bin • a contribution to the derivative of free energy is accumulated into a bin • ABF force calculated from accumulated free energy derivative is applied to the system • accumulated free energy derivative very rapidly converges C7790 Introduction to Molecular Modelling -20ABF, example DD (difference of distances)  range from -6.5 to 6.5 Å, 130 bins 2 ns long simulation 300 K, vacuum, GAFF force field, time step 1 fs Show movie? C7790 Introduction to Molecular Modelling -21Multiple Walkers Approach Server Client 1 Client 2 Client 3 server collects information about free energy surface (FES) and redistributes it among clients Applicable to: ➢Metadynamics ➢Adaptive biasing force Advantages: ➢ Faster convergence ➢ Easy to implement ➢ Parallel scaling is almost linear FES data are exchanged every Nupdate MD steps Nupdate ~ 500-1000 steps C7790 Introduction to Molecular Modelling -22Multiple Walkers Approach Two walkers (from reactants and products) One walker (from reactants) Nucleophilic substitution reaction (test case) C7790 Introduction to Molecular Modelling -23Molecular Dynamics and Reactions 9th November 2010, Bratislava -23- C7790 Introduction to Molecular Modelling -24Description of Chemical Reactions EVB X QM/MM X QM theory precision complexity, sampling problems too complex for biomolecular systems problem with boundaries between QM and MM parts requires parametrizations QM MM QM products reactants EVB C7790 Introduction to Molecular Modelling -25- CPMD Car-Parrinello Molecular Dynamics 9th November 2010, Bratislava -25- C7790 Introduction to Molecular Modelling -26- Method Equations of motion: fictitious mass of wavefunction (ca 300-700 a.u. , typical value is about 600 a.u. ions wavefunction constraints due to orthonormality of wavefunction C7790 Introduction to Molecular Modelling -27CPMD versus BOMD CPMD • no SCF procedure • motion of ions in time • motions of wavefunction in time • time step ~ 0.1 fs (5 a.u.) • DFT only (in CPMD) • hybrid functional possible but very slow • dispersion correction available • planewaves wavefunction (periodicity!} • wavefunction quality is determined by cutoff (single value) • pseudopotentials required (core electrons) BOMD • SCF procedure • motion of ions in time only • wavefunction follows ions by SCF (BO) • time step max ~ 1 fs • gradients require very tight convergence of wavefunction optimization