Tutorial Series: Ab Initio Molecular Dynamics - FHI Conference System .

Transcription

Tutorial Series:Ab initio molecular dynamicsMariana RossiMax Planck Institute for the Structure and Dynamics of MatterFritz Haber Institute of the MPG

Addressing nuclear degrees of freedomSchrödinger Equation̂ EΨHΨ Ψ(t)̂iℏ HΨ(t) t2

Addressing nuclear degrees of freedomSchrödinger Equation̂ EΨHΨ Ψ(t)̂iℏ HΨ(t) tBorn-Oppenheimer approximation: Ψ ψe χn2

Addressing nuclear degrees of freedomSchrödinger Equation̂ EΨHΨ Ψ(t)̂iℏ HΨ(t) tBorn-Oppenheimer approximation: Ψ ψe χn2

Addressing nuclear degrees of freedomSchrödinger Equation̂ EΨHΨ Ψ(t)̂iℏ HΨ(t) tBorn-Oppenheimer approximation: Ψ ψe χn2

Addressing nuclear degrees of freedomSchrödinger Equation̂ EΨHΨ Ψ(t)̂iℏ HΨ(t) tBorn-Oppenheimer approximation: Ψ ψe χn2

Nuclear Hamiltonian (classical) and the canonical ensembleH(p, q) N2pI 2MII 1 VBO(q1, , qN ) Probability in the canonical ensemble (N, V, T fixed) given by𝒫(p, q) e βH(p,q)Zβ 1/kBTZ dp dqe βH(p,q)3

Probabilities and averages In classical mechanics, positions and momenta are not correlated𝒫(p, q) 𝒫(p)𝒫(q) e2 βp /2M2/2M βp dpee βV(q) dqe βV(q)4

Probabilities and averages In classical mechanics, positions and momenta are not correlated𝒫(p, q) 𝒫(p)𝒫(q) e2 βp /2M2/2M βp dpee βV(q) dqe βV(q) The blue term is easy to calculate, while the red term tends to be very difficult.Important because:⟨A(q)⟩ dpdq𝒫(p, q)A(q) dqA(q)𝒫(q) ⟨A(q0)B(qt)⟩ dpdq𝒫(p, q)A(q0)B(qt) dq𝒫(q)A(q0)B(qt) 4

Static and equilibrium properties at finite T⟨A(q)⟩ dpdq𝒫(p, q)A(q) dqA(q)𝒫(q) q2Energyq15

Static and equilibrium properties at finite T⟨A(q)⟩ dpdq𝒫(p, q)A(q) dqA(q)𝒫(q) q2Energyq15

Static and equilibrium properties at finite T⟨A(q)⟩ dpdq𝒫(p, q)A(q) dqA(q)𝒫(q) Calculating such integrals on a grid scales as G(G grid points, N atoms)Instead, it is better to perform some sort of importance samplingq2q1Energy 3N5

Static and equilibrium properties at finite T⟨A(q)⟩ dpdq𝒫(p, q)A(q) dqA(q)𝒫(q) Calculating such integrals on a grid scales as G(G grid points, N atoms)Instead, it is better to perform some sort of importance samplingGenerate manyconfigurations qsaccording to 𝒫(q) andcalculateq2q1M1⟨A(q)⟩ A(qs)M s 1Energy 3N5

The Hamiltonian time evolutionpIH(p, q) VBO(q1, , qN ) 2MIp HI·q I pI MI H·p I V(q)I qI6

The Hamiltonian time evolutionpIH(p, q) VBO(q1, , qN ) 2MIp HI·q I pI MI H·p I V(q)I qIConfigurations created by this time evolution are consistent with and conserve 𝒫(p, q)d𝒫(p, q)dH βH e 0dtdt𝒫(p0, q0) 𝒫(pt, qt)6

Static and equilibrium properties at finite TErgodic hypothesis: ensemble average equal to time averageτ1⟨A(q)⟩ dpdq𝒫(p, q)A(q) limτ A[q(t)]dt τ 0q2q1Energy 7

Static and equilibrium properties at finite TErgodic hypothesis: ensemble average equal to time averageτ1⟨A(q)⟩ dpdq𝒫(p, q)A(q) limτ A[q(t)]dt τ 0q2q1Energy 3.%7

Static and equilibrium properties at finite TErgodic hypothesis: ensemble average equal to time averageτ1⟨A(q)⟩ dpdq𝒫(p, q)A(q) limτ A[q(t)]dt τ 0Can be calculated through “molecular dynamics”q2q1Energy 7

How to get (classical) trajectories?Molecular dynamics(q1, p1)(q2, p2)(q3, p3)1. Assign initial q (position) and p (momenta)8

How to get (classical) trajectories?Molecular dynamics(q1, p1)(q2, p2)(q3, p3)1. Assign initial q (position) and p (momenta)2. Evolve (numerically) Newton’s equation of motionfor a finite time increment8

