Car Parrinello Molecular Dynamics Nicholas Walker Jrg Hutter

Car Parrinello Molecular Dynamics Nicholas Walker Jürg Hutter, “Introduction to Ab Initio Molecular Dynamics. ” Physical Chemistry Institute, University of Zurich. Winterthurerstrasse 190 Dominik Marx and Jürg Hutter, “Ab initio molecular dynamics: Theory and Implementation. ” Modern Methods and Algorithms of Quantum Chemistry, J. Grotendorst (Ed. ), John von Neumann Institute for Computing, Jülich, NIC Series, Vol. 1, ISBN 3 -00 -005618 -1, pp. 301 -449, 2000.

Introduction – Ab Initio Molecular Dynamics • Chemically complex systems are not well-suited for classical MD • Many different types of atoms • Qualitative changes in electronic structure • Ab initio MD relies on DFT (Kohn-Sham) • Electronic variables explicitly considered • Not integrated out beforehand • Treated as active degrees of freedom • Emergent properties can be observed easily • Tracing back behavior to a specific mechanism is difficult

Introduction – Why Do We Care? • Ab initio molecular dynamics is accurate, but slow • Electronic structure problem is difficult • Smaller timesteps are used • Born-Oppenheimer method computationally complex • Recalculate electronic structure problem at every timestep • Car Parrinello method avoids recalculating the electronic structure • Considerable speedup is gained • But at what cost and with what challenges?

Relevant Kohn-Sham Equations • Kohn-Sham Energy • Charge Density • Exchange-Correlation Energy

Born-Oppenheimer MD • The BO approximation allows for the separation of the ionic and electronic systems • The ionic interaction energy is the same in BOMD as in Kohn-Sham • The Lagrangian of the system is thus expressed as

BOMD Scheme • Verlet algorithm for dynamics v[: ] = v[: ]+dt/(2*m[: ])*f[: ] r[: ] = r[: ]+dt*v[: ] Optimize Kohn-Sham orbitals Calculate forces v[: ] = v[: ]+dt/(2*m[: ])*f[: ] dt – timestep m – ion mass r – ion position v – ion velocity f – force acting on ion

Car Parrinello MD • Exploit time-scale separation of fast electronic and slow nuclear motion • Classical mechanical adiabatic energy-scale separation • Map two component quantum/classical problem to two component classical problem • Separate energy scales • Lose explicit time-dependence of quantum subsystem

CPMD Lagrangian • Extended energy functional to introduce orthonormality constraint • Orbitals considered as classical fields in Lagrangian • Resulting equations of motion

CPMD Temperature • The nuclei evolve in time at an instantaneous physical temperature • Proportional to sum of nuclear kinetic energies (equipartition theorem) • The electrons evolve in time at a fictitious temperature • Proportional to sum of fictitious electronic kinetic energies (equipartition theorem) • Electrons are “cold” – close to instantaneous minimized energy (BO surface) • Ground state wavefunction optimized for initial configuration will stay close to the ground state during time evolution if it is at a sufficiently low temperature

CPMD Adiabacity • Separate nuclear and electronic motion • Electronic subsystem must stay cold for a long time • Electronic subsystem must follow slow nuclear motion adiabatically • Nuclei still kept at higher temperature • Achieved through decoupling of the two subsystems and adiabatic time evolution • Power spectra of both dynamics must not have too much overlap in the frequency domain • Energy transfer between “hot” nuclei and “cold” electrons becomes practically impossible

CPMD Controlling Adiabacity • Adiabatic separation satisfied by large frequency gap • Frequency spectrum of orbital classical fields close to the minimum (ground state) • Both the nuclei frequency spectrum and the smallest energy gap are determined by the system • Only control parameter is the fictitious electronic mass • Decreasing the mass shifts the frequency spectrum up, but also stretches it

CPMD Timescale • The maximum possible frequency determined by the cutoff frequency is also shifted up by lowering the electronic mass • This imposes an arbitrary condition on the maximum possible molecular dynamics time step • Because of this, compromises must be made on the control parameter • One work-around is to increase the nuclear masses to depress the nuclear vibrational frequency (at the cost of renormalization)

CPMD Forces and Constraints • The forces on the orbital fields are calculated as the action of the Kohn-Sham Hamiltonian • The forces with respect to the nuclear positions are the same as in BOMD • The constraint forces arise from the extended energy functional

CPMD Integrator • Constraint terms prohibit direct application of Verlet method • Assuming there is no nuclear position dependent overlap of the wavefunctions, only the orbital fields are constrained • Constraints are incorporated using the RATTLE algorithm v[: ] = v[: ]+dt/(2*m[: ])*f[: ] rp[: ] = r[: ]+dt*v[: ] dt – timestep Calculate Lagrange multiplier m – mass r[: ] = rp[: ]+L*r[: ] r – position Calculate forces rp – temporary position v[: ] = v[: ]+dt/(2*m[: ])*f[: ] v – velocity Calculate Lagrange multiplier f – force v[: ] = v[: ]+L*r[: ]

CPMD Performance Comparison With BOMD • Simulation of 8 Si atoms for 8 ps at 360 -370 K with energy cutoff of 10 Ry at the gamma point of a simple cubic supercell and fictitious mass of 400 au

Conclusions • CPMD is a useful method for running faster ab initio molecular dynamics simulations • • Fictitious electron dynamics Adiabatic energy-scale separation Electrons kept near ground state to retain accurate ionic forces Fictitious mass control parameter • BOMD can still be made as fast or faster than CPMD • Sacrifices energy conservation • CPMD provides far better energy conservation while still being fast • Small timesteps