Numerical integration schemes for the equations of motion, constants - PowerPoint PPT Presentation

Numerical integration schemes for the equations of motion, constants of motion, control of stability, accuracy Ari Paavo Seitsonen École Normale Supérieure CECAM QM/MM School, Lausanne April 8 th -12 th , 2019 apsi QM/MM 2019 1 / 53

Molecular dynamics apsi QM/MM 2019 2 / 53

Molecular dynamics Propagation of Newton’s equation of motion (with discrete equations of motion) F I = M I a = M I ¨ R I Alternative derivation from the Lagrange formalism: N 1 � R N � R N , ˙ R N � � 2 M I ˙ R 2 I − U � L = , I = 1 U is the interaction potential between the particles. The Euler-Lagrange equation d ∂ L = ∂ L dt ∂ ˙ ∂ R I R I Most common algorithm: Verlet algorithm (in a few variations) apsi QM/MM 2019 3 / 53

Verlet algorithm: Velocity Verlet discretisation of Newton’s equation of motion M I ¨ R I = F I i) Propagate ionic positions R I ( t ) according to R I ( t + ∆ t ) = R I ( t ) + ∆ t v I ( t ) + (∆ t ) 2 F I ( t ) 2 M I ii) Evaluate forces F I ( t + ∆ t ) at R I ( t + ∆ t ) iii) Update velocities v I ( t + ∆ t ) = v I ( t ) + ∆ t [ F I ( t ) + F I ( t + ∆ t )] 2 M I apsi QM/MM 2019 4 / 53

Velocity Verlet: Advantages Other algorithms provides can have better short time stability and allow larger time steps, but . . . simple and efficient; needs only forces, no higher energy derivatives still correct up to and including third order, (∆ t ) 3 explicitly time reversible sympletic: conserves volume in phase space superior long time stability (energy conservation) of the Verlet algorithm apsi QM/MM 2019 5 / 53

Velocity Verlet: Choice of time step The time step is in general chosen as large as possible . . . “Possible” = stable dynamics = energy conserved; or, drift in energy acceptable Rule of thumb: 6-10 times smaller than the fastest period in the system; otherwise sampling of that mode is impossible Time step can be changed during simulation(!) apsi QM/MM 2019 6 / 53

Velocity Verlet: Choice of time step AlCl 3 dimer Example of a good/bad choice of time step Highest vibrational frequency 595 cm − 1 ⇒ period T = 56 fs Divergence between δ t = 400 .. 500 atu = 9.6-12.0 fs ≈ 1/5 T apsi QM/MM 2019 7 / 53

Equations of motion: Alternative derivation Propagation methods Define phase space vector Γ = ( x , p ) and commutator { A , H } = ∂ A ∂ H ∂ p − ∂ A ∂ H ∂ x ∂ p ∂ x Hamilton’s equations of motion: d Γ dt = { Γ , H } Define ˆ L so that i ˆ L Γ = { Γ , H } Γ = i ˆ ˙ L Γ ⇒ Γ( t ) = e i ˆ L t Γ( 0 ) Such formalism has been used by Mark Tuckerman et al to derive new integrators apsi QM/MM 2019 8 / 53

Tricks Simulated annealing Multiple time scales / RESPA Periodic boundary conditions Ewald summation Thermodynamic integration Cell lists etc apsi QM/MM 2019 9 / 53

Molecular dynamics: Summary Molecular dynamics can be used to perform real-time dynamics in atomistic systems Maximum time step ∆ t ≈ 1 fs (highest ionic frequency 2000 − 3000 cm − 1 ) Temperature can be controlled via rescaling – (initial) equilibration – and thermostats ( e. g. Nosé-Hoover thermostat chains) for NVT ensemble Constraints can be used to pose restrictions on the atoms They can be used to direct reactions, however in complicated (potential/free) energy landscapes they might not yield the correct reaction path (in reasonable simulation time, at least) Metadynamics looks like a promising method for finding reaction paths and (potential/free) energy surfaces apsi QM/MM 2019 10 / 53

Ab initio molecular dynamics apsi QM/MM 2019 11 / 53

Realistic MD simulations M I ¨ R I = −∇ R E ( { R J } ) Classical molecular dynamics: E ( { R J } ) given e. g. by pair potentials How about estimating E ( { R J } ) directly from electronic structure method? What is needed is −∇ R E ( { R J } ) = − dE d R I apsi QM/MM 2019 12 / 53

Classical vs MD simulations When is electronic structure needed explicitly, when is classical treatment sufficient? ◮ Chemical reactions: Breaking and creation of chemical bonds ◮ Changing coordination ◮ Changing type of interaction ◮ Difficult chemistry of elements Combination of both: QM/MM apsi QM/MM 2019 13 / 53

Born-Oppenheimer molecular dynamics apsi QM/MM 2019 14 / 53

Born-Oppenheimer Ansatz Separate the total wave function to quickly varying electronic and slowly varying ionic wave function: N BO � ˜ Φ BO ( { r i } , { R I } ; t ) = χ ( { R I } ; t ) Ψ k ( { r i } , { R I } ) ˜ k = 0 Leads to a Schrödinger-like equation for the electrons and a Newton-like equation for the ions (after some assumptions for the ionic wave function): H e ˜ E e { R I } ˜ Ψ k ( { r i } , { R I } ) = Ψ k ( { r i } , { R I } ) M I ¨ R I = F I Electrons always at the ground state when observed by the ions Usually valid, however there are several cases when this Ansatz fails apsi QM/MM 2019 15 / 53