How to get (classical) trajectories?Molecular dynamics(q1, p1)(q2, p2)(q3, p3)1. Assign initial q (position) and p (momenta)2. Evolve (numerically) Newton’s equation of motionfor a finite time increment(q1(t dt), p1(t dt))(q2(t dt),p2(t dt))3. Assign new position and momenta(q3(t dt), p3(t dt))8

How to get (classical) trajectories?Molecular dynamics(q1, p1)(q2, p2)(q3, p3)1. Assign initial q (position) and p (momenta)2. Evolve (numerically) Newton’s equation of motionfor a finite time increment(q1(t dt), p1(t dt))(q2(t dt),p2(t dt))3. Assign new position and momenta(q3(t dt), p3(t dt))8

Integrating the equations of motion Taylor expansion of q123···q(t Δt) q(t) Δtq(t) q(t)Δt 𝒪(Δt )29

Integrating the equations of motion Taylor expansion of q 123···q(t Δt) q(t) Δtq(t) q(t)Δt 𝒪(Δt )2123···q(t Δt) q(t) Δtq(t) q(t)Δt 𝒪(Δt )224··q(t Δt) q(t Δt) 2q(t) q(t)Δt 𝒪(Δt ) V(q(t)) F(t) MM9

Integrating the equations of motion Taylor expansion of q 123···q(t Δt) q(t) Δtq(t) q(t)Δt 𝒪(Δt )2123···q(t Δt) q(t) Δtq(t) q(t)Δt 𝒪(Δt )224··q(t Δt) q(t Δt) 2q(t) q(t)Δt 𝒪(Δt ) V(q(t)) F(t) MM"Velocity Verlet"F(t)·q(t Δt) q(t)· ΔtM“Verlet” is time reversible and symplectic9

Symmetric splitting of velocity Verlet Implemented in most codes (including FHI-aims)One StepΔtF(t)Δt·q(t ) q(t)· 2M 2Δt·q(t Δt) q(t) q(t )Δt2Evaluate F(t Δt) from q(t Δt)ΔtF(t Δt)Δt·q(t Δt) q(t· ) 2M2Compute any instantaneous properties at step Δt10

The ab initio forces“Evaluate F(t Δt) from q(t Δt)”F(q) q [minρE[ρ]] Calculating these forces tends not to be trivial within an electronic structure architecture11

The ab initio forces“Evaluate F(t Δt) from q(t Δt)”F(q) q [minρE[ρ]] Calculating these forces tends not to be trivial within an electronic structure architectureIn FHI-aims we have (at least!):F FHFP F FNSC FmpHartreepot.PulayNon ncl. GGA, mGGA,Hartree pot.and relativity)11

The ab initio forces“Evaluate F(t Δt) from q(t Δt)”F(q) q [minρE[ρ]] Calculating these forces tends not to be trivial within an electronic structure architectureIn FHI-aims we have (at least!):F FHFP F FNSC FmpHartreepot.PulayNon ncl. GGA, mGGA,Hartree pot.and relativity) FvdW Fexx 11

The Hellman-Feynman and Pulay forces The Harris functional1E[ρ] fiεi ρ(r)νH dr ρ(r)νxcdr Exc[ρ] ENN 2i12

The Hellman-Feynman and Pulay forces The Harris functional1E[ρ] fiεi ρ(r)νH dr ρ(r)νxcdr Exc[ρ] ENN 2iZI ZJ1 RI RJ 2 IJ12

The Hellman-Feynman and Pulay forces The Harris functional1E[ρ] fiεi ρ(r)νH dr ρ(r)νxcdr Exc[ρ] ENN 2i The forces Hellman-Feynman forcesZI ZJ1 RI RJ 2 IJδE[ρ(r; R), 0 -9R]F δRϵi ⟨ψi ĥ KS ψi⟩- E νNNextHFF ρKS(r)dr RI RI12

The Hellman-Feynman and Pulay forces Pulay forces come from i fiδεi ifi⟨ψi δĥ KS ψi⟩ ifi⟨δψi ĥ KS ψi⟩ fi⟨ψi ĥ KS δψi⟩iThey would be zero in the complete basis set limit even for atom-centred basis sets13

The Hellman-Feynman and Pulay forces Pulay forces come from ifiδεi fi⟨ψi δĥ KS ψi⟩ i fi⟨δψi ĥ KS ψi⟩ i fi⟨ψi ĥ KS δψi⟩i They would be zero in the complete basis set limit even for atom-centred basis sets By using ψi cμiφμ and ciμhμν ϵi ciμsμν, plus imposing δ⟨ψi ψi⟩ δ( sμνc*iμciν) 0μμPFc 2Reμ iμν μνfic*c⟨iνiμ φμ R ĥ KS ϵi φν⟩“Beyond LDA” functionals and relativistic corrections add further terms to the Pulayforces. See Blum et al., CPC (2009)13

