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 20 Potential Energy Surface - IV (Transition States) C7790 Introduction to Molecular Modelling -2- Context microworldmacroworld equilibrium (equilibrium constant) kinetics (rate constant) free energy (Gibbs/Helmholtz) partition function phenomenological thermodynamics statistical thermodynamics microstates (mechanical properties, E) states (thermodynamic properties, G, T,…) microstate ≠ microworld Description levels (model chemistry): • quantum mechanics • semiempirical methods • ab initio methods • post-HF methods • DFT methods • molecular mechanics • coarse-grained mechanics Structure EnergyFunction Simulations: • molecular dynamics • Monte Carlo simulations • docking • … C7790 Introduction to Molecular Modelling -3We only need to know the properties of individual components involved in the reaction at standard conditions (or at different conditions, which are well defined). Revision: Activation energy and modelling aA + bB cC+ dD 𝑘 RT G B e h Tk k   − = Eyring equation rate constant via TS A B TS easier for modelling Δ𝐺≠ = Δ𝐺𝑓,𝑇𝑆 − (𝑎Δ𝐺𝑓,𝐴 + 𝑏Δ𝐺𝑓,𝐵) activation free energy transition state C7790 Introduction to Molecular Modelling -4Revision: Partition function and modelling Helmholtz energy F: QTkF B ln−== − = K j E j eQ 1  Canonical partition function: Consider only the most important microstate approximation E1, E2, E3, … The most important microstate is the microstate with the lowest energy. 𝐹 = 𝐸1 It is very often used for qualitative consideration or when computationally demanding methods are employed (typically quantum chemical calculations). C7790 Introduction to Molecular Modelling -5Revision: Potential energy calculations QM (Quantum mechanics) MM (Molecular mechanics) CGM (Coarse-grained mechanics) )(RE )(RE )(RE R - position of atom nuclei R - position of atoms R - position of beads approximations approximations Potential energy surface can be calculated by various method (model chemistry)! C7790 Introduction to Molecular Modelling -6Revision: Transition State on PES E (x, y) x y local minima local maxima saddle points transition states C7790 Introduction to Molecular Modelling -7Elementary reaction and PES A+B C+D Elementary reaction is the transformation of reactants and products represented by prereaction and post-reaction complexes separated by just one transition state. TS states (reaction coordinate) pre-RC post-RC individual non-interacting components at standard state stationary points on PES pre-RC (reaction complex) TS post-RC reaction pathway potentialenergy potentialenergy C7790 Introduction to Molecular Modelling -8Finding Transition States on PES If the initial geometry is already close to searched transition state: ➢ local optimalization methods If the reactants, products or both are known: ➢ semi-global optimization methods ▪ Single Coordinate Driving Method ▪ Synchronous Transit Methods ▪ Nudged Elastic Band Methods C7790 Introduction to Molecular Modelling -9Reaction Coordinate ➢ The reaction coordinate describes the geometrical changes, which are necessary for transforming reactants into products via the transition state and vice versa. ➢ Its correct form is not a priori known. ➢ Approximate reaction pathways can be obtained by: ➢ nudged elastic band method ➢ string method ➢ If the transition state is known, then the intrinsic reaction coordinate (IRC) can be calculated. ➢ IRC describes geometrical changes resulting from going "downhill" from TS in both directions in such a way that the energy gradient in all perpendicular directions is zero. Since the potential energy is a state function, the exact knowledge of reaction pathway is not needed. Only proper characterization of states (reactants, products, transition states) is necessary. However, at least some crude representation of reaction pathway is necessary to find transition state. C7790 Introduction to Molecular Modelling -10Local Methods Locating Transition Structures C7790 Introduction to Molecular Modelling -11Local Methods LM geometry optimizers: ➢ move geometry "downhill" only TS geometry optimizers: ➢ move geometry "uphill" along negative eigenvector ➢ at the same time, move geometry "downhill" along the other degrees of freedom ➢ Transition state (TS) is a saddle point on PES. ➢ The saddle point is a stationary point: ➢ PES curvature at the saddle point provides one negative Hessian H eigenvalue: ➢ Searching for a transition state structure is more complex procedure than searching for a local minima because: 𝜕𝐸(𝐑) 𝜕𝐑 = 0 𝐇𝐜 𝑘 = 𝜆 𝑘 𝐜 𝑘 all positive lk - local minimum (LM) one negative, other positive lk - first order saddle point transition state - TS TS LM geometry energy C7790 Introduction to Molecular Modelling -12Local Methods, cont. Transition state geometry optimizers: ➢ pseudo-second order methods (quasi-Newton methods) need to be employed ➢ starting geometry (structure) must be very close to the saddle point ➢ initial Hessian must be of high quality ➢ empirical Hessians employed in local geometry optimization are not usable ➢ Hessian calculated at the same level theory as used during geometry optimization is good choice. However, this is computationally very demanding. ➢ for difficult cases, Hessian need to be recalculated at the same level of theory at each geometry optimization step (quasi-Newton -> full Newton method). This is extremely computationally demanding and suitable only for small molecular systems. Starting structure for SP (TS) search can be found by pseudo-global optimalization methods: ➢ single coordinate driving method ➢ linear and quadratic synchronous transit methods ➢ nudged elastic methods ➢ string methods ➢ others …. C7790 Introduction to Molecular Modelling -13Local Methods - Practical Realizations %Chk=ts.chk # PM3 Opt(CalcFC,TS,NoEigenTest,MaxCycle=25) transition state optimization XX X .... ..... ..... request optimization of the transition state CalcFC force constants (Hessian) are calculated on the input geometry CalcAll force constants (Hessian) are calculated on each geometry during optimization (full Newton-Raphson method) Transition structure search in Gaussian do not waste computational time on bad structures skip very tight requirements on starting structures if the starting structure is not of good quality, the optimization can fail C7790 Introduction to Molecular Modelling -14Validation of TS ➢ It is necessary to validate the nature of stationary state representing the transition state employing vibrational analysis. ➢ Only ONE normal mode must be imaginary, the other normal modes must be represented by frequencies, which are real numbers. ➢ Atom movement along imaginary vibrational mode must follow the intended reaction change (formation/breaking bonds). C7790 Introduction to Molecular Modelling -15Semi-global Methods Single Coordinate Driving (SCD) Method Synchronous Transit (ST) Methods Nudged Elastic Band (NEB) Methods Locating Transition Structures C7790 Introduction to Molecular Modelling -16Single Coordinate Driving Method Single coordinate driving (SCD) method is a simple approach for finding initial structures for the transition state geometry optimization. SCD is a special case of potential energy surface scan methods (limited to a single geometry parameter). Driven geometry parameter should represent the geometry change occurring during the reaction. Suitable coordinates can be: ➢ distance of formed bond ➢ distance of broken bond ➢ dihedral angle representing rotation around bond during a conformational change ➢ others … ➢ The SCD method can suffer from driving hysteresis. Two different driving pathways can be observed when driving from the reactant or product state. C7790 Introduction to Molecular Modelling -17Linear synchronous transit (LST) method ➢ Geometry parameters (Cartesian or internal coordinates) vary linearly between initial (reactants) and final (products) states: ➢ The order of atoms in 𝑹 𝑟𝑒𝑎𝑐𝑡𝑎𝑛𝑡𝑠 and 𝑹 𝑝𝑟𝑜𝑑𝑢𝑐𝑡𝑠 must be the same! ➢ The guess of transition state geometry is taken as a geometry with the maximum energy. ➢ The structure with maximum energy is usually far away from the correct transition state. ➢ Thus, the geometry obtained by LST must be optimized using a local method. Synchronous Transit Methods 𝑹𝑖 = 𝑁 − 𝑖 𝑁 𝑹 𝑟𝑒𝑎𝑐𝑡𝑎𝑛𝑡𝑠 + 𝑖 𝑁 𝑹 𝑝𝑟𝑜𝑑𝑢𝑐𝑡𝑠 i - 0, 1, 2, …, N Improved versions: ➢ QST2 - quadratic synchronous transit between reactants and products ➢ QST3 - quadratic synchronous transit between reactants and products via crude approximation of transition state C7790 Introduction to Molecular Modelling -18Nudged Elastic Band Method Nudged elastic band (NEB) method is an approach for finding reaction paths and transition structures. ➢ Reaction path is represented by beads. Each bead represent a specific system geometry. ➢ NEB ensures equidistant spacing between beads by springs. ➢ Potential energy is minimized in all perpendicular degrees of freedom. ➢ Due to path discrimination, only guess of transition state geometry is obtained, which need to be refined by local optimization. https://www.scm.com/doc/AMS/Tasks/NEB.html The string method (SM) is similar as NEB. SM uses splines to achieve equidistant bead separation. C7790 Introduction to Molecular Modelling -19- Summary ➢ Finding transition states is necessary step for quantification of reaction kinetics. ➢ Broad range of methods is available to locate either correct TS or at least its estimate. ➢ In comparison with characterization of reactants and products, finding transition structures is much difficult task. ➢ Nature of found transition states needs to be verified by vibrational analysis. ➢ Transition state must be a stationary state with only one imaginary molecular vibrations.