Transcription
Ab initio Molecular DynamicsJürg HutterPhysical ChemistryUniversity of ZurichAugust 25, 2015
Overview Equations of motion Integrators, Sources of error General Lagrangian for AIMD Car–Parrinello MD Born–Oppenheimer MD ab initio Langevin MD Stability and efficiency
Equations of Motion (EOM)Newton’s EOM for a set of classical point particles in apotential.MI R̈ I dV (R)dR IThese EOM generate for a given number of particles N in avolume V the micro canonical ensemble (NVE ensemble).The total energy E is a constant of motion!
Total EnergyTotal Energy Kinetic energy Potential energyKinetic energy T (Ṙ) NXMII 12Ṙ2Potential energy V (R)We will use the total energy as an indicator for the numericalaccuracy of simulations.
Lagrange Equationd L L dt Ṙ RwithL(R, Ṙ) T (Ṙ) V (R)Equivalent to Newton’s EOM in Cartesian coordinates, but ismore general and flexible. Extended systems Constraints
Integration of EOMDiscretization of timeR(t) R(t τ ) R(t 2τ ) · · · R(t mτ )V (t) V (t τ ) V (t 2τ ) · · · V (t mτ ) Efficiency: minimal number of force evaluations, minimalnumber of stored quantities Stability: minimal drift in constant of motion (energy) Accuracy: minimal distance to exact trajectory
Sources of Errors Type of integratorpredictor-corrector, time-reversible, symplectic Time step τshort time accuracy measured as O(τ n ) Consistency of forces and energye.g. cutoffs leading to non-smooth energy surfaces Accuracy of forcese.g. convergence of iterative force calculations (SCF,constraints)
Velocity Verlet IntegratorR(t τ ) R(t) τ V (t) V (t τ ) V (t) τ2f (t)2Mτ[f (t) f (t τ )2M Efficiency: 1 force evaluation, 3 storage vectors Stability: time reversible Accuracy: O(τ 2 ) Simple adaptation for constraints (shake, rattle, roll) Simple adaptation for multiple time steps and thermostats
Test on Required Accuracy of ForcesClassical Force Field Calculations, 64 molecules, 330 KTIP3P (flexible), SPME (α 0.44, GMAX 25),
Stability: Accuracy of ForcesStdev. fHartree/BohrStdev. �10 1010 0810 0710 0610 0510 04DriftµHartree/ns -0.140.01-0.10-0.6315.211599.31
Born–Oppenheimer MD: The Easy WaydV (R)EOMdR IV (R) min [EKS ({Φ(r )}; R) const.] Kohn–Sham BO potentialMI R̈ I ΦForcesd minΦ EKS ({Φ(r)}; R)dR EKS const. X (EKS const.) Φi R R Φi Ri {z}fKS (R) 0
Computational Details System 64 water molecules density 1gcm 3 Temperature 330K Timestep 0.5fs DFT Calculations GPW, TZV2P basis (2560 bsf), PBE functional Cutoff 280 Rydberg, default 10 12 OT-DIIS, Preconditioner FULL SINGLE INVERSE Reference trajectory (1ps), SCF 10 10
Stability in BOMDUnbiased initial guess; Φ(t) Φ0 (R(t)) SCFMAE EKSHartreeMAE fHartree/BohrDriftKelvin/ns10 0810 0710 0610 0510 041.2 · 10 119.5 · 10 106.9 · 10 087.4 · 10 063.3 · 10 045.1 · 10 095.6 · 10 084.8 · 10 075.6 · 10 065.9 · 10 050.00.10.42.350Consistent with results from classical MDNote accuracy of forces!
Efficiency: Initial Guess of Wavefunction4th order Gear predictor (PS extrapolation in CP2K)Method SCFIterationsDrift (Kelvin/ns)Guess10 0614.380.4Gear(4)10 076.475.7Gear(4)10 065.2211.8Gear(4)10 054.6086.8What is the problem?Time reversibility has been broken!
Generalized LagrangianL(q, q̇, x, ẋ) 11M q̇2 µẋ2 E(q, y) k µG(kx yk)22y F (q, x)wavefunction optimizationG(kx yk)wavefunction retention potentialEquations of motion E E F G F kµ q y q y q E F G G Fµẍ kµ y x x y xM q̈ Extension of Niklasson Lagrangian, PRL 100 123004 (2008)
Dynamical System
Car–Parrinello Molecular Dynamicsy x G(kx yk) 0LagrangianL(q, q̇, x, ẋ) 11M q̇2 µẋ2 E(q, x)22Equations of motion E q Eµẍ xM q̈
Properties of CPMD Accuracy: MediumDistance from BO surface controlled by mass µRequires renormalization of dynamic quantities (e.g.vibrational spectra) Stability: ExcellentAll forces can be calculated to machine precision easily Efficiency: GoodEfficiency is strongly system dependent (electronic gap)Requires many nuclear gradient calculationsNot implemented in CP2K!
BOMDy Minx E(q, x) and µ 0Lagrangian1M q̇2 E(q, y)21L(x, ẋ) ẋ2 kG(kx yk)2L(q, q̇) Equations of motion M q̈ E qdecoupled equations ẍ k G x
BOMD with Incomplete SCF Convergencey F (q, x) Minx E(q, x)and µ 0Equations of motion E E F q y q G G Fẍ k x y xM q̈ SCF Error: Neglect of force terms E F y qand G F y x .EOM are coupled through terms neglected!
ASPC IntegratorIntegration of electronic DOF (x) has to be accurate: good wavefunction guess gives improvedefficiency stable: do not destroy time-reversibility of nuclear trajectoryASPC: Always Stable Predictor Corrector J. Kolafa, J. Comput Chem. 25: 335–342 (2004) ASPC(k): time-reversible to order 2k 1
Orthogonality ConstraintWavefunction extrapolation for non-orthogonal basis setsCinit KXBj 1 Ct jτ CTt jτ St jτ Ct τ{z} j 0PSCoefficients Bi are given by ASPC algorithmx Cinity[t] CtS overlap matrixJ. VandeVondele et al. Comp. Phys. Comm. 167: 103–128 (2005)
Importance of Time-ReversibilityMethod SCFIterationsDrift (Kelvin/ns)Guess10 0614.380.4ASPC(3)10 065.010.2ASPC(3)10 053.024.5Gear(4)10 076.475.7Gear(4)10 065.2211.8Gear(4)10 054.6086.8
Efficiency and DriftMethod SCFIterationsDrift (Kelvin/ns)Guess10 0614.380.4ASPC(4)10 065.010.2ASPC(4)10 053.024.5ASPC(4)10 041.621742.4ASPC(4)10 021.0321733.2
Efficiency and DriftMethod SCFIterationsDrift (Kelvin/ns)ASPC(4)10 041.621742.4ASPC(5)10 041.631094.0ASPC(6)10 041.79397.4ASPC(7)10 041.97445.8ASPC(8)10 042.0624.1
Langevin BOMDStarting point: BOMD with ASPC(k) extrapolationAnalysis of forcesfBO (R) fHF (R) fPulay (R) fnsc (R) {z}f (R) fBO : correct BO force fHF : Hellmann–Feynman force fPulay : Pulay force fnsc : non-self consistency error force f : approximate BO force
Forces in Approximate BOMDApproximate fnsc by Zf̃nsc dr Vxc (ρi ) ρi ρ VH ( ρ) I ρiwith ρ ρo ρi , ρo final (output) density, ρi initial (predicted)density.Now assumef (R) f̃nsc (R) fBO (R) γD Ṙwhere γD is a constant friction parameter.
Langevin EOMM R̈ fBO (R) (γD γL )Ṙ Θwith Θ a Gaussian random noise term andhΘ(0)Θ(t)i 6(γD γL )MkB T δ(t)Given temperature T determines γD2kB3 hM Ṙ iand an arbitrary γL this
Langevin BOMDT. D. Kühne et al., Phys. Rev. Lett. 98 066441 (2007) single SCF step plus force correction neededextremely efficient for systems with slow SCF convergence γD is small: correct statistics and dynamics difficult to stabilize in complex systems
Summary: BOMD in CP2K Born–Oppenheimer MD with ASPC(3) is default in CP2K SCF convergence criteria depends on system10 5 10 6 is a reasonable starting guess Best used together with OT and theFULL SINGLE INVERSE preconditioner Langevin dynamics is an optionneeds some special care, mostly for uniform systems
ab initio Langevin MD Stability and efficiency. Equations of Motion (EOM) Newton's EOM for a set of classical point particles in a potential. MIR I dV(R) dRI These EOM generate for a given number of particles N in a volume V the micro canonical ensemble (NVE ensemble). The total energy E is a constant of motion! Total Energy Total Energy Kinetic energy Potential energy Kinetic energy .