Hartree potential non-self-consistent force correction Still not the whole story. The energy functional at a given m SCF step is actually1(m) (m)(m 1)(m 1)(m)(m 1)(m 1)(m 1)E [ρ] fi εi ρ(r)νH dr ρ(r)νxc dr Exc[ρ] ENN 2iMixed densityat step m-1 After some math we getFNSC (ρ (m 1)(r) (m)ρKS (r))(m 1)δνHUnmixed densityat step mδRIdr (ρ (m 1)(r) (m)ρKS (r))(m 1)δνxcδRIdr14

Hartree potential non-self-consistent force correction Still not the whole story. The energy functional at a given m SCF step is actually1(m) (m)(m 1)(m 1)(m)(m 1)(m 1)(m 1)E [ρ] fi εi ρ(r)νH dr ρ(r)νxc dr Exc[ρ] ENN 2iMixed densityat step m-1 After some math we getFNSC (ρ (m 1)(r) (m)ρKS (r))(m 1)δνHUnmixed densityat step mδRIdr (ρ (m 1)(r) (m)ρKS (r))(m 1)δνxcδRIdrImplemented in 2021 by Herzain Rivera & Mariana Rossi14

Need of accurate self-consistency Hamiltonian dynamics (microcanonical ensemble): Energy should be conserved!How close to self consistency do we need to be?15

Need of accurate self-consistency Hamiltonian dynamics (microcanonical ensemble): Energy should be conserved!How close to self consistency do we need to be?large thresholds forenergy and density convergence(cheaper!)BOMD: C2H6small thresholds forenergy and density convergence15

Hartree potential non-self-consistent force correctionImpact of non-self-consistent force correction term from Hartree potential0.50.4Energy (eV) Zundel cation— 1e-6, no NSC— 1e-3, with NSC— 1e-2, with NSC0.30.20.10.00246Time (ps)81016

Time step and the accuracy of integrationΔt 3fs17

Time step and the accuracy of integrationΔt 3fs E.17

Time step and the accuracy of integrationn en0Δt 3fsĤG G0kb T(1) E(2)X p2IH(R, p) V (R)2MI(3)rI V (R) ! MI R̈I FI(4)IṗI "@H @RItime step? What is a goodṘI pI /MI Depends on the highest(5)vibrational frequency (thus mass) ofEyour VBO (R)( r;R)(6)systemDF Tp( ! k/M )(7)1 X2L MI ṘI µ2IZXichoose a time step Typically,# Zmax)corresponding to 1/(10ω2 dr i(femtosecond(r, t) E[n(r;timeR)] 2 ijdr i (r, t) j (r, t)scale)Z 17

ṘI pI /MIA bit about Car-Parrinello MD (5)R. Car and M.Parrinello, Phys. Rev. Lett. 55, 2471-2474 (1985) VBO (R) EDF T ( r; R)(6)p(fictitious) degreesof freedom for the electronsin the Lagrangian(7)! k/MExtended Lagrangian: addsolve coupled equations of motion"1 X2L MI ṘI µ2IZX2 dr i (r, t) i#V( , ; R) I2 ijif Zµ i i (r, t) j (r, t)drZdr i (r, t) j (r, t)(8)ij andij(9) MI R̈I rI V ( , ; R)(10) 1 V ( , ; R) X j ji (11) 2ijZ 1 X A eZEii B(t t)A(t ) 1 XeEiB (t t)A (t ) 1Z1Ai TTXZeEi(12)iTdtA(t)(13)0dt0 B(t t0 )A(t0 )(14)18

ṘI pI /MIA bit about Car-Parrinello MD (5)R. Car and M.Parrinello, Phys. Rev. Lett. 55, 2471-2474 (1985) VBO (R) EDF T ( r; R)(6)p(fictitious) degreesof freedom for the electronsin the Lagrangian(7)! k/MExtended Lagrangian: addsolve coupled equations of motion"1 X2L MI ṘI µ2IZX2 dr i (r, t) i#andLagrange multipliersV( , ; R) 2ijKohn-ShamorbitalsFictitiouselectron mass Zµ i i (r, t) j (r, t)drZdr i (r, t) j (r, t)(8)ij ij(9) MI R̈I rI V ( , ; R)(10) 1 V ( , ; R) X j ji (11) 2ijZ 1 X A eZEii B(t t)A(t ) 1 XeEiB (t t)A (t ) 1Z1Ai TTXZeEi(12)iTdtA(t)(13)0dt0 B(t t0 )A(t0 )(14)18

