I am new to the world of Molecular Dynamics, and am curious to know exactly what is considered to be an ab initio Molecular Dynamics (AIMD) method, and how do they work?

The seminal work by Roberto Car and Michele Parrinello, titled "Unified Approach for Molecular Dynamics and Density-Functional Theory" was published 35 years ago!

I have recently come across the following "types" of MD methods while reading some research papers:

If I am missing methods, which I am sure I am, please feel free to add them through an answer! Also, it would be appreciated if one method is explained per answer, and is summarized in 2-3 paragraphs.

  • 2
    $\begingroup$ This might need to be split into two questions. You're mixing "AIMD" with "quantum MD". In one case we are going beyond Newton's laws and trying to include quantum forces for the MD. In the other case we're actually trying to solve the dynamics quantum mechanically, meaning solving the time-dependent Schroedinger equation. Classical MD on a quantum forcefield is different from quantum MD. I've commented our the "quantum MD" topics, they can be asked separately :) $\endgroup$ Commented Jul 15, 2020 at 19:38
  • $\begingroup$ Agree with @NikeDattani, Ehrenfest also introduces non-adiabatic effects so this further complicates things $\endgroup$
    – Cody Aldaz
    Commented Jul 15, 2020 at 19:43
  • $\begingroup$ Ahh ok. This is exactly what confused me! I am new to MD as a whole so idk what is ab-initio and why it is (or not) quantum MD under these contexts? That was my confusion if Quantum MD was for solving QD. Thanks for the edit @NikeDattani. I am willing to ask the second question as well. $\endgroup$
    – epalos
    Commented Jul 15, 2020 at 19:47
  • $\begingroup$ @NikeDattani I would want to save BOMD for someone else, because it should be fairly trivial to include. Would ab initio Ring Polymer/Path Integral MD be worth including or fit with this question? Technically, the ab initio part is basically the same as in BO/CPMD, but the nuclear part is done differently to approximate quantum effects. $\endgroup$
    – Tyberius
    Commented Jul 22, 2020 at 18:38
  • 1
    $\begingroup$ @EtiennePalos can you please make a hyperlink to a paper, for last three bulleted items in the list? It may help people find the papers and explain them. $\endgroup$ Commented Jul 24, 2020 at 17:23

3 Answers 3


CPMD: Car-Parrinello Molecular Dynamics

An approximation of BOMD (Born-Oppenheimer MD) where fictitious dynamics is used on the electrons to keep them close to their ground state, so that we do not have to keep solving for their ground state at every single step. We start with Newton's 2nd law (as does classical MD), but instead of the force being calculated by a full-fledged ab initio calculation at every step, the force itself has an EOM (equation of motion) which below is given by Eq. \eqref{eq:fictitious}. For one nucleus with position $\vec{r}$ and several electrons with orbitals $\{\psi_i\}$ we get:

\begin{align} \tag{1} \vec{F} &= m\vec{\ddot{r}} \\ - \nabla \, E\left[\{ \psi_i \} , \vec{r} \right] &= m\vec{\ddot{r}}\tag{2}\\ \mu \ddot{\psi}_i(\vec{r},t) &= - \frac{\delta E}{\delta \psi_i^*(\vec{r},t)} + \sum_j \Lambda_{ij} \psi_j(\vec{r},t),\tag{3}\label{eq:fictitious} \end{align}

where $\Lambda_{ij}$ is a matrix of Lagrange multipliers to allow satisfaction of the constraint that the wavefunctions $\psi_i$ must be orthogonal; and $E[\{\psi_i\},\vec{r}]$ is an energy functional (usually a Kohn–Sham energy one). For several nucleii, just make a new subscript for $\vec{r}$ and change the functional to $E[\{\psi_i\},\{\vec{r_I}\}]$, then the equations are exactly the same.

Implemented in:

  • CPMD: Literally named after the method! (open-source)
  • CP2K: Might also be named after the method! (open-source)
  • NWChem (open-source)
  • 2
    $\begingroup$ Interestingly, CP2K doesn't (or at least no longer) implement traditional CP, only 2nd generation CP. $\endgroup$
    – Tyberius
    Commented Jul 31, 2020 at 22:34
  • $\begingroup$ But does CP2K stand for "Car-Parinello 2000" ? $\endgroup$ Commented Jul 31, 2020 at 22:43
  • $\begingroup$ I just found where they explain the name: cp2k.org/faq:name $\endgroup$
    – Tyberius
    Commented Jul 31, 2020 at 22:44
  • $\begingroup$ I guessed correctly! Okay so CP2K can be removed from my list, if you suggest. $\endgroup$ Commented Jul 31, 2020 at 22:45
  • 1
    $\begingroup$ @AbdulMuhaymin-FreePalestine sorry you're right, I was accidentally looking at the CP2K links. Please feel free to edit my answer. $\endgroup$ Commented Apr 29 at 12:08

2nd Generation CPMD

Car-Parrinello MD avoids repeatedly solving the electronic problem by propagating the orbitals as if they were particles governed by Newton's equations. This is much more efficient than having to solve at each time step as is done in Born-Oppenheimer MD, though at the cost of decreasing the maximum timestep for the dynamics (too large a step will lose the ground state), slightly reduced accuracy (not exactly at the ground state for each time step) and introducing a spurious "mass parameter" to describe the electronic motion.

