Path integral molecular dynamics

Path integral molecular dynamics (PIMD) is a method of incorporating quantum mechanics into molecular dynamics simulations using Feynman path integrals. Such simulations are particularly useful for studying nuclear quantum effects in light atoms and molecules such as hydrogen, helium, neon and argon, as well as in quantum rotators such as methane and hydrogen bonded systems such as water. In PIMD, one uses the Born–Oppenheimer approximation to separate the wavefunction into a nuclear part and an electronic part. The nuclei are treated quantum mechanically by mapping each quantum nucleus onto a classical system of several fictitious particles connected by springs (harmonic potentials) governed by an effective Hamiltonian, which is derived from Feynman's path integral. The resulting classical system, although complex, can be solved relatively quickly. There are now a number of commonly used condensed matter computer simulation techniques that make use of the path integral formulation including centroid molecular dynamics (CMD), and ring polymer molecular dynamics (RPMD). The same techniques are also used in path integral Monte Carlo (PIMC).

Theory behind the method

In the path integral formulation the canonical partition function (in one dimension) is written as

$Q(\beta ,V)=\int {\mathrm {d} }x_{1}\int _{x_{1}}^{x_{1}}Dx(\tau )e^{-S[x(\tau )]}$

where $S[x(\tau )]$ is the Euclidean action, given by

$S[x(\tau )]=\int _{0}^{\beta \hbar }d\tau H(x(\tau ))$

where $x(\tau )$ is the path in time $\tau$ and $H$ is the Hamiltonian. One can identify the inverse temperature, $\beta$ with an imaginary time $it/\hbar$, (a technique known as a Wick rotation) thus connecting this path integral to standard path integrals appearing elsewhere in physics. Trotter's formula is given by:

$e^{-\beta (H_{0}+V)}=\lim _{P\rightarrow \infty }(e^{-\beta H_{0}/P}e^{-\beta V/P})^{P}$

Where $P$ is called the Trotter number. In the limit where $P\rightarrow \infty$ this equation becomes exact. In the case where $P=1$ one obtains classical dynamics. Path integral molecular dynamics exploits this formula by choosing a finite $P$, which usually ranges between 5 and 100.

Application of this formula leads to:

$Q_{P}=\left({\frac {mP}{2\pi \beta \hbar ^{2}}}\right)^{P/2}\int \cdots \int {\mathrm {d} }x_{1}\cdots {\mathrm {d} }x_{P}e^{-\beta \Phi _{P}(x_{1}\cdots x_{P};\beta )}$

where the Euclidean time is discretized in units of

$\varepsilon ={\frac {\beta \hbar }{P}},P\in {\mathbb {Z} }$

$x_{t}=x(t\beta \hbar /P)$

$x_{P+1}=x_{1}$

and

$\Phi _{P}(x_{1}\cdots x_{P};\beta )={\frac {mP}{2\beta ^{2}\hbar ^{2}}}\sum _{t=1}^{P}(x_{t}-x_{t+1})^{2}+{\frac {1}{P}}\sum _{t=1}^{P}V(x_{t}).$

It has long been recognized that there is an isomorphism between this discretized quantum mechanical description, and the classical statistical mechanics of polyatomic fluids, in particular flexible ring molecules, due to the periodic boundary conditions in imaginary time. It can be seen from the first term of the above equation that each particle $x_{t}$ interacts with is neighbors $x_{t-1}$ and $x_{t+1}$ via a harmonic spring. The second term provides the internal potential energy.

In three dimensions one has partition function in terms of the density operator

$Q(\beta ,V)=\operatorname {Tr} ({\hat {\rho }}(\beta ))=\operatorname {Tr} \lbrace \exp \left[-\beta {\hat {H}}\right]\rbrace$

which, thanks to the Trotter's formula, allows one to tease out $\exp \left[-\beta (U_{\mathrm {spring} }+U_{\mathrm {internal} })\right]$, where

$U_{\mathrm {spring} }={\frac {mP}{2\beta ^{2}\hbar ^{2}}}\sum _{t=1}^{P}|\mathbf {r} _{t}-\mathbf {r} _{t+1}|^{2}$

and

$U_{\mathrm {internal} }={\frac {1}{P}}\sum _{t=1}^{P}V(\mathbf {r} _{t})$

The internal energy is given by

$\langle U\rangle ={\frac {3NP}{2\beta }}-\langle U_{\mathrm {spring} }\rangle +\langle U_{\mathrm {internal} }\rangle$

The average kinetic energy is known as the primitive estimator, i.e.

$\langle K_{P}\rangle ={\frac {3NP}{2\beta }}-\langle U_{\mathrm {spring} }\rangle$

Rotational degrees of freedom

In the case of systems having $d$ rotational degrees of freedom the Hamiltonian can be written in the form:

${\hat {H}}={\hat {T}}^{\mathrm {translational} }+{\hat {T}}^{\mathrm {rotational} }+{\hat {V}}$

where the rotational part of the kinetic energy operator is given by:

$T^{\mathrm {rotational} }=\sum _{i=1}^{d^{\mathrm {rotational} }}{\frac {{\hat {L}}_{i}^{2}}{2\Theta _{ii}}}$

where ${\hat {L}}_{i}$ are the components of the angular momentum operator, and $\Theta _{ii}$ are the moments of inertia.

Combination with other simulation techniques

Applications

The technique has been used to calculate time correlation functions.

References

Portions of this page come come from page Path Integral Formulation on the SklogWiki, licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License.

Further reading

• Feynman, R. P. (1972). "Chapter 3". Statistical Mechanics. Reading, Massachusetts: Benjamin. ISBN 0-201-36076-4.

External links

• "Density matrices and path integrals" (computer code). SMAC-wiki.
• John Shumway; Matthew Gilbert (2008). "Path Integral Monte Carlo Simulation".{{cite web}}: CS1 maint: multiple names: authors list (link)