Ab-Initio Molecular Dynamics

Thomas D. Kühne

Source-linked AI summary

Ab-Initio Molecular Dynamics

Thomas D. Kühne

arXiv:1201.5945v2Jan 28, 2012physics.chem-phcond-mat.mtrl-scicond-mat.softcond-mat.stat-mechphysics.comp-ph

TL;DR

Classical molecular dynamics is limited by empirical-force transferability, while ab-initio molecular dynamics improves predictive power at substantial computational cost. This review develops and surveys Born-Oppenheimer, Car-Parrinello, and second-generation Car-Parrinello molecular dynamics, reporting efficient simulations on medium-sized systems over nanosecond timescales.

Problem
Empirical force fields can lack transferability beyond the systems and phase-diagram regions used for parametrization, motivating accurate on-the-fly electronic-structure forces despite their computational cost.
Method
The review derives Born-Oppenheimer molecular dynamics and presents Car-Parrinello and second-generation Car-Parrinello approaches, including orbital transformations that preserve idempotency during electronic minimization.
Results
Second-generation Car-Parrinello molecular dynamics combines Born-Oppenheimer accuracy and long time steps with Car-Parrinello efficiency, achieving system-dependent speed-ups of one to two orders of magnitude.
Takeaways & Limitations
The method enables ab-initio molecular dynamics for medium-sized systems of up to a few thousand atoms over timescales of a couple of nanoseconds.
Takeaways & Limitations
Density-functional-theory calculations require approximations because the universal kinetic-plus-electron-electron functional is not known explicitly.

Abstract

from arXiv · show

Computer simulation methods, such as Monte Carlo or Molecular Dynamics, are very powerful computational techniques that provide detailed and essentially exact information on classical many-body problems. With the advent of ab-initio molecular dynamics, where the forces are computed on-the-fly by accurate electronic structure calculations, the scope of either method has been greatly extended. This new approach, which unifies Newton's and Schrödinger's equations, allows for complex simulations without relying on any adjustable parameter. This review is intended to outline the basic principles as well as a survey of the field. Beginning with the derivation of Born-Oppenheimer molecular dynamics, the Car-Parrinello method and the recently devised efficient and accurate Car-Parrinello-like approach to Born-Oppenheimer molecular dynamics, which unifies best of both schemes are discussed. The predictive power of this novel second-generation Car-Parrinello approach is demonstrated by a series of applications ranging from liquid metals, to semiconductors and water. This development allows for ab-initio molecular dynamics simulations on much larger length and time scales than previously thought feasible.

INTRODUCTION

Molecular dynamics provides equilibrium and dynamical information by numerically solving Newton’s equations, but its force models and classical treatment limit predictive scope. AIMD addresses these limitations by computing forces from electronic structure calculations, although finite computational resources still restrict accessible scales.

MD numerically solves Newton’s equations to compute equilibrium thermodynamic and dynamical properties at finite temperature.
Empirical potentials are difficult to transfer beyond their fitting systems and cannot reliably describe many chemical-bonding processes.
AIMD computes interatomic forces on-the-fly from accurate electronic structure calculations, potentially removing these empirical-potential limitations.
MD can evaluate thermal averages through temporal averaging under the ergodicity hypothesis, while also exposing nuclear real-time evolution.
Classical nuclei are usually adequate, but very light atoms or low temperatures may require quantum treatments such as imaginary-time path integrals.
Finite simulation resources restrict length and time scales; periodic systems capture only correlations much smaller than L and relaxation times much shorter than T.

AN AB-INITIO POTENTIAL

The Born-Oppenheimer potential treats electrons as being in instantaneous equilibrium with fixed nuclei, reducing nuclear dynamics to an electronic ground-state energy problem. Its direct many-body formulation is computationally prohibitive because the wavefunction depends on all electronic coordinates.

Under the Born-Oppenheimer approximation, the electronic Hamiltonian depends parametrically on nuclear positions while electrons remain in instantaneous equilibrium with the nuclei.
The electronic ground state is obtained from a high-dimensional eigenvalue problem involving the many-body wavefunction and its energy.
A real-space grid with 100 points per coordinate would require 10^(6Ne) points for Ne electrons, making direct many-body solutions impractical.

Density Functional Theory

Density Functional Theory replaces the many-electron wavefunction with the electron density and maps the interacting problem onto a fictitious single-particle Kohn-Sham system. Its practical use depends on approximating the unknown universal energy functional and exchange-correlation contributions.