To address these problems, Thomas Kuhne et al. developed the "Car-Parrinello like approach to BOMD", also referred to as second-generation CPMD. The key differences of this approach are:

  • Rather than propagating the orbitals (or rather the MO coefficients $\mathbf{C}$), 2nd-CPMD propagates the density $\mathbf{P}$ (or $\mathbf{PS}$ for nonorthogonal orbitals). The density seems to evolve more smoothly than the coefficients, making it easier to work with.
  • A predictor-corrector method (in their paper, the Always Stable Predictor Corrector (ASPC) method, but in principle any such method) is used to propagate the density. This generates a prediction of the next coefficients $\mathbf{C}^p(t_n)$ based on the previous $K$ density matrices. A corrected set of coefficients $\mathbf{C}(t_n)$ is then formed as $$\mathbf{C}(t_n)=\omega \text{MIN}[\mathbf{C}^p(t_n)]+(1-\omega)\mathbf{C}^p(t_n)$$ $$\omega=\frac{K}{2K-1}$$ where $\text{MIN}$ is a minimization and $K$ is a parameter which determines the accuracy $O(\Delta t^{2K-2})$. This update procedure eliminates the need for the mass parameter.
  • The nuclear dynamics become dissipative, possibly due to the nonsymplectic electron dynamics. This is corrected by performing a short validation run of dynamics to compute a damping coefficient $\gamma$ for the system.

There are some finer details to make this scheme work (specific parameterization of $\mathbf{C}$, how $\gamma$ is obtained) which are available in the original papers (both of which have arXiv preprint versions available).


  • CP2K


  1. Thomas D. Kühne, Matthias Krack, Fawzi R. Mohamed, and Michele Parrinello Phys. Rev. Lett. 98, 066401 DOI: 10.1103/PhysRevLett.98.066401
  2. Kühne, T.D. (2014), Second generation Car–Parrinello molecular dynamics. WIREs Comput Mol Sci, 4: 391-406. DOI: 10.1002/wcms.1176
  • $\begingroup$ I forgot about the bounty and didn't get a reminder as I thought I would. Barely remembered, and manged to award it with a total of 11 minutes to spare! $\endgroup$ Commented Jul 28, 2020 at 15:59

ab initio Ehrenfest Dynamics

From Li et.al.,2005, JCP

"The Born Oppenheimer (BO) and extended Lagrangian (EL) trajectories are founded on the assumption that a single electronic potential surface governs the dynamics. .. A major limitation of adiabatic trajectories is that they are not applicable to reactions involving nonadiabatic electronic processes, i.e., multiple potential-energy surfaces." (Ex: Conical Intersections)

Potential Energy Surfaces of multiple electronic states

To account for electronic adiabaticity, we solve the full time-dependent Schrödinger equation for both nuclear and electronic degrees of freedom. In Ehrenfest methodology, the adiabatic potential energy surface

$$ E_{eff} = \langle\Phi|\hat{H_{el}}|\Phi\rangle = \sum_i |a_i|^2 E_i^2$$ Thus, the atoms evolve on an effective potential representing an average over the adiabatic states weighted by their state populations $|a_i|^2$. The method is also therefore referred to as mean-field approach. As a comparison, for BOMD/Ehrenfest dynamics:

  • BOMD

$$\hat{H}_{el} (\mathbf{r}; \mathbf{R}) \Phi_k(\mathbf{r}; \mathbf{R}) = E^{el}_k(\mathbf{R})\Phi_k(\mathbf{r}; \mathbf{R})$$

$$M_I\ddot{\mathbf{R}}_I =-\nabla_IE^{el}_k(\mathbf{R})=-\mathop{\nabla_I}_{\text{min }\Phi_k}\langle\Phi_k|\hat{H}_{el}|\Phi_k\rangle$$

The electronic wavefunction $\Phi_k(\mathbf{r}; \mathbf{R})$ is static (only implicit time-dependence) and the nuclear degrees of freedom are handled classically. The nuclear degrees of freedom are decoupled from electronic degrees of freedom, while for each MD step the electronic wavefunction has to be optimized for the ground state.

  • Ehrenfest dynamics

$$i\hbar\frac{\partial \Phi(\mathbf{r};\mathbf{R},t)}{\partial t}= \hat{H}_{el} (\mathbf{r}; \mathbf{R}) \Phi(\mathbf{r};\mathbf{R},t) $$

$$M_I\ddot{\mathbf{R}}_I =-\nabla_I\langle\hat{H}_{el}(\mathbf{r};\mathbf{R})\rangle$$

Here we have an explicit time dependence of the electronic wavefunction. Electronic and nuclear time evolutions are propagated with a three-time-step integrator. The electronic wavefunction is evolved via TD-SCF approach.

  • $\begingroup$ +1. A heroic effort for sure! Did you not consider answering this question: mattermodeling.stackexchange.com/q/1676/5 ? You could have got a 50-point bounty! $\endgroup$ Commented Aug 4, 2020 at 21:36
  • $\begingroup$ @NikeDattani I am not too knowledgeable in the Ehrenfest dynamics. I have studied and done some basic exercises, but I am not a professional or even ametuer user of Ehrenfest dynamics. I hope the contrasts between BOMD and Ehrenfest dynamics makes it a little bit easier to understand the basic needs and speciality of the theory. Plus I didn't see it earlier enough! $\endgroup$
    – mykd
    Commented Aug 4, 2020 at 21:42
  • 1
    $\begingroup$ No problem. You can always see the active bounties available here: mattermodeling.stackexchange.com/?tab=bounties and you can also see all previous bounties here: chat.stackexchange.com/rooms/109905/bounty-land $\endgroup$ Commented Aug 4, 2020 at 21:44

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .