Explain ab initio molecular dynamics like I'm five

Question

In Kashiwaya, S., Shi, Y., Lu, J. et al. Synthesis of goldene comprising single-atom layer gold. Nat. Synth (2024) monolayers of hexagonal Au(111) planes are produced and released into solution where they are observed to be stable in some cases. This followed theoretical work that showed they would likely be stable as long as supported by holes in sheets of graphene:

Inspired by the experimental realization of free-standing monolayer iron suspended in graphene pores, theoretical studies based on density-functional theory (DFT) and ab initio molecular dynamics (AIMD) predicted stable 2D Au membranes with or without confinement by graphene pores(13–17), specifically, stable Au membranes suspended in graphene pores with diameters up to 20 nm. Notably, the DFT structural optimization indicated that the favoured 2D crystal would be of densely packed hexagonal structure that corresponds to the bulk monolayer of Au(111) face-centred cubic (fcc) lattice.

I don't have a strong background in simulation. My current understanding is (in a nutshell):

• Molecular dynamics (in cases where actual dynamical behavior is to be modeled) is essentially implementing Newton's $F=ma$ for many bodies with an ODE-like solver, using an interatomic force fields usually obtained from some library, together with some thermal agitation or other forces of interest.
• DFT (in cases where static structural and electronic nature of (usually) periodic materials is to be modeled) provides an approximate quantum mechanical solution for the electron wave functions, allowing a calculation of potential energy. The nuclei are then moved slightly, energy is calculated again, and an algorithm finds the nuclear positions which minimize potential energy.

Before I can go off and dig in to the literature to find out exactly what ab initio molecular dynamics (AIMD) is and how it works, and how it differs from other techniques, it would be particularly helpful if I had some basics.

I'm guessing that the potential energy is still calculated with a DFT-like quantum mechanical approximation, but instead of a energy minimization routine for the nuclear positions, an ODE-like solver is used to calculate the trajectories of the nuclei as in molecular dynamics, and probably in this case a fluid bath at a finite temperature is included.

Have I got it right so far? Are there any suggested references written with an emphasis on clarity and understanding rather than mathematical rigor? (I know that for some those are synonymous, but for me it's much more helpful if I know where the math is going before I read it).

---

In a quick search I found Q-Chem 4.4 User’s Manual and while helpful, its target audience are those with a high degree of familiarity, not myself.

9.7 Ab Initio Molecular Dynamics

Q-Chem can propagate classical molecular dynamics trajectories on the Born-Oppenheimer potential energy surface generated by a particular theoretical model chemistry (e.g., B3LYP/6-31G* or MP2/aug-cc-pVTZ). This procedure, in which the forces on the nuclei are evaluated on-the-fly, is known variously as “direct dynamics”, “ab initio molecular dynamics” (AIMD), or “Born-Oppenheimer molecular dynamics” (BOMD). In its most straightforward form, a BOMD calculation consists of an energy + gradient calculation at each molecular dynamics time step, and thus each time step is comparable in cost to one geometry optimization step. A BOMD calculation may be requested using any SCF energy + gradient method available in Q-Chem, including excited-state dynamics in cases where excited-state analytic gradients are available. As usual, Q-Chem will automatically evaluate derivatives by finite-difference if the analytic versions are not available for the requested method, but in AIMD applications this is very likely to be prohibitively expensive.

While the number of time steps required in most AIMD trajectories dictates that economical (typically SCF-based) underlying electronic structure methods are required, any method with available analytic gradients can reasonably be used for BOMD, including (within Q-Chem) HF, DFT, MP2, RI-MP2, CCSD, and CCSD(T). The RI-MP2 method, especially when combined with Fock matrix and $Z$-vector extrapolation (as described below) is particularly effective as an alternative to DFT-based dynamics.

Answer 1

As you suspected, molecular dynamics (MD) is simply propagating Newton's II equation of motion ${\bf F}_i = m_i {\bf a}_i$ for the nuclei with masses $m_i$ at positions ${\bf r}_i$, velocities ${\bf v}_i = \dot{\bf r}_i$ and accelerations ${\bf a}_i = \dot{\bf v}_i = \ddot{\bf r}_i$. To solve the equations of motion, you need to know the initial positions ${\bf r}_i$ and velocities ${\bf v}_i$ at the time $t=0$. The solution is typically found by the leapfrog algorithm.

In both classical and ab initio MD, the forces can be computed from the (full set of) atomic positions, ${\bf F}_i = {\bf F}_i(\{{\bf r}_j\})$, as the negative gradient of the potential energy ${\bf F}_i = -\nabla_i E$, where we are differentiating with respect to the coordinates of the $i$th particle.

In classical MD, the energy (and thereby the forces as well) are given by some analytical expression; for example, the harmonic bond term $E(r) = \frac 1 2 k (r-r_0)^2$.

In ab initio MD, the energy is computed by describing the electrons quantum mechanically (Born-Oppenheimer approximation), which requires solving an approximate Schrödinger equation for them. What you get out is again an energy, which depends on the set of the classical coordinates of the nuclei $\{{\bf r}_i\}$, and efficient expressions exist for evaluating the force from the solution of the Schrödinger equation.

Note: MD simulations can be performed in various ensembles, which represent different physical environments; e.g. NVE is constant number of particles in constant volume and total energy (microcanonical ensemble), while NPT has constant pressure and temperature (isothermal-isobaric ensemble). Since systems are usually studied in interaction with the environment at constant pressure and temperature, the NPT ensemble is the typical choice. The ensemble is implemented in MD by adaptive adjustment of the velocities; various algorithms for temperature (thermostat) and pressure (barostat) control are available in the literature; see https://www.compchems.com/thermostats-in-molecular-dynamics/ and https://www.compchems.com/barostats-in-molecular-dynamics/, for example.

For further reading on MD, I suggest e.g. the GROMACS manual.