„Deriving“ Molecular Dynamics

We start from the Schrödinger equation for a system of $M$ nuclei/ions:

$i\hbar {\frac {\partial }{\partial t}}{{\Lambda }_{n}}\left(\left\{{{\mathbf {R} }_{I}}\right\},t\right)=\left[-\sum \limits _{I}{{\frac {{\hbar }^{2}}{2{{M}_{I}}}}{{\Delta }_{{\mathbf {R} }_{I}}}+{{E}_{n}}\left(\left\{{{\mathbf {R} }_{I}}\right\}\right)}\right]\cdot {{\Lambda }_{n}}\left(\left\{{{\mathbf {R} }_{I}}\right\},t\right)$

We write the complex wave-function in terms of a real-valued amplitude and phase

${{\Lambda }_{n}}\left(\left\{{{\mathbf {R} }_{I}}\right\},t\right)={{A}_{n}}\left(\left\{{{\mathbf {R} }_{I}}\right\},t\right)\cdot \exp \left[{\frac {i}{\hbar }}{{S}_{n}}\left(\left\{{{\mathbf {R} }_{I}}\right\},t\right)\right]$

Here ${{A}^{2}}={{\rho }_{M}}$ is the $M$ -particle density function, and the index is neglected as $n$ is kept fixed from now on. $S$ is called the ‘action’ and defines the classical momentum associated with particle $I$ :

${\begin{aligned}\left\langle {{\hat {P}}_{I}}\right\rangle \left(t\right)=&\int _{V\times V\times ...\times V}{\Lambda _{n}^{*}\left(\left\{{{\mathbf {R} }_{I}}\right\},t\right)}{\frac {\hbar }{i}}{{\nabla }_{{\mathbf {R} }_{I}}}{{\Lambda }_{n}}\left(\left\{{{\mathbf {R} }_{I}}\right\},t\right){{d}^{3M}}R\\=&\int _{V\times V\times ...\times V}{{{\rho }_{M}}\left(\left\{{{\mathbf {R} }_{I}}\right\}\right){{\nabla }_{{\mathbf {R} }_{I}}}S\left(\left\{{{\mathbf {R} }_{I}}\right\}\right)}{{d}^{3M}}R+\underbrace {\int _{V\times V\times ...\times V}{A\left(\left\{{{\mathbf {R} }_{I}}\right\}\right){\frac {\hbar }{i}}{{\nabla }_{{\mathbf {R} }_{I}}}A\left(\left\{{{\mathbf {R} }_{I}}\right\}\right)}{{d}^{3M}}R} _{0}\\=&\int _{V\times V\times ...\times V}{{{\rho }_{M}}\left(\left\{{{\mathbf {R} }_{I}}\right\}\right){{P}_{I}}\left(\left\{{{\mathbf {R} }_{I}}\right\}\right)}{{d}^{3M}}R\end{aligned}}$

The first line is the momentum of the particle indexed by $I$ in quantum mechanics. Substituting into ${{\Lambda }_{n}}\left(\left\{{{\mathbf {R} }_{I}}\right\},t\right)$ and evaluating the integral gives the second line. The second term gives 0, because $A$ decays fast enough, i.e. $\int _{V}{A\nabla A}{{d}^{3}}R=\int _{V}{{\frac {1}{2}}\nabla {{A}^{2}}}{{d}^{3}}R={\frac {1}{2}}\int _{\partial V}{{{A}^{2}}dF}=0$ .

One can get the equations of motion for the amplitude and phase by solving the Schrödinger equation (which is an equation for complex quantities) for its real and imaginary parts.

${\begin{aligned}&{\frac {\partial S}{\partial t}}+\sum \limits _{I}{{\frac {1}{2{{M}_{I}}}}{{\left({{\nabla }_{{\mathbf {R} }_{I}}}S\right)}^{2}}+E\left(\left\{{{\mathbf {R} }_{I}}\right\}\right)}=\sum \limits _{I}{\frac {{\hbar }^{2}}{2{{M}_{I}}}}{\frac {{{\Delta }_{{\mathbf {R} }_{I}}}A}{A}}\\&{\frac {\partial A}{\partial t}}+\sum \limits _{I}{{\frac {1}{2{{M}_{I}}}}\left[2\left({{\nabla }_{{\mathbf {R} }_{I}}}A\right)\left({{\nabla }_{{\mathbf {R} }_{I}}}S\right)+A{{\Delta }_{{\mathbf {R} }_{I}}}S\right]}=0\\\end{aligned}}$

By multiplying the second line by $2A$ , one can get a continuity equation (cf. fluid dynamics) for the square amplitude. First move $A$ into the derivationː