The Hohenberg-Kohn theorem establishes a one-to-one mapping between ground-state density and external potential, making density the central variable.
The ground-state problem can be solved by self-consistent diagonalization or by minimizing the quantum expectation value.
The density minimization must satisfy N-representability, whereas general v-representability has no known solution.
DFT expresses kinetic, electron-electron, and electron-ion energies as functionals of the density, but the universal kinetic-plus-interaction functional is unknown.
Thomas-Fermi approximations neglect many-body correlation effects, while exchange-correlation functionals are introduced to account for them.
Kohn-Sham orbitals are fictitious orbitals without strict physical meaning, except in specific exact-functional isolated-system cases.
The Kohn-Sham scheme maps the interacting many-body problem onto a fictitious single-particle system with an effective potential.

The Exchange and Correlation Functional

The exchange–correlation functional is formally required in DFT but is not known exactly beyond the uniform electron gas, so practical AIMD uses approximations. Exchange can be treated explicitly with orbitals, whereas correlation remains without an exact orbital or density expression.

The exact exchange–correlation functional is unknown except for the uniform electron gas, requiring approximate treatments in practical DFT.
Exchange energy can be calculated exactly using an explicit orbital functional, including with Kohn–Sham orbitals.
Hartree–Fock exchange is computationally costly for periodic systems because its nonlocal form requires four-center integrals.
Correlation energy has no exact expression in terms of either orbitals or densities.
Because exchange–correlation energy is usually smaller than the remaining known energy terms, simple approximations may still yield qualitatively correct ground-state energies without adjustable parameters.

Born-Oppenheimer Molecular Dynamics

Born–Oppenheimer molecular dynamics minimizes the electronic energy at every nuclear configuration under orbital orthonormality constraints. Its force expression contains Hellmann–Feynman, Pulay, and non-self-consistent contributions, making accurate force evaluation more demanding than energy evaluation.

BOMD minimizes the electronic energy at every molecular-dynamics step subject to orbital orthonormality.
The BOMD force includes Hellmann–Feynman, Pulay, and non-self-consistent contributions.
Pulay forces arise from orbital orthonormality constraints when basis functions depend explicitly on nuclear positions.
Neglecting Pulay or non-self-consistent forces produces inconsistent forces because numerical calculations are not exactly self-consistent.
Force errors depend linearly on electronic-density errors, whereas energy errors do not, making accurate forces more demanding than accurate energies.
Born–Oppenheimer decoupling permits integration steps up to the nuclear resonance limit and, in principle, allows metals to be treated regardless of band gap.

Car-Parrinello Molecular Dynamics

Car–Parrinello molecular dynamics propagates electronic orbitals and nuclei together using a fictitious electronic inertia, avoiding full self-consistent minimization at every step. Its efficiency depends on maintaining adiabatic separation, while metallic systems require additional treatment and the BOMD–CPMD choice remains application-dependent.

CPMD couples electron and ion dynamics by treating electronic degrees of freedom as classical variables with fictitious mass parameter µ.
CPMD reduces nuclear-force cost because no self-consistent-field cycle is required at every step.
The electronic and ionic frequencies must remain separated, requiring the highest ionic phonon frequency to be much smaller than the lowest electronic frequency.
The maximum joint timestep depends on the fictitious inertia µ and must remain below the inverse electronic frequency.
Choosing µ trades computational efficiency against deviation from the instantaneous Born–Oppenheimer surface.
CPMD eliminates non-self-consistent forces by evaluating the instantaneous electronic state without fully minimizing the energy.
CPMD timesteps are roughly two orders of magnitude larger than Ehrenfest timesteps but still about one order shorter than BOMD timesteps.
Metallic CPMD requires an electronic thermostat or fractional occupations, and preferring CPMD over BOMD depends on accuracy criteria and application.

SECOND GENERATION CAR-PARRINELLO MOLECULAR DYNAMICS

DFT-based AIMD remains limited by high computational cost and short integration timesteps, constraining the length and time scales accessible to simulations.

High DFT-based AIMD cost limits attainable simulation length and time scales despite substantial progress.
Short integration timesteps remain a limitation of AIMD techniques.

An Efficient and Accurate Car-Parrinello-like

Second-generation CPMD combines BOMD’s accuracy and long time steps with CPMD’s efficiency by using coupled electron-ion dynamics that keeps electrons near the instantaneous Born-Oppenheimer surface without an artificial mass parameter.

Second-generation CPMD combines BOMD accuracy and long time steps with CPMD efficiency.
The method replaces fictitious electronic Newtonian dynamics with coupled electron-ion dynamics that keeps electrons close to the instantaneous BO surface.It does not require an artificial mass parameter.
Its efficiency improves by one to two orders of magnitude depending on the system.The passage reports demonstrations across a wide range of applications.
The approach propagates the smoother density kernel PS instead of the more rapidly varying coefficient matrix C.This formulation applies within mean-field electronic structure theories using a nonorthogonal basis.

