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 Suprieure CECAM QM/MM School, Lausanne April 8 th -12 th , 2019 apsi QM/MM 2019 1 / 53
apsi QM/MM 2019 1 / 53
apsi QM/MM 2019 2 / 53
Propagation of Newton’s equation of motion (with discrete equations of motion) FI = MIa = MI ¨ RI Alternative derivation from the Lagrange formalism: L
RN =
N
1 2MI ˙ R2
I − U
, U is the interaction potential between the particles. The Euler-Lagrange equation d dt ∂L ∂ ˙ RI = ∂L ∂RI Most common algorithm: Verlet algorithm (in a few variations)
apsi QM/MM 2019 3 / 53
discretisation of Newton’s equation of motion MI ¨ RI = FI i) Propagate ionic positions RI(t) according to RI(t + ∆t) = RI(t) + ∆t vI(t) + (∆t)2 2MI FI(t) ii) Evaluate forces FI(t + ∆t) at RI(t + ∆t) iii) Update velocities vI(t + ∆t) = vI(t) + ∆t 2MI [FI(t) + FI(t + ∆t)]
apsi QM/MM 2019 4 / 53
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
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;
Time step can be changed during simulation(!)
apsi QM/MM 2019 6 / 53
AlCl3 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
Propagation methods Define phase space vector Γ = (x, p) and commutator {A, H} = ∂A ∂x ∂H ∂p − ∂A ∂p ∂H ∂x Hamilton’s equations of motion: dΓ dt = {Γ, H} Define ˆ L so that i ˆ LΓ = {Γ, H} ˙ Γ = i ˆ LΓ ⇒ Γ(t) = ei ˆ
LtΓ(0)
Such formalism has been used by Mark Tuckerman et al to derive new integrators
apsi QM/MM 2019 8 / 53
Simulated annealing Multiple time scales / RESPA Periodic boundary conditions Ewald summation Thermodynamic integration Cell lists etc
apsi QM/MM 2019 9 / 53
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
apsi QM/MM 2019 11 / 53
MI ¨ RI = −∇RE ({RJ}) Classical molecular dynamics: E ({RJ}) given e. g. by pair potentials How about estimating E ({RJ}) directly from electronic structure method? What is needed is −∇RE ({RJ}) = − dE
dRI
apsi QM/MM 2019 12 / 53
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
apsi QM/MM 2019 14 / 53
Separate the total wave function to quickly varying electronic and slowly varying ionic wave function: ΦBO ({ri} , {RI} ; t) =
NBO
˜ Ψk ({ri} , {RI}) ˜ χ ({RI} ; t) 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): He ˜ Ψk ({ri} , {RI}) = Ee
{RI} ˜
Ψk ({ri} , {RI}) MI ¨ RI = FI 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
Lagrangean LBO
R
N
1 2MI ˙ R2
I − min {ψi} EKS
equations of motion: MI ¨ RI = −∇R
= − d dRI
{ψi} EKS
If the right-hand side can be evaluated analytically it can be plugged directly to the Verlet algorithm
apsi QM/MM 2019 16 / 53
What is needed is − d dRI
{ψi} EKS
with the constraint that the orbitals remains orthonormal; this is achieved using Lagrange multipliers in the Lagrangean EKS = EKS +
Λij (ψi| ψj − δij) Forces dEKS dRI = ∂EKS ∂RI +
Λij ∂ ∂RI ψi| ψj +
∂ ψi| ∂RI ∂EKS ∂ ψi| +
Λij |ψj When |ψi optimal FKS (RI) = −∂EKS ∂RI +
Λij ∂ ∂RI ψi| ψj
apsi QM/MM 2019 17 / 53
The error in the forces depends on the convergence criterion set for the electronic structure in BOMD:
apsi QM/MM 2019 18 / 53
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 2MI ˙ R2
I + min {ψi} EKS
◮ NVT:
N
1 2s2 MI ˙ R2
I + min {ψi} EKS
+ 1 2Qs ˙ s2 + gkBT ln(s)
apsi QM/MM 2019 19 / 53
apsi QM/MM 2019 20 / 53
Roberto Car & Michele Parrinello, Physical Review Letters 55, 2471 (1985) They postulated Langangean LCP
ψi
R
M
1 2µ
ψi | ˙ ψi
+
N
1 2MI ˙ R2
I
Reminder: EKS contains the Lagrange multipliers for orthonormality of
Fictitious or fake dynamics of electrons µ = fictitious mass or inertia parametre Simultaneous dynamics of ions and electrons
apsi QM/MM 2019 21 / 53
Euler-Lagrange equations d dt ∂LCP ∂
ψi
∂LCP ∂ ψi| d dt ∂LCP ∂
RI
∂LCP ∂ RI| Equations of motion µ ¨ ψi = − ∂EKS ∂ ψi| +
Λij |ψj MI ¨ RI = −∂EKS ∂RI +
Λij ∂ ∂RI ψi| ψj
apsi QM/MM 2019 22 / 53
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
Constant of motion Econserved =
M
1 2µ
ψi | ˙ ψi
+
N
1 2MI ˙ R2
I
Note: instantaneous value of EKS
, not minimum Thus no need to optimise the orbitals at each step
apsi QM/MM 2019 24 / 53
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
Pastore, Smargiassi & Buda, PRA 1991 Vibrational spectra of electrons and ions do not overlap: Triangle = highest ionic frequency f e (ω) = ∞ cos (ωt) ˙ ψ (t)
ψ (0)
apsi QM/MM 2019 26 / 53
Thus there’s no efficient mechanism for exchange of energies: The two subsystems are adiabatically decoupled Triangle = highest ionic frequency f e (ω) = ∞
t=0
cos (ωt)
ψi (t)
ψi (0)
apsi QM/MM 2019 27 / 53
Physical and conserved energy: Ephysical = EKS
+
N
1 2MI ˙ R2
I
Econserved =
M
1 2µ
ψi | ˙ ψi
+
N
1 2MI ˙ R2
I
= Ekin,fict + Ephysical The difference, Ekin,fict = M
i=1 1 2µ
ψi | ˙ ψi
changes in the physical energy
apsi QM/MM 2019 28 / 53
Model system: Two-atom Si-fcc Energy components Ekin,f
apsi QM/MM 2019 29 / 53
Deviation of forces in CP dynamics from the true BO forces small and/but oscillating Fx(Si) Fx,CP(Si)-Fx,BO(Si)
apsi QM/MM 2019 30 / 53
Harmonic analysis: ωe
ij =
µ εi occupied, εj unoccupied (virtual) orbitals Lowest frequency ωe
min ∝
µ Highest frequency ωe
max ∝
µ Thus maximum possible time step (∆te)max ∝ µ Ecut
apsi QM/MM 2019 31 / 53
Lowest frequency has to be well above ionic frequencies ωe
min ∝
µ Highest frequency limits the maximum possible time step ωe
max ∝
µ (∆te)max ∝ µ Ecut If ∆t fixed and µ chosen
◮ too small: Electrons too light and adiabacity will be lost ◮ too large: Electrons too heavy, the slowest electronic motion starts
to overlap with the ionic frequencies and adiabacity will be lost
apsi QM/MM 2019 32 / 53
Vacancy in hot 64-atom Si cell
apsi QM/MM 2019 33 / 53
Sn2: Degeneracy of HOMO and LUMO at short distances
apsi QM/MM 2019 34 / 53
Two-level, two-electron model Wave function ψ =
2
2
θ is the electronic degree of freedom Constant gap Opening-closing gap G
apsi QM/MM 2019 35 / 53
One way to (try to) overcome the problem in coupling of electronic and ionic dynamics is to thermostat also the electrons [Blöchl & Parrinello, PRB 1992] Thus electrons cannot heat up; if they try to, thermostat will adsorb the excess heat Target fictitious kinetic energy Ekin,0 instead of temperature “Mass” of thermostat to be selected appropriately:
◮ Too light: Adiabacity violated (electrons may heat up) ◮ Too heavy: Ions dragged excessively
Please remember: The conserved quantity changed
apsi QM/MM 2019 36 / 53
Example: Aluminium Dependence of the heat transfer on the choice of Ekin,0 in solid Al
apsi QM/MM 2019 37 / 53
64 atoms of molten aluminium (a): Without thermostat (b): With thermostat
apsi QM/MM 2019 38 / 53
Check: Radial pair correlation function
◮ Solid line: CP-MD with thermostat ◮ Dashed line: Calculations by Jacucci et al apsi QM/MM 2019 39 / 53
The fictitious electronic mass exerts an extra “mass” on the ions and thereby modifies the equations of motion: MI ¨ RI = FI + µ
¨ RI ∂φi ∂r ∂φi ∂r The new equations of motion: (MI + dMI) ¨ RI = FI where dMI = 2 3µEI
kin
is an unphysical “mass”, or drag, due to the fictitious kinetics of the electrons
Example: Vibrations in water molecule mode harmonic BOMD 50 100 200 400 dM/M [%] bend 1548 1543 1539 1535 1529 1514 0.95×10−3µ sym. 3515 3508 3494 3478 3449 3388 1.81×10−3µ asym. 3621 3616 3600 3585 3556 3498 1.71×10−3µ apsi QM/MM 2019 40 / 53
Equations of motion µ ¨ ψi = − ∂EKS ∂ ψi| +
Λij |ψj In principle differential equations, however after discretisation difference equations (Verlet algorithm) Therefore the algorithm for the constraints Λij depends on the integration method
apsi QM/MM 2019 41 / 53
Define Xij = ∆t2 2µ Λp
ij
Yij = ∆t2 2µ Λv
ij
C wf coefficients Equations of type XX† + XB + B†X† = I − A Y = 1 2
A, B, Q of type Aij =
G c∗ GicGj
Solve iteratively: X(n+1) = 1 2
apsi QM/MM 2019 42 / 53
apsi QM/MM 2019 43 / 53
In larger molecules or crystals the structural optimisation might be difficult, especially the closer to the minimum one is CPMD can be used to perform the optimisation by simulated annealing: Rescaling the (atomic and possibly also electronic) velocities: ˙ R′
I = α ˙
RI Easy to incorporate into the velocity Verlet algorithm Optimised structure when all velocities (temperature) are ≈ 0
◮ Check by calculating the ionic forces
The ionic masses are “unphysical”: Select to “flatten” the vibrational spectrum (e. g. high mass on hydrogens) Faster convergence due to the “global” optimisation
apsi QM/MM 2019 44 / 53
µ can be chosen to be dependent on the basis set: µ (G) = µ0 , H (G, G) ≤ α (µ/α) 1
2G2 + V (G, G)
Kind of “pre-conditioning” of the equation of motion Allows for larger time step However, leads to much larger corrections on the ionic frequencies and no analytical formula can be used
apsi QM/MM 2019 45 / 53
apsi QM/MM 2019 46 / 53
Born-Oppenheimer MD Car-Parrinello MD Exactly on BO surface Always slightly off BO surface ∆t ≈ ionic time scales, ∆t ≪ ionic time scales, maximum time step possible (much) shorter time step necessary Expensive minimisation Orthogonalisation only, at each MD step less expensive per MD step Not stable against deviations Stable against deviations from BO surface from BO surface ⇒ Energy/temperature drift, thermostatting of ions necessary Same machinery in zero-gap Thermostatting of electrons systems to prevent energy exchange Many applications in solids Used in liquids, ... apsi QM/MM 2019 47 / 53
apsi QM/MM 2019 48 / 53
apsi QM/MM 2019 49 / 53
The error in the forces depends on the convergence criterion set for the electronic structure in BOMD:
apsi QM/MM 2019 50 / 53
Effect of µ: Too large value leads to loss of adiabacity Thermostatting the electrons recovers the correct behaviour
apsi QM/MM 2019 51 / 53
The radial distribution functions are correct and independent of the method used
apsi QM/MM 2019 52 / 53
Car-Parrinello method can yield very stable dynamical trajectories, provided the electrons and ions are adiabatically decoupled The method is best suited for e. g. liquids and large molecules with an electronic gap The speed of the method is comparable or faster than using Born-Oppenheimer dynamics — and still more accurate (i. e. stable)
apsi QM/MM 2019 53 / 53