Transcription
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsBasics III: Ionic Relaxation, Stress & Cell Shapes,Phonons and Molecular DynamicsDoris VogtenhuberRennes, 1. September 2016Faculty of Physics, AG-CMP, University of ViennaDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingOutline1234Ionic RelaxationIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingLattice RelaxationCell Volume OptimizationINCAR parameters in VASPPhononsIntroductionINCAR Parameters, Problem HandlingMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingIntroductionBasic ConsiderationsHermiticity of Ĥ forces on the atoms can be calculatedvia the Hellmann-Feynman theorem hΨ0 He (R) Ψ0 i hΨ0 (R) I He (R) Ψ0 (R)i I 0 (R) I RForces acting on the ions are given by the expectation value ofthe gradient of the electronic Hamiltonian in the ground-stateatomic coordinates in a cell with fixed cell shape:Hellmann-Feynman forcesgeometry of the unit cell (volume, shape):Hellmann-Feynman stressesDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingIntroduction. . . Basic Considerations V , cellshape . . . ) min.in equilibrium: E (R, minimizing E(1): find the atoms’ positions R f ( x ) with, f search for the (local) minimum of E (R)0expanded around equilibrium xf ( x ) B 11a b x x B x ā ( x x 0 )B( x x 0 )22 2fBij Hessian matrix xi xjat a stationaryPpoint the gradient of f ( gi ( x )) vanishes: f gi ( x ) x xj xj0 ) 0j Bij ( iat a minimum: B: has to be positive definiteDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingIntroduction. Basic Considerations: Newton Algorithm123start with an arbitrary point x 1calculate the gradient of f at x 1 : fx 1 x 0 )g ( x 1 ) x B( perform a step x 2 x 1 B 1 g ( x 1 )in practice: B is approximated by the largest eigenvalue of theHessian matrix, Γmax (B)steepest descent algorithmDoris Vogtenhuber
Steepest descentIonic RelaxationIntroductionLattice RelaxationAlgorithms used in VASPPhononsINCAR parameters in VASP, Problem HandlingMolecularDynamicssteepest descentHessianmatrixmate B by the largest eigenvalue of them (Jacobi algorithm for linear equations)Introduction al guess x1 ulate the gradient g x1e a step into the direction of the steepest descent 1 Γmax B g x1 x1 x2Steepest Descent Algorithmat step 2 and 3 until convergence is reachedions with long steep valleys convergence can be very slowΓ maxΓ minG. K RESSE , I ONIC1guess x 12calculate g ( x 1 )3Page 5OPTIMISATION4Doris Vogtenhuberstep along the steepestdescent direction1 x 2 x 1 Γmax g ( x 1 )repeat 2 3 converged geometry
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingIntroductionConvergence of the Steepest Descent Algorithmminimize the number of steps requested to reach the affordedaccuracy in the ion positions: step-widths along g ( x 1 )Eigenvalues of B: vibrational modes of the systemΓmax : “hardest mode” maximum stable step widthΓmin : “softest mode” slowest convergenceto reduce the error in all components to a certain fraction, thenumber of steps can be estimated from ΓΓmaxminuse preconditioning of B to speed up convergenceDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingAlgorithms used in VASPOverviewaims:12reach asymptotic convergence ratesmaintain the relaxation historyQuasi-Newton Schemes (DIIS): direct inversion in the iterativesubspaceConjugate Gradient (GC): search directions are conjugated tothe previous seach directionsDamped Molecular Dynamics (MD): minimization problem iscast into a simulated annealing approachDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingAlgorithms used in VASPThe Quasi-Newton Algorithmsimple Quasi-Newton Scheme: for a set of points x i andgradients g i (i 1, . . . , N)find a linear combination of x i which minimizes g iPconstraint: i αi 1:X g i (αi x i ) B(X B(Xαi x i x 0 )iii XXαi x 0 )iαi B( x i x 0 ) igradient: linear in its argumentsDoris Vogtenhuberαi x i Xiαi g i
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingAlgorithms used in VASPThe Full DIIS Algorithm1start with a single initial point x i2steepest descent step along gradient g ( x 1 ) : x 2 x 1 λ g 13 new gradient g 2 g ( x 2 )456search for the minimalgradient in the subspace spanned byP g i : gopt i αi g iPcalculate the corresponding position xopt i αi x i x 3 xopt λ goptDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingAlgorithms used in VASP steepest descent step from x0 to x1 (arrows correspond to gradients g0 and g1 ) gradient along indicated red line is now know, determine optimal position x 1opt g x1opt another steepest descent step form x1opt along goptfull DIIS calculate gradient x2now the gradient is known in the entire 2 dimensional space(linearity condition) and the function can be minimised exactlystart with single initial point x 0a 0x 0 a1x,1 a0 a1 1x0steepest descent step (sds)1opt. position xoptx1x21sds from xoptalong gopt x 2x0 g ( x 2 )x1optG. K RESSE , I ONICPage 13OPTIMISATIONlinearity gradient is knownin 2Dminimize f exactlyDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingAlgorithms used in VASPThe Conjugate-Gradient Algorithm (CG)search directions: conjugated to the previous seach directionsstart with x 0123steepest descent step along gradient with line minimizationgradient at the current position g ( x N )conjugate g ( x N ) to the previous search direction: s ( x N ) g ( x N ) γ g ( x N 1 ),45γ ( g ( x N ) g ( x N 1 )) · g ( x N ) g ( x N 1 ) · g ( x N 1 )line minimization along s Nif g is not sufficiently small: continue with 1search directions are orthogonal (step 3): s N B s M N, MCG finds the min. of a quadratic function with k DOF ink 1 steps exactlyDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingAlgorithms used in VASP teepest descent step from x0 , search for minimum along g0 by performing several trialteps (crosses, at least one triastep is required) x1 etermine new gradient g1 g x1 and conjugate it to get s1 (green arrow)or 2d-functions the gradient points now directly to the minimumminimisation along search direction s1 The CG Algorithmsds from x 0 along g 0x0trial step(s) x,Nx 1), x 1x1 new g 1 g ( x 1 )x2x0conjugate g 1 : , s 1x1 s 1 points directly towardsthe minimums1x1G. K RESSE , I ONICOPTIMISATIONDoris VogtenhuberPage 15minimization along s 1
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingAlgorithms used in VASPDamped MD (MD)atoms’ positions x are regarded as dynamic degrees of freedomforces ( gradients) accelerate the motion of the atoms µ x equation of motions of the atoms: x 2αFintroduce an additional friction term µintegration of this eqation: simple velocity Verlet algorithm vN 1/2 xN 1 N /(1 µ/2)(1 µ/2) vN 1/2 2αF xN vN 1/2Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingAlgorithms used in VASPaves like a rolling ball with a frictionll accelerate initially, and then deaccelerate when close to the minimume optimal friction is chosen the ball will glide right away into the minimumtoo small friction it will overshoot the minimum and accelerate backtool large friction relaxation will also slow down (behaves like a steepestDampedent)x0MD (MD)“rolling ball” with friction(µ)µ too small: minimumovershot,back-accelerationG. K RESSE , I ONICOPTIMISATIONDoris Vogtenhuberµ too large: relaxationslows downPage 18
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingINCAR Parameters in VASPOverviewAlgorithmDIISCGdamped MDmain flagIBRION 1IBRION 2IBRION 3additional flagsPOTIM, NFREEPOTIMPOTIM, SMASSterminationEDIFFGEDIFFGEDIFFGEDIFFG “convergence criterium”:EDIFFG 0 :EDIFFG 0 : (E N E N 1 ) EDIFFG N EDIFFG i 1, Nions FiDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingINCAR Parameters in VASPParameter Usage by the Algorithms of VASPDIISPOTIM ( 0.5) generally determines the step size (no lineminimizations)NFREE # of ionic steps stored in the iteration history: for theset of points x i and gradients g i (i 1, . . . , N)NFREE ( 5) max(N)CGPOTIM ( 0.5) : size of the first trial step, the subsequent lineminimization is performed using Brent’s algorithm damped MD: in vN 1/2 (1 µ/2) vN 1/2 2αF N /(1 µ/2)POTIM α, good choices: 0.15 POTIM p 0.4SMASS ( 0.4) µ, which should be 2 Γmin /ΓmaxDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingChoice of the WhymostAppropriateAlgorithmso manyalgorithms :-(.decision chartyesCGReally, this is too complicatednoyesyesclose to minimum1 3 degrees of freedomnonoyesdamped MD or QUICKMINvery broad vib. spectrum 20 degrees of freedomnoDIISG. K RESSE , I ONIC OPTIMISATIONDoris VogtenhuberPage 25
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingProblem Handling(Some Other) Reasons for Bad Convergenceunreasonable starting geometry (POSCAR)lattice parameters, atomic positionscheck OUTCAR for interatomic distances, forces of the inputgeometry, external pressure, ( if ISIF 0)sub-optimal settings of (some) INCAR parametersbad electronic convergence of (one of) the ionic steps wrong forcescheck OSZICAR for the convergence of each ionic step: dE ,charge density convergenceincrease NELM, decrease the (spin density) mixing parameterschoose a different BZ-integration method ISMEAR, SIGMAchoose a different electronic relaxation algorithm ALGObasis sets too small ( aliasing errors)is the k-mesh appropriate? (modify KPOINTS)Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingProblem HandlingAliasing Errorsrelated to errors caused by the truncated FFT gridRfolding theorem: ρ ψn ψn (VGH ) contain components up to2n 2Gcutoff after back-transformation from ρr to ρG (V )residual Vector (V ψ): components up to 3GcutoffFourier grid has to include all wave-vectors up to 2Gcutoff .if this is not the case: aliasing (“wrap around”) errors:components of ρ are wrapped around from the other side ofthe box due to the periodicityhigh frequency components are aliased to low-frequencycomponentsDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingProblem HandlingDrifts in the Forcesimpact of aliasing errors on the results::in a lattice with perfect translation symmetry: if all atoms ofthe cell are shifted by the same translation vector,E has to remain exactlysamePNatthe i 0forces sum up to 0: i 1Faliasing errors destroy the translational invariance: atoms equivalent by symmetry are equivalent no longer drifts explicitely unless ISYM 0BUT VASP symmetrizes ρ and FDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingProblem HandlingDrifts in the Forcesto reduce drifts in the forces (written in OUTCAR)bulk & surfaces: increase the precision ENCUT, PRECsurfaces: in 3D periodic cell, the origin of the cell is arbitrary,i.e. the slab may start drifting through the vacuumkeep (at least) one layer fixed (Selective Dynamics optionin POSCAR)polar surfaces: include dipole corrections (IDIPOL, LDIPOL)to avoid artificial electrostatic forces across the vacuum layerDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsCell Volume OptimizationINCAR parameters in VASPOutline1234Ionic RelaxationIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingLattice RelaxationCell Volume OptimizationINCAR parameters in VASPPhononsIntroductionINCAR Parameters, Problem HandlingMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsCell Volume OptimizationINCAR parameters in VASPCell Volume OptimizationIntroductionthe equilibrium volume Veq and shape of a crystal calculatedfrom ab initio depend on the XC-type used:LDA: overbinding a0 too smallPBE, PW91: underbinding a0 too largeresults are improved using specially designed functionals(PBEsol, HSE),. accurate calculations should always be performed for thecell at equilibrium for the respective XC-type to avoid artifacts(unless there is good reason not to do so)Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsCell Volume OptimizationINCAR parameters in VASPCell Volume OptimizationStrategies to obtain Veq“by hand”: series of calculations at different cell volumes, Veq V (min(E (V )):very old-fashioned, almost impracticable for non-cubic cells“by hand”: fitting to thermodynamic equations:eg. Birch-Murnaghan fitVASP: automatic optimization, based on the calculatedHellmann-Feynman stressesthe automatic geometry optimization sensitively depends onthe quality of the used basis sets:E-cutoffs (completeness of the basis set), FFT-grids k-meshesDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsCell Volume OptimizationINCAR parameters in VASPCell Volume OptimizationEnergy Cutoff: Basis Setsat each k, the plane waves that are included in the basis have 2 to fulfill the criterium 2m G k 2 EcutoffeEcutoff defined by ENCUT: default: max(ENMAX), given inPOTCAR for each element 2 changes of cell volume and -shapeEcutoff G the default cutoff should only be used for calculations withfixed cell-shape and -volume, eg.frozen phononssurface and adsorption calculationsMD (NVT ensemble)Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsCell Volume OptimizationINCAR parameters in VASPCell Volume OptimizationFixed basis-set calculationsGcutthe cutoff decreases by a factor τ1 τ1when the lattice is expanded fromexplanation:τ1 τ1 τ2 τ10b2b 2π/τ11 latticeexpandedτ1set1 for theexpandedlatticeτthebasiscutoff decreasesbya alowera factorcorrespondseffectively tocut- ττ101off Gcut and therefore a lower0effective cutoff Gcutisquality,lowerthe energy is overestimated atE is overestimated for larger V slarger volumesthevolumeapparentVeq is too smalltheis underestimatedforfixed basis-set calculationsG’cut b2τ1b 2π/τ11Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsCell Volume OptimizationINCAR parameters in VASPnsteadfixed cutoffCell ofVolumeOptimization-3.4Improvement using fixed basis sets?start from WAVECAR with ISTART 2E (eV)240 eV, basis set fixed240 eV, cutoff fixed270 eV, basis set fixed270 eV, cutoff fixedNO!!, becauseeffective decreases with increasing VEcutoff quality of the basis set becomesworse with increasing V-3.6 min(E (V )) is shifteddense k meshes necessary to obtain 2 )!smooth curves ( k G1111.51212.5133volume V (A )Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsCell Volume OptimizationINCAR parameters in VASPCell Shape relaxationsess tensor%,pressure (kBar)onStress TensorVASP does not adopt the basis set ina run0u-stress tensor σij : implicitely calculatedwith a fixed-basis-set setup-500for Cu (270eV): contraction predictedby error (p -50 kB)increase ENCUT (by 30%) todefault cutoffceRACY AND-1000200250300350400cutoff energy E(eV)VALIDAIONOF RESULTSPage 9Doris Vogtenhuberperform lattice relaxationscalculate stress tensors and pressure(P 13 Trσij )
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsCell Volume OptimizationINCAR parameters in VASPCell Shape relaxations“recipe” for Determining Cell Shapesalways use an increased cutoff: ENCUT 1.3*max(ENMAX)do it step-wise:1234start with 1-2 steps (NSW) from your guessed input geometry(coarse pre-relaxation)delete WAVECARcontinue from CONTCAR with slightly more stepsrepeat 1-3 until the remaining pressure (and stress tensorcomponents) are in accordance with the afforded accuracyif the space-group of the system is known, use ISYM 2 toavoid symmetry violations due to numerical errorsDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionINCAR Parameters, Problem HandlingOutline1234Ionic RelaxationIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingLattice RelaxationCell Volume OptimizationINCAR parameters in VASPPhononsIntroductionINCAR Parameters, Problem HandlingMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionINCAR Parameters, Problem HandlingIntroductionBasicsvibrations of the crystal lattice influenceelasticthermodynamicopticalelectronic transport properties“soft modes” indicatephase transitions (bulk) ordissociation (dissociative adsorbtion processes)Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionINCAR Parameters, Problem HandlingIntroductionBasics 0 (lm) is displaced by uif atom m in cell l of a crystal R R(lm) R0 (lm) u (lm)Pkinetic energy: T 12 lmα Mm u̇α2 (lm), potential energy: expanded V (R(lm)) 0 (lm)) V0 (R {z } V0 0 12X V (R(lm)) lmα Xlmα,l 0 m0 β Rα (lm){zuα (lm) 0 in equilibrium 2 V (R(lm)) Rα (lm) Rβ (l 0 m0 ) {z}Φαβ (ll 0 mm0 ) force constant 0 (lm)Φαβ (ll 0 mm0 ): derivative taken at R(lm) RDoris Vogtenhuber}uα (lm)uβ (l 0 m0 )
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionINCAR Parameters, Problem HandlingIntroductionBasicsΦαβ (ll 0 mm0 ): component α of the force acting on atom (lm),caused by the displacement of atom (l 0 m0 ) in direction βequations of motionMm üα (lm) X V Φαβ (ll 0 mm0 )uβ (l 0 m0 ) uα (lm)0 0l m βuse symmetryharmonic ansatz: uα (lm, t) ω 2 eα (m) X Mm eα (m)e i qRl e iωtX1 eβ (m0 ) ( (Mm Mm0 ) 2 Φαβ (ll 0 mm0 )e i q(Rl 0 Rl ) )β,m0l0 Doris Vogtenhuber{zDαβ (mm0 , q ) dynamical matrix}
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionINCAR Parameters, Problem HandlingPhononsBulkVASP calculates phonons at the zone-center supercell approachthe elements of the Hessian Matrix are calculated either byfinite displacement of the ions:IBRION 5,6; NFREE, POTIMassume: displacements are within the harmonic limitusing density functional perturbation theory IBRION 7,8the tags making use of the symmetry (IBRION 6,8) can beused in vasp.5.2 onlyDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionINCAR Parameters, Problem HandlingPhononsVibrational Modes of Moleculesselect vibrational modes of interest using SelectiveDynamics in POSCAR (eg. change of the modes of amolecules upon adsorption on a surface)vibrational frequencies of adsorbates: usually calculatedaccurately ifonly the adsorbate itself and the NN substrate atoms are notkept fixedin any case: test how many “shells” have to be included toconverge the frequenciessaves computing timeDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionINCAR Parameters, Problem HandlingProblem Handlingpossible sources of errorsnegative frequencies (imaginary modes):may indicate structural instabilities (mode softening),or: calculation not properly converged increase EDIFFfrom the default valueVASP.4.6 only: POTIM has to be set explicitely:recommended POTIM 0.015 or smaller, the default value(0.5) certainly is not within the harmonic limit. unreasonable frequenciesVASP can not continue from an unfinished run. for thecalculation of eg vibration frequencies of adsorbates (largenumber of atoms in the unit cell): reduce the calculatedvibration modes to a reasonable numberDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPOutline1234Ionic RelaxationIntroductionAlgorithms used in VASPINCAR parameters in VASP, Problem HandlingLattice RelaxationCell Volume OptimizationINCAR parameters in VASPPhononsIntroductionINCAR Parameters, Problem HandlingMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPIntroductionGeneral Remarksclassical equations of motion (EOM) for atoms in amicrocanocical NVE ensemble (p: momenta , q: positions)H(p, q) NXp i 2i 1dp H(p, q) dt q,mi V (q1 , .qn )dq H(p, q) dt pergodic hypothesis: ensembe and time averages are related:RhAiH HdpdqA(q)e kB T1 R k HTτdpdqe BZτdtA(t)0 MD can be used to compute observables A.Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPStandard Versionstandard MD: on the Born-Oppenheimer surface,Hellmann-Feynman forces, Thermostat: Nosé i (t) Newtonian EOM for the set of atoms i, Mi R E i (t) R coupled set of equations, wavefunctions kept areorthonormal via a Lagrangian multiplier λijµψ̈i ( r , t) δEδψi ( r , t) Xλij ψj ( r , t)jVerlet algorithm with damping factor (friction term) µ vN 1/2 xN 1 N )/(1 µ/2)((1 µ/2) vN 1/2 2αF xN vN 1/2 Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPIntroduction: (Chemical) Reactionsfor any reaction, according to Arrhenius’ law:‡dc0 (i) E kc0 (i) k Ae kB T ,dtA : Arrhenius prefactor A‡Eyring-Polanyi theory: k kBhT e kB T A‡ . . . free energy difference between the transition state (‡)and the initial state (0).the free energy A can be evaluated via statisticalthermodynamics:Ai kB T log Qi k Doris VogtenhuberkB T Q ‡· 0hQQ : partition function
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPIntroductionSome Statistics’ BasicsQ tot for species i Qitot Qitrans Qirot Qivib Qielecronicin extended systems with translational symmetry:Qitrans and Qirot are constant and cancel outQivib : Qivib Qielectronic :QMQielhνi2kB Thνi 1 e kB Tei 1 e k EiBT (harmonic approx.)Qiel,‡Qiel,0 e‡ EiTB kthe reaction constant k is given as‡k ikB T Q vib,‡ E· vib,0 · e kB ThQDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniquesstandard version of MD:uses Carthesian coordinatestransition states: obtained using the Nudged Elastic Band(NEB) method inefficient, slow for chemical reactionsimprovement: Advanced MD Techniquesinstead of cartesian coordinates: use a more clever choice ofdelocalized, internal coordinates ξ (bond lenghts, -angles,. . . )ergodic hypothesis used to calculate hAi via its time averageimplemented in VASP.5 by Tomas Buckocompile VASP with -Dtbdyn to replace standard MD byadvanced MD techniquesDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniquesin systems with richly structured Potential EnergyHypersurfaces (PES): forces on the atoms might not drag thesystem over an energy barrier of the PES the system gets stuck in a basin of the PESmethods to avoid this behavior:add a bias potential Ṽ (ξ) to enhance the sampling in regionsof the PES with low probability P(ξi ) (eg transition stateregions):“umbrella sampling”constrain the MD by adding geometrical constraints viaadditional terms in the Lagrangian, enforcing the constraint“blue moon sampling”Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniques: Biased MDa bias potential Ṽ (ξ) is used to enhance sampling of theinternal coordinate ξ(q)H̃(p, q) P̃(ξi ) H(p, q) Ṽ (ξ),ξ ξ(q)R H̃δ(ξ(q) ξi )e kB T dpdqhδ(ξ(q) ξi )iH̃ R H̃e kB T dpdq2. recover the correct distribution of A at the end by usinghAiH ṼkB TA(q)eṼe kB TH̃H̃Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniques: Metadynamicsadditional DOFs (α) driving the reaction: ξα , ξ α (velocity)and mass µα , are coupled to the relevant geometricalparameters (collective variables Ξα (x)) via harmonic springswith force constants kα :L L0 X1X12µα ξ α kα (Ξα (x) ξα )2 Ṽ (t, ξ)22alphaalpha ξ(t) ξ(itG ) 2Pt/t2w 2Ṽ (t, ξ) h i 1G e sum of Gaussian hills (hi , wi ) updated at every time-step tGduring the calculationtG : 1-2 orders of magnitude than t of the MDA(ξ)t limt Ṽ (t, ξ) constDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniques: Constrained MDmodify the Lagrange multiplier L by adding a term includingall geometric constraints r : L(q, q̇) L(q, q̇) rXλ i σii 1with σi ξi (q) ξi , ξi . . . fixed variable1 standard leap-frog MD to obtain qi (t δt)2 use new positions to compute λi constraints3 update v and q by adding a contribution due to the restoringforce ( λ) qi (t δt)4 repeat 1-3 until σ(q) matches the convergence criterium(SHAKE algorithm)Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniques: TD integration of A-gradientsT. Bucko, J.Phys.Cond.Matt. 20, 064211 (2008)the #DOF dynamical variables of Ĥ are split intothe active reaction variable ξ (pξ , qξ ), defining the reactionpath 1 2, (slow modes)inactive set q {q1 , . . . , qM 1 }, pq (fast modes; not frozen,but do not contribute to the minimum A-path as their thermalmotions are nearly harmonic)Z ξ(2) A A1 2 dξ ξ ξ ξ(1)ξ is constrained to remain constant to ξ ξ 0,also, pξ is not sampled in the MDDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniques: TD integration of A-gradientsthe Hamiltonians of the constrained (Hξc ) and unconstrained (H)ensembles are:Hξc 1 tp Xp V (q, ξ)2H 1Hξc pξt (Y · pq ) (pξt Z pξ )2withXα,β 2MMMXXX1 qα qβ1 ξ qβ1 ξ, Yα ,Z mi xi ximi xi ximi xii 1i 1(α, β 1, . . . M 1)Doris Vogtenhuberi 1
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniques: TD integration of A-gradientsconstrained and unconstrained ensemble averages of aquantity O are related via a “blue moon” correctionE.A.Carter et.al., Chem.Phys.Lett 156, 472 (1989)DE1OZ 2ξ hOi D 1 E Z 2 ξthe constraints on the system to remain on the reaction pathare included via the Lagrangian multiplier λ (accounting forthe reaction coordinate, calculated using the SHAKEalgorithm) in the modified LagragianL (x, ξ, ẋ) L(x, ẋ) λ(ξx ξ)Doris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniques: TD integration of A-gradientsthe free energy gradients can then be calculated: A ξ ξ 11 Z 2 ξ *"Z 21 λξ kB TZ 1MX1 ξ Zmi xi xii 1# ξ crucial for blue moon ensemble techniques: appropriate choiceof the parameter ξDoris Vogtenhuber
Ionic RelaxationLattice RelaxationPhononsMolecular DynamicsIntroductionMD Algorithms implemented in VASPThermostats implemented in VASPMD Algorithms implemented in VASPAdvanced MD Techniques: Sl
Damped Molecular Dynamics (MD): minimization problem is cast into a simulated annealing approach Doris Vogtenhuber. Ionic Relaxation Lattice Relaxation Phonons Molecular Dynamics Introduction Algorithms used in VASP INCAR parameters in VASP, Problem Handling Algorithms used in VASP The Quasi-Newton Algorithm simple Quasi-Newton Scheme: for a set of points xi and gradients gi (i 1;:::;N) nd .