MI R I F I(4)ṘI pI /MI(5) R) MṘI pI /MIA bit about Car-Parrinello MDVBO (R) EDF T ( r; R)(6)p R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471-2474 (1985) k/M VBO(R)EDF T ( r; R)(7) (6)! Extended Lagrangian:(6) add (fictitious)solve coupled equationsof motion(7)"ZX#pdegreesof! k/M#1 XZ2 "X2 i (r, t) XL MI Ṙ1I µdr V(,22 2L ZMI ṘI µdr i (r, t) V( ,(5) I; R) 22drijZµ i iiI t)iji (r, t) j (r,Fictitious-(8)electron mass dr i (r, t) j (r, t) ij (9) MI R̈I rI V ( , ; R)(10)1 V ( , ; R) X (11)jji 2I freedom for the electronsin the Lagrangian and(7); R) 2V( ,ij Z; R) 2e(12)ij(8)jXXZ EeiiEi Z T 1Ai eZiTEi(12)(12)1 XX A edtA(t)(13)11iEiZT0 A eA dtA(t)(13)Z Ti iZT1ZEiT 0XiAi dtA(t)(13) 11Ei TB(t t)A(t ) eB (t t)A (t Z) dt0 B(t t0 )A(t0 )(14)Z 1 XeZi (r, t) j (r, t)-Z Eidrij(8)ijZZ Kohn-Sham dr i (r, t) j (r, t) ij (9)dr i (r, t) j (r, t) ij (9)orbitals MI R̈Mr V ( , ; R)(10)I ; R)(10)I R̈I I rI V ( , X X1V(,;R)1V(,;R)µ i µ i j jij (11)ji (11) 2 i2ijjXLagrange multipliers Z dr i (r, t) j (r, t)18

MI R I F I(4)ṘI pI /MI(5) R) MṘI pI /MIA bit about Car-Parrinello MDVBO (R) EDF T ( r; R)(6)p R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471-2474 (1985) k/M VBO(R)EDF T ( r; R)(7) (6)! Extended Lagrangian:(6) add (fictitious)solve coupled equationsof motion(7)"#ZXpdegreesof! k/M#1 XZ2 "X2 i (r, t) XL MI Ṙ1I µdr V(,22 2L ZMI ṘI µdr i (r, t) V( , I; R) 2ij2drµ i iiI t)iji (r, t) j (r,Fictitious(8)electron mass dr i (r, t) j (r, t) ij (9)Z MI R̈I rI V ( , ; R)(10)1 V ( , ; R) X (11)jji 2I freedom for the electronsin the Lagrangian and(7); R) 2V( ,ij Z; R) 2Lagrange multipliers Z dr i (r, t) j (r, t)ijdri (r, t) j (r, t)EieEi1Ai dtA(t)(13) 1 X TB(t t)A(t ) eij(8)jXXEiZ (12)Eeitobe verysmallfictitious mass of the electrons needZ e(12) Adiabatic separation:XZZ e(12)1 X1Z 1AX eA at lfconsistencycalculationstep1 XZT A eA dtA(t)(13)ZiT(8)ijZZ Kohn-Sham dr i (r, t) j (r, t) ij (9)dr i (r, t) j (r, t) ij (9)orbitals MI R̈Mr V ( , ; R)(10)I ; R)(10)I R̈I I rI V ( , X X1V(,;R)1V(,;R)µ i µ i j jij (11)ji (11) 2 i2ijj1Z(5)EiEiZiiismaller time stepiTiT0TT 0Zi1EiB (t t)A (t Z) dt0 B(t t0 )A(t0 )(14)18

Other Acceleration TechniquesJ. Kolafa JCC 25, 335 (2004);T. Kühne, et al. PRL 98, 066401 (2007); Steneteg et al., PRB 82, 075110 (2010); Niklasson JCTC 16 (2020) Idea: Use a very good guess for the wave function of the “next step”. Ensuring time-reversibility iskey for stability of propagation.19

THE AUTHORwhere Other Acceleration zz ) and T r[] 0.TheRamanintensityI(!)isthTechniques1re 3 ( xx as yy )andTr[] 0.TheRamanintensityI(!)isthenexpressedzzTHE AUTHORJ. Kolafa JCC 25, 335 (2004);T. Kühne, et al. PRL 98, 066401 (2007); Steneteg4 et al., PRB 82, 075110 (2010); Niklasson JCTC 16 (2020)(194)I(!) Iiso the“nextstep”.Ensuringtime-reversibilityis )andTr[] 0.TheRamanintensityI(!)isthenexpressedyyzz Z 1)I(!) Iiso Ianisokey for stability of propagation. 3Ni!tZ 1(195)Iiso dte (0) (t) N42 nsity from previous stepi!t1 ) I(!) Iiso IIanisodte (0) 1(t) iso Z32 1Ni!tZMoreinvolved: 1Z(196)I dte Tr[(0)·(t)] .aniso 1Ni!tN2 1 i!tI dte (0) (t) iso) ExtrapolateI dte Tr[(0)·(t)] edon information ofaniso2 12 X1previousZtimesteps 1 c'Niijji!tX IanisoFormulate dte Tr[(0)·(t)] ariables(density,wavej cij 'j2 1i functions) that move in a harmonic potential around the self-consistentsolution C'jXTcij 'j C'P C Ci TjP C C S Density(197)!s h' 'iijijMatrixBasis Sets C')S ! sij h'i 'j iTP C C1( xxyy319

THE AUTHORwhere Other Acceleration zz ) and T r[] 0.TheRamanintensityI(!)isthTechniques1re 3 ( xx as yy )andTr[] 0.TheRamanintensityI(!)isthenexpressedzzTHE AUTHORJ. Kolafa JCC 25, 335 (2004);T. Kühne, et al. PRL 98, 066401 (2007); Steneteg4 et al., PRB 82, 075110 (2010); Niklasson JCTC 16 (2020)(194)I(!) Iiso the“nextstep”.Ensuringtime-reversibilityis )andTr[] 0.TheRamanintensityI(!)isthenexpressedyyzz Z 1)I(!) Iiso Ianisokey for stability of propagation. 3Ni!tZ 1(195)Iiso dte (0) (t) N42 nsity from previous stepi!t1 ) I(!) Iiso IIanisodte (0) 1(t) iso Z32 1Ni!tZMoreinvolved: 1Z(196)I dte Tr[(0)·(t)] .aniso 1Ni!tN2 1 i!tI dte (0) (t) iso) ExtrapolateI dte Tr[(0)·(t)] edon information ofaniso2 12 X1previousZtimesteps 1 c'Niijji!tX IanisoFormulate dte Tr[(0)·(t)] ariables(density,wavej cij 'j2 1i functions) that move in a harmonic potential around the self-consistentsolution C'jXTcij 'j C'P C Ci TjP C C S Density(197)!s h' 'iijijMatrixBasis Sets C')S ! sij h'i 'j mswithlimitedsuccess.Ongoingefforts! P C C1( xxyy319

Temperature control: the canonical ensemble The idea: couple the system to a thermostat (heat bath) Interesting because: Experiments are usually done at constant temperature Better modeling of conformational changesBathEnergy isconservedSystemEnergy isnot conserved20

1Nosé-Hoover thermostatṘ p/m(42)S. Nosé, J. Chem. Phys. 81, 511 (1984) & W. G. Hoover, Phys. Rev. A 31, 1695 (1985).Extended Hamiltonian (or Lagrangian):HN H2X p2p I V (R) 3N kB T 2MI2Q(43)IṗI FI p I(44)21

1Nosé-Hoover thermostatṘ p/m(42)S. Nosé, J. Chem. Phys. 81, 511 (1984) & W. G. Hoover, Phys. Rev. A 31, 1695 (1985).Extended Hamiltonian (or Lagrangian):HN H2X p2p I V (R) 3N kB T 2MI2Q(43)IOriginal systemṗI FIOscillator p IFictitious(44)21

1000dt A(t )B(t t ) (40)1 hA(0)B(t)i T0Nosé-Hoover thermostatṘ p/mZ 81,S. Nosé, J. Chem. Phys.t 511 (1984) & W. G. Hoover, Phys. Rev. A 31, 1695 (1985).ṗ Fd K(t )p( ) (t)Extended Hamiltonian(orLagrangian):1HN H 2Ṙ p/mX p2p I V (R) 3N kB T 2MI2QIX2pI2p HN H Original system V (R) FictitiousṗI 3NFIkOscillator T p IB2MI2QIp pIMomenta are damped by fictitious oscillator: ṗI FIQF I r RI VFI dT 1T̄(42)(41)(42)(43)(43) (44)(44)(45)ṘI (t)(46)h (t) (0)i 2kB T (t)rT (t) dtT (t)2 (t) T̄3T̄ N (47)(48)21

1000dt A(t )B(t t ) (40)1 hA(0)B(t)i T0Nosé-Hoover thermostatṘ p/mZ 81,S. Nosé, J. Chem. Phys.t 511 (1984) & W. G. Hoover, Phys. Rev. A 31, 1695 (1985).ṗ Fd K(t )p( ) (t)Extended Hamiltonian(orLagrangian):1HN H 2Ṙ p/mX p2p I V (R) 3N kB T 2MI2QIX2p 2pIHN H Original system V (R) FictitiousṗI 3NFIkOscillator T p IB2MI2QIp pIMomenta are damped by fictitious oscillator: ṗI FIQF I r RI Vof phase space Ergodicity problems - system may be stuck inFa regionṘI (t)I Possible solution: Nosé-Hoover chainsh (t) (0)i 2k T (t)(42)(41)(42)(43)(43) (44)(44)(45)(46)(47) to theAttach another fictitious oscillator to the first, and anotherto the second, and anotherBr third, . (chain of fictitiousoscillators)dTT (t) dtT (t)T̄ 1T̄ 2 (t)Martyna, Klein,Tuckerman, J. Chem. Phys. 97,(48)2635 (1992)3T̄ N 21

How to model a thermostat: first ideas Temperature rescaling: Berendsen “thermostat” Rescale velocities by a factor containing the ratio of target and instant temperatureDoes not sample the canonical ensemble (wrong temperature distribution)“Flying ice-cube” effect: rotations and translations acquire high Ekin and vibrations are frozenH. J. C. Berendsen, et al. J. Chem. Phys. 81 3684 (1984)22

How to model a thermostat: first ideas Temperature rescaling: Berendsen “thermostat” Rescale velocities by a factor containing the ratio of target and instant temperatureDoes not sample the canonical ensemble (wrong temperature distribution)“Flying ice-cube” effect: rotations and translations acquire high Ekin and vibrations are frozenH. J. C. Berendsen, et al. J. Chem. Phys. 81 3684 (1984) Simple stochastic idea: Andersen thermostat At each nth time-step, replace velocity of a random particle by one drawn from a MaxwellBoltzmann distribution at target temperature Not very efficient, no conserved quantityVery sensitive on nH. C. Andersen, J. Chem. Phys. 72, 2384 (1980)22

HN H II2MI V (R) 2Q 3N kB T Stochastic Velocity Rescaling(43)p G. Bussi, D. Donadio, and M. Parrinello,126, 014101ṗIJ. Chem. FPhys.pI (2007).(44)I QCombine concepts from velocity rescaling (fast!) with concepts from stochastic thermostatsF I r RI V(45)(accurate!)FI ṘI (t)(46)h (t) (0)i 2kB T (t)rT (t) dtT (t)2 (t) T̄3T̄ N (47)Target temperature follows a stochastic differential equation: dT 1T̄(48)323

HN H II2MI V (R) 2Q 3N kB T Stochastic Velocity Rescaling(43)p G. Bussi, D. Donadio, and M. Parrinello,126, 014101ṗIJ. Chem. FPhys.pI (2007).(44)I QCombine concepts from velocity rescaling (fast!) with concepts from stochastic thermostatsF I r RI V(45)(accurate!)FI ṘI (t)(46)h (t) (0)i 2kB T (t)rT (t) dtT (t)2 (t) T̄3T̄ N (47)Target temperature follows a stochastic differential equation: dT 1T̄Temperaturerescaling(48)White noise323

HN H II2MI V (R) 2Q 3N kB T Stochastic Velocity Rescaling(43)p G. Bussi, D. Donadio, and M. Parrinello,126, 014101ṗIJ. Chem. FPhys.pI (2007).(44)I QCombine concepts from velocity rescaling (fast!) with concepts from stochastic thermostatsF I r RI V(45)(accurate!)FI ṘI (t)(46)h (t) (0)i 2kB T (t)rT (t) dtT (t)2 (t) T̄3T̄ N (47)Target temperature follows a stochastic differential equation: dT 1T̄Temperaturerescaling (48)White noiseVery successful thermostat, weakly dependent on relaxation time τPseudo-Hamiltonian is conserved3Bussi, Parrinello, Phys. Rev. E 75, 056707 (2007)23

FI r RI VLangevin (stochastic)thermostatFI ṘI (t)(45)(46)22Xpp S. A. Adelman and J. D. Doll,J. Chem. Phys. 64, 2375 (1976).I 2kT (t) 3N kB T (47)HN H h (t) (0)i VB(R)2MIr2Q IdTT (t)dt via the LangevinT (t) equation:Modeldynamicsp 12 (t)(48)ṗ F pIII T̄T̄3T̄ N QMI R̈I FII ṘI (t)FI rRI V (49)FI ṘI (t)h (t) (0)i 2kB T (t)r dTT (t) dtT (t) 12 (t) T̄T̄3T̄ N 3MI R̈I FII ṘI (t)24

FI r RI VLangevin (stochastic)thermostatFI ṘI (t)(45)(46)22Xpp S. A. Adelman and J. D. Doll,J. Chem. Phys. 64, 2375 (1976).I 2kT (t) 3N kB T (47)HN H h (t) (0)i VB(R)2MIr2Q IdTT (t)dt via the LangevinT (t) equation:Modeldynamicsp 12 (t)(48)ṗ F pIII T̄T̄3T̄ N QMI R̈I FII ṘI (t)FI rRI V (49)ṘI (t)Original systemFriction andFWhiteI Noiseh (t) (0)i 2kB T (t)r dTT (t) dtT (t) 12 (t) T̄T̄3T̄ N 3MI R̈I FII ṘI (t)24

FI r RI VLangevin (stochastic)thermostatFI ṘI (t)(45)(46)22Xpp S. A. Adelman and J. D. Doll,J. Chem. Phys. 64, 2375 (1976).I 2kT (t) 3N kB T (47)HN H h (t) (0)i VB(R)2MIr2Q IdTT (t)dt via the LangevinT (t) equation:Modeldynamicsp 12 (t)(48)ṗ F pIII T̄T̄3T̄ N QMI R̈I FI I ṘI (t)FI rRI V (49)ṘI (t)Original systemFriction andFWhiteI Noiseh (t) (0)i 2kB T (t)r Sensitive on γdTT (t) dtT (t) 12 achievethe“best”critical T̄T̄3T̄N damping?3MI R̈I FIDisturbs dynamics considerablyI ṘI (t)24

Colored noise thermostatsM. Ceriotti, G. Bussi, M. Parrinello, JCTC 2010, 6, 1170-1180 (http://gle4md.org/index.html) Extremely flexible class of thermostats based on theGeneralized Langevin Equation (GLE)MarkovianṘ(no memory)processinhighdimensions p/m 0ṗV (R) ṡ0 pAp Bp ,sn array of uncorrelated Gaussian noises, V (R) is the gradient of the potential (that cae Ap and Bp are matrices that obey the relation0Ap C p TC p Ap TBp Bp ,the covariance matrix defined as Cp h(p, s) (p, s)i. By integrating out the s25dT

Colored noise thermostatsM. Ceriotti, G. Bussi, M. Parrinello, JCTC 2010, 6, 1170-1180 (http://gle4md.org/index.html) Extremely flexible class of thermostats based on theGeneralized Langevin Equation (GLE)MarkovianṘ(no memory)processinhighdimensions p/m 0ṗV (R) ṡ0 pAp Bp ,sextraOriginal systemFriction and White Noisefictitious degrees0of freedom Gaussian noises, V (R) is the gradient of theuncorrelatedn array ofe Ap and Bp are matrices that obey the relationAp C p TC p Ap potential (that caTBp Bp ,the covariance matrix defined as Cp h(p, s) (p, s)i. By integrating out the s25dT

p(t t) p(t t/2) F (t t)2Colored noise thermostats(36)ZZM. Ceriotti, G. Bussi, M. Parrinello, JCTC 2010,(http://gle4md.org/index.html)1 6, 1170-11803N3NH/kB ThAi dRdpeA(p,R)(37)Extremely flexible classofthermostatsbasedontheZZZGeneralizedLangevin Equation (GLE)13N3NH/kB ThA(0)B(t)i y)processinhighdimensions Ṙ p/mZ Z0ṗV (R)p1 T000 A B, p p A(p(t ), R(t )) (39)dthAiṡ0sT0extraOriginal systemFriction andZ WhiteNoiseT1fictitious degrees0000hA(0)B(t)i dtA(t)B(t t)(40)offreedomn array of uncorrelated Gaussian noises, V (R) is thegradientofthepotential(thatcaT 0e Ap and Bp arematrices m (integrating out s):Z tṗ Fd K(t )p( ) (t)(41)TTAp C1p Cp Ap Bp Bp ,Ṙ p/m(42)Tthe covariance matrix defined as Cp h(p, s) (p, s)i. By integrating out the s25d

p(t t) p(t t/2) F (t t)2Colored noise thermostats(36)ZZM. Ceriotti, G. Bussi, M. Parrinello, JCTC 2010,(http://gle4md.org/index.html)1 6, 1170-11803N3NH/kB ThAi dRdpeA(p,R)(37)Extremely flexible classofthermostatsbasedontheZZZGeneralizedLangevin Equation (GLE)13N3NH/kB ThA(0)B(t)i y)processinhighdimensions Ṙ p/mZ Z0ṗV (R)p1 T000 A B, p p A(p(t ), R(t )) (39)dthAiṡ0sT06THE AUTHORextraOriginal systemFriction andZ WhiteNoiseT1fictitious degrees0000hA(0)B(t)i dtA(t)B(t t)(40)offreedomn array of uncorrelated Gaussian noises, V (R) is thegradientofthepotential(thatcaT 0e Ap (67)and Bp arematrices m (integrating ẋout s):p/mZ tt Vṗ F ṗ d K(t )p( ) (t)(41)(68)K(t)p( ) (t)TTAp C1p xCA BB,ppppMemory KernelColored NoiseṘ p/m(42)T(t)s)(0) kB TByK(t)Dissipation:the (69)covarianceFluctuationmatrix definedasH(t)Cp h(p,(p, s)i.integrating out the s25d

Pressure control: Isobaric-isothermic ensemble Definition of instantaneous pressure:Stress TensorSystem26

Pressure control: Isobaric-isothermic ensemble Definition of instantaneous pressure:Stress Tensor Similar schemes as thermostats: pressure rescaling, extended Lagrangian, stochasticpressure rescalingParinello and Rahman, J. Appl. Phys 52, 7182 (1981);Bussi, Zykova-Timan, Parrinello, J. Chem. Phys. 130, 074101 (2009)System26

Pressure control: Isobaric-isothermic ensemble Definition of instantaneous pressure:Stress Tensor Similar schemes as thermostats: pressure rescaling, extended Lagrangian, stochasticpressure rescalingParinello and Rahman, J. Appl. Phys 52, 7182 (1981);Bussi, Zykova-Timan, Parrinello, J. Chem. Phys. 130, 074101 (2009) Use thermostat together with a barostat to control pressure and temperatureBathSystem26

The i-PI wrapper code: sockets interface The i-PI code is a python interface for the calculation of (path integral ab initio) molecular dynamics Interfaced to many electronic structure codes (check our webpage!)i-PI communicates with the electronic structure codes through internet (or UNIX) sockets, making itextremely flexible27

The i-PI wrapper code: sockets interface The i-PI code is a python interface for the calculation of (path integral ab initio) molecular dynamics Interfaced to many electronic structure codes (check our webpage!)i-PI communicates with the electronic structure codes through internet (or UNIX) sockets, making itextremely flexiblehttp://ipi-code.org/27

Thei-PIwrappercode:socketsinterfaceA Python interface for (path integrals) (ab initio) moleculardynamicsExtreme implementation of the “black box” concept. Positions andforces are communicated over internet (or UNIX) smo.github.io/gle4md/index.html?page ipi7Michele Ceriotti - EPFL - COSMONuclear quantum effects made easy with i-PI28

Thei-PIwrappercode:socketsinterfaceA Python interface for (path integrals) (ab initio) moleculardynamicsExtreme implementation of the “black box” concept. Positions andforces are communicated over internet (or UNIX) smo.github.io/gle4md/index.html?page ipi7Michele Ceriotti - EPFL - COSMONuclear quantum effects made easy with i-PI28

The whole pointof using sockets is owever, one needs to make sure that the client and server can talkThe pointtoofeachusing other.sockets is to enable remote execution — server must talk to .In HPC systems it might not be trivial13Michele Ceriotti - EPFL - COSMONuclear quantum effects made easy with i-PI29

The whole pointof using sockets is owever, one needs to make sure that the client and server can talkThe pointtoofeachusing other.sockets is to enable remote execution — server must talk to .In HPC systems it might not be trivial13Michele Ceriotti - EPFL - COSMONuclear quantum effects made easy with i-PI29

The whole pointof using sockets is owever, one needs to make sure that the client and server can talkThe pointtoofeachusing other.sockets is to enable remote execution — server must talk to .In HPC systems it might not be trivial13Michele Ceriotti - EPFL - COSMONuclear quantum effects made easy with i-PI29

Molecular Crystals: Free Energy DifferencesΔ Ecryst/N-EmolForce FieldRossi, Gasparotto, Ceriotti, PRL 117, 115702 (2016)30

Molecular Crystals: Free Energy Differences(178)(179)(180) Z1el0" Eelel ikB TXln [kB T /( !i )]Ndof kB T2#E0i0C(182)F( ) ZF (183)(184)(185)Force FieldĤel 0X2 !i g0dgcoth( !i g/2)g40CH(181)elnucleiRossi, Gasparotto, Ceriotti, PRL 117, 115702 (2016)Δ Ecryst/N-EmolharmFc!q CHCA1baln Z@Fd@V ( ) VF F (1)VDebyeZ 1 Z 1@V ( ) d hVF F VDebye i d@0030

Z1"0X2 !i g0dgcoth( !i g/2)g4(180)MolecularCrystals: FreeEnergyDifferences0harmFc!qΔ Ecryst/N-EmolCH0(181)iNdof kB T2Rossi, Gasparotto, CeriottiX, PRL 117, 115702 (2016) ECHE0 kB Tln [kB T /( !i )]#E0i(182)F( ) ZF (183)(184)0C(185)CHCA(186)(187)(188)Force Field CA1baln Z@Fd@V ( ) VF F (1)VDebyeZ 1 Z 1@V ( ) d hVF F VDebye i d@00V ( ) VDF T (1)VF FZ 1 Z 1@V ( ) d hVDF T VF F i d@0030

(183)F d@Molecular Crystals: Free Energy Differencesa(184)Δ Ecryst/N-EmolV ( ), VF F115702 (1(2016))VDebyeRossi, Gasparotto, CeriottiPRL 117,CHCA(185)(186)(187)(188)0C(189) CA Z0 1@V ( )@Force Field10hVF FVDebye i dV ( ) VDF T (1)VF FZ 1 Z 1@V ( ) d hVDF T VF F i d@00(190)(191)d ZCAQA Z1m0 !0 m!1hK(µ)iµ 3kB T /2dµµ30

Molecular Crystals: Free Energy DifferencesΔ Ecryst/N-EmolRossi, Gasparotto, Ceriotti, PRL 117, 115702 (2016)0CForce FieldAb initio (DFT)30

11@V ( )(185) d hVF F VDebye i d@00Molecular Crystals: Free Energy Differences(186)Rossi, Gasparotto, Ceriotti, PRL 117,

Tutorial Series: Ab initio molecular dynamics Mariana Rossi Max Planck Institute for the Structure and Dynamics of Matter Fritz Haber Institute of the MPG