Coupled Electron-Ion Dynamics

The coupled electron-ion scheme predicts electronic states from previous density-kernel matrices and corrects them efficiently, maintaining time reversibility while closely approaching the instantaneous ground state.

The electronic dynamics is specified by a predictor-corrector integrator rather than fictitious electron dynamics from a modified Lagrangian.Electronic short-term integration requires a highly accurate and efficient algorithm.
The predictor uses previous PS matrices because the density kernel evolves more smoothly and is easier to predict than C.The predictor is followed by a corrector that reduces error and deviation from the instantaneous ground state.
The electron dynamics is accurate and time reversible up to O(∆t^2K−2), with ω chosen for stable relaxation toward the instantaneous ground state.
The scheme requires only one preconditioned electronic gradient calculation per AIMD step in general.Repeated prediction can approach the ground state more closely but requires additional electronic force calculations.
It avoids both the SCF cycle and iterative wavefunction optimization while permitting time steps as large as in standard BOMD.

Electronic Forces by Orbital Transformations

Orbital transformations parameterize occupied orbitals through an auxiliary variable while preserving idempotency, enabling efficient electronic minimization and approximate force evaluation near the ground state.

The orbital transformation method introduces an auxiliary variable X to parameterize the occupied orbitals.
The constraint X^T S C^p(t_n)=0 makes every finite step satisfy the idempotency condition when the reference orbitals are orthonormal.Minimization therefore occurs in a linear auxiliary tangent space.
The approximate energy functional uses the predicted density ρ^p(r) and provides nuclear forces from its analytic ionic-coordinate gradient.An additional force term appears because the corrected and predicted densities differ.
A single preconditioned minimization step leaves C(t_n) as an approximate eigenfunction, producing insignificant force error when it remains close to the ground state.

Modified Langevin Equation

The electronic propagation can make nuclear dynamics dissipative, so a modified Langevin equation introduces damping and fluctuation noise to recover canonical sampling, with the damping bootstrapped from temperature.

Despite proximity to the instantaneous ground state, the nuclear dynamics is dissipative because the electron propagation scheme is likely not symplectic.
The modified Langevin equation uses damping and additive white noise obeying the fluctuation-dissipation theorem to sample the canonical distribution.
The dissipative force is modeled with an intrinsic damping coefficient γ_D under the assumption that energy dissipation is exponential.
The unknown damping coefficient need not be known beforehand and can be bootstrapped by requiring the correct average temperature from equipartition.
Illustrative Examples: Liquid Silicon, Silica and: For liquid SiO2, the energy shift is 4.16×10−4 Hartree per atom with one corrector step and 3.5×10−5 Hartree per atom with two.The average mean-force deviation is unbiased.

Water

The method is tested on liquid silicon, silica, and water using a single electronic corrector step, with results close to reference calculations. It accurately reproduces structural and dynamical behavior while substantially improving simulation efficiency.

Water: Calculations cover liquid metallic silicon, silica, and water to test performance across different band gaps, system sizes, and liquid types.The systems represent liquid metals, complex polarizable ionic liquids, and hydrogen-bonded fluids.
Water: The reported simulations use experimental liquid densities, TZV2P basis sets, norm-conserving pseudopotentials, GGA exchange-correlation, and Γ-point sampling.The Langevin friction coefficient values were determined using the stated integration procedure.
Water: A single corrector step keeps energies slightly and nearly constantly above the electronic ground-state surface while allowing the deviation to be controlled by additional corrector steps.The reported production simulations use only one preconditioned electronic gradient calculation.
Water: The method remains close to Born-Oppenheimer molecular dynamics reference results for liquid systems, including metallic liquid silicon, where ordinary Car-Parrinello schemes are problematic.A single preconditioned gradient calculation is sufficient, and difficult cases show a speed-up of two orders of magnitude over pure extrapolation.
Water: Velocity autocorrelation functions and their Fourier transforms agree with reference calculations and experiment, demonstrating accurate dynamical-property simulations.The calculations are shown for water at 325 K.
Water: A 1 ns liquid Si64 trajectory produces a Maxwellian kinetic-energy distribution, supporting canonical sampling despite the stochastic dynamics.The long trajectory reduces noise and samples the tails of the kinetic-energy distribution.
Water: Across the reviewed applications, second-generation CPMD enables simulations of medium-sized systems up to a few thousand atoms for as long as a couple of nanoseconds.The review also reports successful structure relaxation through dynamic annealing and geometry optimization.