Born-Oppenheimer MD Lagrangean N 1 � � { ψ i } , R N � R , ˙ � 2 M I ˙ R 2 { ψ i } E KS � L BO R = I − min I = 1 equations of motion: = − d � { ψ i } , R N �� Ψ , R N �� M I ¨ � E KS � { ψ i } E KS � R I = −∇ R min d R I If the right-hand side can be evaluated analytically it can be plugged directly to the Verlet algorithm apsi QM/MM 2019 16 / 53

Forces in BOMD What is needed is − d � { ψ i } , R N �� { ψ i } E KS � min d R I with the constraint that the orbitals remains orthonormal; this is achieved using Lagrange multipliers in the Lagrangean E KS = E KS + � Λ ij ( � ψ i | ψ j � − δ ij ) ij Forces   d E KS = ∂ E KS  ∂ E KS ∂ ∂ � ψ i | � � � + Λ ij � ψ i | ψ j � + ∂ � ψ i | + Λ ij | ψ j � d R I  ∂ R I ∂ R I ∂ R I ij ij j When | ψ i � optimal F KS ( R I ) = − ∂ E KS ∂ � + Λ ij � ψ i | ψ j � ∂ R I ∂ R I ij apsi QM/MM 2019 17 / 53

BOMD: Error in forces The error in the forces depends on the convergence criterion set for the electronic structure in BOMD: apsi QM/MM 2019 18 / 53

BOMD: Observations The energy needs to be minimal in order to estimate the forces The accuracy of the forces depends on the level of self-consistency Thus a competition between accuracy and computational cost Constant of motion: ◮ NVE: N 1 { ψ i } , R N � � 2 M I ˙ R 2 { ψ i } E KS � I + min I = 1 ◮ NVT: N 1 + 1 s 2 + gk B T ln( s ) � 2 s 2 M I ˙ { ψ i } , R N � R 2 { ψ i } E KS 2 Q s ˙ � I + min I = 1 apsi QM/MM 2019 19 / 53

Car-Parrinello method apsi QM/MM 2019 20 / 53

Car-Parrinello method Roberto Car & Michele Parrinello, Physical Review Letters 55 , 2471 (1985) They postulated Langangean M 1 � � � � � � ˙ ; R , ˙ � ψ i | ˙ ˙ L CP { ψ i } , ψ i R = 2 µ ψ i i = 1 { ψ i } , R N � − min { ψ i } E KS � N 1 � 2 M I ˙ R 2 + I I = 1 Reminder: E KS contains the Lagrange multipliers for orthonormality of orbitals Fictitious or fake dynamics of electrons µ = fictitious mass or inertia parametre Simultaneous dynamics of ions and electrons apsi QM/MM 2019 21 / 53

Car-Parrinello method: Equations of motion Euler-Lagrange equations d ∂ L CP ∂ L CP = dt � � ∂ � ψ i | ˙ ∂ ψ i � � d ∂ L CP ∂ L CP = dt � � ∂ � R I | ˙ ∂ R I � � Equations of motion − ∂ E KS µ ¨ � ψ i = ∂ � ψ i | + Λ ij | ψ j � j − ∂ E KS ∂ M I ¨ � R I = + Λ ij � ψ i | ψ j � ∂ R I ∂ R I ij apsi QM/MM 2019 22 / 53

Car-Parrinello method: Simultaneous dynamics Unified Approach for Molecular Dynamics and Density-Functional Theory Electronic and ionic structure evolve simultaneously Whereas in BOMD first the electronic structure is optimised, then the ions are moved apsi QM/MM 2019 23 / 53

Car-Parrinello method: Constant of motion Constant of motion M N 1 1 � � � ψ i | ˙ ˙ { ψ i } , R N � � 2 M I ˙ E conserved = + E KS R 2 � 2 µ ψ i + I i = 1 I = 1 { ψ i } , R N � Note: instantaneous value of E KS � , not minimum Thus no need to optimise the orbitals at each step apsi QM/MM 2019 24 / 53

Magic Car-Parrinello method Does the Car-Parrinello method yield physical results even if the orbitals are not at the Born-Oppenheimer surface? ◮ Yes — provided that the electronic and ionic degrees of freedom remain adiabatically separated and the electrons close to the Born-Oppenheimer surface ◮ Why? — dynamics of the electrons is artificial, or unphysical and thus has to average out during the time scale of ionic movement Another way of viewing: The electrons are slightly above the BO surface but remain there and average out the effects on the ions (to be considered with care) apsi QM/MM 2019 25 / 53

Adiabatic separation Pastore, Smargiassi & Buda, PRA 1991 Vibrational spectra of electrons and ions do not overlap: Triangle = highest ionic frequency � ∞ apsi QM/MM 2019 26 / 53 f e ( ω ) = � � � � ˙ � ˙ cos ( ω t ) ψ ( t ) dt ψ ( 0 ) �

Adiabatic separation Thus there’s no efficient mechanism for exchange of energies: The two subsystems are adiabatically decoupled Triangle = highest ionic frequency � ∞ f e ( ω ) = � � � cos ( ω t ) � ψ i ( t ) ˙ � ˙ dt ψ i ( 0 ) � t = 0 i apsi QM/MM 2019 27 / 53