$2A\cdot {\frac {\partial A}{\partial t}}={\frac {\partial {{A}^{2}}}{\partial t}}\quad A\cdot 2\left({{\nabla }_{{\mathbf {R} }_{I}}}A\right)\left({{\nabla }_{{\mathbf {R} }_{I}}}S\right)=\left({{\nabla }_{{\mathbf {R} }_{I}}}{{A}^{2}}\right)\left({{\nabla }_{{\mathbf {R} }_{I}}}S\right)$ , then let's notice that

$\left({{\nabla }_{{\mathbf {R} }_{I}}}{{A}^{2}}\right)\left({{\nabla }_{{\mathbf {R} }_{I}}}S\right)+{{A}^{2}}{{\Delta }_{{\mathbf {R} }_{I}}}S={{\nabla }_{{\mathbf {R} }_{I}}}\left({{A}^{2}}\cdot {{\nabla }_{{\mathbf {R} }_{I}}}S\right)$ , then introducing the quantity ${{\mathbf {J} }_{I}}={{A}^{2}}\cdot {{\nabla }_{{\mathbf {R} }_{I}}}S/{{M}_{I}}$ , one can get

${\frac {\partial {{A}^{2}}}{\partial t}}+\sum \limits _{I}{{{\nabla }_{{\mathbf {R} }_{I}}}{{\mathbf {J} }_{I}}=0}$ .

What can we state for the first line? Only one term depends on the action quantum $\hbar$ . When can this term be neglected? Let's compare with second term on left-hand side:

${\frac {1}{2{{M}_{I}}}}{{\left({{\nabla }_{{\mathbf {R} }_{I}}}S\right)}^{2}}\gg {\frac {{\hbar }^{2}}{2{{M}_{I}}}}{\frac {{{\Delta }_{{\mathbf {R} }_{I}}}A}{A}}\quad {\text{if}}\quad {{C}_{x}}\gg {\frac {\hbar }{{P}_{I}}}={{\lambda }_{I}}$

The ‘quantum mechanical’ right-hand term can be neglected if the spatial scale of the description is large as compared to the de-Broglie wavelength of the particles. For nuclei with energies of the order of E this is fulfilled if

${{\left({\frac {m}{M}}\right)}^{1/2}}{{C}_{r}}<{{10}^{-12}}m\ll {{C}_{x}}$

In this case the equation for S reduces to the Hamilton-Jacobi equation of classical mechanics

${\frac {\partial S}{\partial t}}+H\left(\left\{{{\mathbf {R} }_{I}}\right\},\left\{{{\nabla }_{{\mathbf {R} }_{I}}}S\right\}\right)=0$

Here, $H$ is the classical Hamilton function:

$H\left(\left\{{{\mathbf {R} }_{I}}\right\},\left\{{{\mathbf {P} }_{I}}\right\}\right)=E\left(\left\{{{\mathbf {R} }_{I}}\right\}\right)+\sum \limits _{I}{\frac {P_{I}^{2}}{2{{M}_{I}}}}$

The Hamilton-Jacobi partial differential equation is equivalent to 3M ordinary differential equations for coordinates and momenta (see classical mechanics texts).

${\begin{matrix}{\frac {d{{q}_{i}}}{dt}}={\frac {\partial H}{\partial {{p}_{i}}}}\quad {\frac {d{{p}_{i}}}{dt}}=-{\frac {\partial H}{\partial {{q}_{i}}}}&{\begin{matrix}{{q}_{i}}\in \left\{{{R}_{1,x}},{{R}_{1,y}},{{R}_{1,z}},...,{{R}_{M,x}},{{R}_{M,y}},{{R}_{M,z}}\right\}\\{{p}_{i}}\in \left\{{{P}_{1,x}},{{P}_{1,y}},{{P}_{1,z}},...,{{P}_{M,x}},{{P}_{M,y}},{{P}_{M,z}}\right\}\\\end{matrix}}\\\end{matrix}}$

These are directly equivalent to Newton’s equations of motionː

${\frac {d{{\mathbf {R} }_{I}}}{dt}}={\frac {{\mathbf {P} }_{I}}{{M}_{I}}}\quad {\frac {d{{\mathbf {P} }_{I}}}{dt}}=-{{\nabla }_{I}}E$

Intermediate Summary:

Dynamics of electrons requires quantum mechanical description, but may evaluated with the nuclei/ions assumed stationary
This leads to an energy functional which acts as a potential energy surface on which nuclei/ions move according to deterministic dynamics.
The energy functional depends on all ion coordinates (there is in general no such thing as e.g. ‘pair interactions’). Its evaluation requires, in principle, to solve the Schrödinger equation for all electrons in a high-dimensional space. This is usually done using efficient single electron approximations (density-functional theory).
Determining the energy surface for all possible nuclear configurations is technically impossible due to dimensionality. Thus in ab-initio MD one evaluates only the configurations that are visited along a trajectory .
Alternatively, one may of course use empirical potentials which may be considered terms in a multi-particle expansion of the true energy surface.