Path integral molecular dynamics

Path integral molecular dynamics (PIMD) is a method of incorporating quantum mechanics into molecular dynamics simulations using Feynman path integrals. 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), ring polymer molecular dynamics (RPMD), and the Feynman–Kleinert quasi-classical Wigner (FK–QCW) method. The same techniques are also used in path integral Monte Carlo (PIMC).

There are two ways to calculate the dynamics calculations of PIMD. The first one is the non-Hamiltonian phase space analysis theory, which has been updated to create an "extended system" of isokinetic equations of motion which overcomes the properties of a system that created issues within the community. The second way is by using Nosé–Hoover chain, which is a chain of variables instead of a single thermostat of variable.

Ring-polymer representation of the partition function

Consider a single distinguishable particle of mass m {\displaystyle m} moving in one dimension, with Hamiltonian

H ^ = T ^ + V ^ = p ^ 2 2 m + V ( q ^ ) . {\displaystyle {\hat {H}}={\hat {T}}+{\hat {V}}={\frac {{\hat {p}}^{2}}{2m}}+V({\hat {q}}).}

Its canonical partition function is

Q = Tr ⁡ ( e − β H ^ ) , β = 1 k B T . {\displaystyle Q=\operatorname {Tr} \left(e^{-\beta {\hat {H}}}\right),\qquad \beta ={\frac {1}{k_{\mathrm {B} }T}}.}

The Boltzmann operator can be divided into n {\displaystyle n} imaginary-time slices,

e − β H ^ = ( e − β n H ^ ) n , β n = β n . {\displaystyle e^{-\beta {\hat {H}}}=\left(e^{-\beta _{n}{\hat {H}}}\right)^{n},\qquad \beta _{n}={\frac {\beta }{n}}.}

Inserting complete sets of position eigenstates between the factors gives

Q = ∫ d q 1 ⋯ ∫ d q n ∏ j = 1 n ⟨ q j | e − β n H ^ | q j + 1 ⟩ , q n + 1 = q 1 . {\displaystyle Q=\int \mathrm {d} q_{1}\cdots \int \mathrm {d} q_{n}\prod _{j=1}^{n}\left\langle q_{j}\left|e^{-\beta _{n}{\hat {H}}}\right|q_{j+1}\right\rangle ,\qquad q_{n+1}=q_{1}.}

The condition q n + 1 = q 1 {\displaystyle q_{n+1}=q_{1}} follows from the trace and makes the discretized path cyclic. For sufficiently large n {\displaystyle n}, the symmetric Suzuki–Trotter factorization gives

e − β n ( T ^ + V ^ ) = e − β n V ^ / 2 e − β n T ^ e − β n V ^ / 2 + O ( β n 3 ) . {\displaystyle e^{-\beta _{n}({\hat {T}}+{\hat {V}})}=e^{-\beta _{n}{\hat {V}}/2}e^{-\beta _{n}{\hat {T}}}e^{-\beta _{n}{\hat {V}}/2}+{\mathcal {O}}(\beta _{n}^{3}).}

The free-particle imaginary-time propagator is obtained by inserting momentum eigenstates:

⟨ q | e − β n T ^ | q ′ ⟩ = 1 2 π ℏ ∫ − ∞ ∞ d p exp ⁡ [ − β n p 2 2 m + i p ( q − q ′ ) ℏ ] = m 2 π β n ℏ 2 exp ⁡ [ − m ( q − q ′ ) 2 2 β n ℏ 2 ] . {\displaystyle {\begin{aligned}\left\langle q\left|e^{-\beta _{n}{\hat {T}}}\right|q'\right\rangle &={\frac {1}{2\pi \hbar }}\int _{-\infty }^{\infty }\mathrm {d} p\,\exp \left[-{\frac {\beta _{n}p^{2}}{2m}}+{\frac {ip(q-q')}{\hbar }}\right]\\[4pt]&={\sqrt {\frac {m}{2\pi \beta _{n}\hbar ^{2}}}}\,\exp \left[-{\frac {m(q-q')^{2}}{2\beta _{n}\hbar ^{2}}}\right].\end{aligned}}}

Defining the ring-polymer frequency

ω n = 1 β n ℏ = n β ℏ , {\displaystyle \omega _{n}={\frac {1}{\beta _{n}\hbar }}={\frac {n}{\beta \hbar }},}

the short-time density matrix becomes

⟨ q | e − β n H ^ | q ′ ⟩ ≃ m 2 π β n ℏ 2 exp ⁡ { − β n [ 1 2 m ω n 2 ( q − q ′ ) 2 + V ( q ) + V ( q ′ ) 2 ] } . {\displaystyle \left\langle q\left|e^{-\beta _{n}{\hat {H}}}\right|q'\right\rangle \simeq {\sqrt {\frac {m}{2\pi \beta _{n}\hbar ^{2}}}}\,\exp \left\{-\beta _{n}\left[{\frac {1}{2}}m\omega _{n}^{2}(q-q')^{2}+{\frac {V(q)+V(q')}{2}}\right]\right\}.}

A set of auxiliary momenta can be introduced using the Gaussian identity

m 2 π β n ℏ 2 e − β n m ω n 2 ( q − q ′ ) 2 / 2 = 1 2 π ℏ ∫ − ∞ ∞ d p e − β n [ p 2 / ( 2 m ) + m ω n 2 ( q − q ′ ) 2 / 2 ] . {\displaystyle {\sqrt {\frac {m}{2\pi \beta _{n}\hbar ^{2}}}}\,e^{-\beta _{n}m\omega _{n}^{2}(q-q')^{2}/2}={\frac {1}{2\pi \hbar }}\int _{-\infty }^{\infty }\mathrm {d} p\,e^{-\beta _{n}[p^{2}/(2m)+m\omega _{n}^{2}(q-q')^{2}/2]}.}

Applying this identity to every imaginary-time slice gives

Q = lim n → ∞ Q n , {\displaystyle Q=\lim _{n\to \infty }Q_{n},}

where

Q n = 1 ( 2 π ℏ ) n ∫ d p ∫ d q e − β n H n ( p , q ) {\displaystyle Q_{n}={\frac {1}{(2\pi \hbar )^{n}}}\int \mathrm {d} {\boldsymbol {p}}\,\int \mathrm {d} {\boldsymbol {q}}\,e^{-\beta _{n}H_{n}({\boldsymbol {p}},{\boldsymbol {q}})}}

and

H n ( p , q ) = ∑ j = 1 n [ p j 2 2 m + 1 2 m ω n 2 ( q j − q j + 1 ) 2 + V ( q j ) ] , q n + 1 = q 1 . {\displaystyle H_{n}({\boldsymbol {p}},{\boldsymbol {q}})=\sum _{j=1}^{n}\left[{\frac {p_{j}^{2}}{2m}}+{\frac {1}{2}}m\omega _{n}^{2}(q_{j}-q_{j+1})^{2}+V(q_{j})\right],\qquad q_{n+1}=q_{1}.}

Thus the quantum canonical partition function is mapped onto the classical partition function of a cyclic polymer containing n {\displaystyle n} beads. Adjacent beads are connected by harmonic springs of frequency ω n {\displaystyle \omega _{n}}, and each bead experiences the physical potential V {\displaystyle V}. The variables p j {\displaystyle p_{j}} are auxiliary sampling momenta and should not be identified with measurements of the quantum momentum. Path integral molecular dynamics samples this ring-polymer distribution using classical molecular-dynamics trajectories.

For the symmetric factorization above, the leading finite-n {\displaystyle n} discretization error scales as O ( n − 2 ) {\displaystyle {\mathcal {O}}(n^{-2})}, subject to the usual regularity conditions on the potential.

Equilibrium estimators and molecular-dynamics sampling

The ring-polymer construction establishes a classical isomorphism for equilibrium statistical mechanics: the quantum system at inverse temperature β {\displaystyle \beta } is represented by a classical ring polymer whose phase-space weight contains β n = β / n {\displaystyle \beta _{n}=\beta /n}. Although the auxiliary ring polymer is sampled with the factor e − β n H n {\displaystyle e^{-\beta _{n}H_{n}}}, the physical temperature remains 1 / ( k B β ) {\displaystyle 1/(k_{\mathrm {B} }\beta )}; its dependence is contained in both β n {\displaystyle \beta _{n}} and the spring frequency ω n = n / ( β ℏ ) {\displaystyle \omega _{n}=n/(\beta \hbar )}.

For an observable represented by an operator A ^ {\displaystyle {\hat {A}}}, the canonical quantum expectation value is

⟨ A ^ ⟩ = 1 Q Tr ⁡ ( e − β H ^ A ^ ) . {\displaystyle \langle {\hat {A}}\rangle ={\frac {1}{Q}}\operatorname {Tr} \left(e^{-\beta {\hat {H}}}{\hat {A}}\right).}

If the observable depends only on position, A ^ = A ( q ^ ) {\displaystyle {\hat {A}}=A({\hat {q}})}, its ring-polymer representation is

⟨ A ^ ⟩ = lim n → ∞ 1 ( 2 π ℏ ) n Q n ∫ d p ∫ d q e − β n H n ( p , q ) A n ( q ) , {\displaystyle \langle {\hat {A}}\rangle =\lim _{n\rightarrow \infty }{\frac {1}{(2\pi \hbar )^{n}Q_{n}}}\int \mathrm {d} {\boldsymbol {p}}\int \mathrm {d} {\boldsymbol {q}}\,e^{-\beta _{n}H_{n}({\boldsymbol {p}},{\boldsymbol {q}})}A_{n}({\boldsymbol {q}}),}

where

A n ( q ) = 1 n ∑ j = 1 n A ( q j ) {\displaystyle A_{n}({\boldsymbol {q}})={\frac {1}{n}}\sum _{j=1}^{n}A(q_{j})}

is the bead-averaged estimator. Thus the observable is first averaged over the beads and then over the canonical distribution of the ring polymer.

Path integral molecular dynamics

Path integral molecular dynamics uses fictitious classical dynamics to sample the ring-polymer canonical distribution. If the sampling dynamics are ergodic and preserve this distribution, the phase-space ensemble average may be evaluated as a long-time average along a trajectory.

Writing the ring-polymer Hamiltonian as

H n ( p , q ) = ∑ j = 1 n p j 2 2 m + U n ( q ) , {\displaystyle H_{n}({\boldsymbol {p}},{\boldsymbol {q}})=\sum _{j=1}^{n}{\frac {p_{j}^{2}}{2m}}+U_{n}({\boldsymbol {q}}),}

where

U n ( q ) = ∑ j = 1 n [ 1 2 m ω n 2 ( q j − q j + 1 ) 2 + V ( q j ) ] , q n + 1 = q 1 , {\displaystyle U_{n}({\boldsymbol {q}})=\sum _{j=1}^{n}\left[{\frac {1}{2}}m\omega _{n}^{2}(q_{j}-q_{j+1})^{2}+V(q_{j})\right],\qquad q_{n+1}=q_{1},}

the corresponding Hamilton equations are

q ˙ = ∂ H n ∂ p = p m , p ˙ = − ∂ H n ∂ q = − ∂ U n ∂ q . {\displaystyle {\dot {\boldsymbol {q}}}={\frac {\partial H_{n}}{\partial {\boldsymbol {p}}}}={\frac {\boldsymbol {p}}{m}},\qquad {\dot {\boldsymbol {p}}}=-{\frac {\partial H_{n}}{\partial {\boldsymbol {q}}}}=-{\frac {\partial U_{n}}{\partial {\boldsymbol {q}}}}.}

Unthermostatted Hamiltonian dynamics conserves H n {\displaystyle H_{n}} and therefore samples a microcanonical rather than a canonical distribution. In addition, energy exchange between weakly coupled ring-polymer modes can be inefficient, and unthermostatted PIMD can be nonergodic for some systems.

Normal-mode representation

Efficient PIMD algorithms frequently transform the free ring polymer into its normal modes. For an even number of beads, relabeled as j = 0 , … , n − 1 {\displaystyle j=0,\ldots ,n-1}, a real orthogonal transformation is

P k = ∑ j = 0 n − 1 C j k p j , Q k = ∑ j = 0 n − 1 C j k q j , {\displaystyle P_{k}=\sum _{j=0}^{n-1}C_{jk}p_{j},\qquad Q_{k}=\sum _{j=0}^{n-1}C_{jk}q_{j},}

where

C j k = { 1 / n , k = 0 , 2 / n cos ⁡ ( 2 π j k / n ) , 1 ≤ k ≤ n / 2 − 1 , 1 / n ( − 1 ) j , k = n / 2 , 2 / n sin ⁡ ( 2 π j k / n ) , n / 2 + 1 ≤ k ≤ n − 1. {\displaystyle C_{jk}={\begin{cases}{\sqrt {1/n}},&k=0,\\[4pt]{\sqrt {2/n}}\cos(2\pi jk/n),&1\leq k\leq n/2-1,\\[4pt]{\sqrt {1/n}}(-1)^{j},&k=n/2,\\[4pt]{\sqrt {2/n}}\sin(2\pi jk/n),&n/2+1\leq k\leq n-1.\end{cases}}}

In these coordinates, the free ring-polymer Hamiltonian is diagonal:

H n ( 0 ) ( P , Q ) = ∑ k = 0 n − 1 [ P k 2 2 m + 1 2 m ω k 2 Q k 2 ] , {\displaystyle H_{n}^{(0)}({\boldsymbol {P}},{\boldsymbol {Q}})=\sum _{k=0}^{n-1}\left[{\frac {P_{k}^{2}}{2m}}+{\frac {1}{2}}m\omega _{k}^{2}Q_{k}^{2}\right],}

with normal-mode frequencies

ω k = 2 ω n sin ⁡ ( k π n ) . {\displaystyle \omega _{k}=2\omega _{n}\sin \left({\frac {k\pi }{n}}\right).}

The mode k = 0 {\displaystyle k=0} has zero spring frequency and corresponds to the ring-polymer centroid,

Q 0 = 1 n ∑ j = 0 n − 1 q j . {\displaystyle Q_{0}={\frac {1}{\sqrt {n}}}\sum _{j=0}^{n-1}q_{j}.}

PILE thermostat

The path integral Langevin equation (PILE) thermostat applies a frequency-dependent Langevin thermostat to the ring-polymer normal modes. Including the physical potential, the continuous-time equations can be written as

Q ˙ k = P k m , {\displaystyle {\dot {Q}}_{k}={\frac {P_{k}}{m}},}

P ˙ k = − m ω k 2 Q k − ∂ V n ∂ Q k − γ k P k + 2 m γ k β n ξ k ( t ) , {\displaystyle {\dot {P}}_{k}=-m\omega _{k}^{2}Q_{k}-{\frac {\partial V_{n}}{\partial Q_{k}}}-\gamma _{k}P_{k}+{\sqrt {\frac {2m\gamma _{k}}{\beta _{n}}}}\,\xi _{k}(t),}

where

V n ( Q ) = ∑ j = 0 n − 1 V ( q j ( Q ) ) {\displaystyle V_{n}({\boldsymbol {Q}})=\sum _{j=0}^{n-1}V(q_{j}({\boldsymbol {Q}}))}

and the independent Gaussian white noises satisfy

⟨ ξ k ( t ) ⟩ = 0 , ⟨ ξ k ( t ) ξ k ′ ( t ′ ) ⟩ = δ k k ′ δ ( t − t ′ ) . {\displaystyle \langle \xi _{k}(t)\rangle =0,\qquad \langle \xi _{k}(t)\xi _{k'}(t')\rangle =\delta _{kk'}\delta (t-t').}

For the internal modes, the PILE choice that minimizes the autocorrelation time of the free ring-polymer energy is

γ k = 2 ω k , k > 0. {\displaystyle \gamma _{k}=2\omega _{k},\qquad k>0.}

Since ω 0 = 0 {\displaystyle \omega _{0}=0}, this prescription does not thermostat the centroid. In the local version, PILE-L, the centroid is assigned an independent friction coefficient

γ 0 = 1 τ 0 , {\displaystyle \gamma _{0}={\frac {1}{\tau _{0}}},}

where τ 0 {\displaystyle \tau _{0}} is a user-selected thermostat time scale. In the global version, PILE-G, the centroid is instead coupled to a global stochastic velocity-rescaling thermostat. The stochastic velocity-rescaling method generates the canonical distribution by rescaling all selected momenta with a common random factor.

Relation to real-time dynamics

The ring-polymer isomorphism is an equilibrium statistical-mechanical relation. The fictitious trajectories used in PIMD are therefore sampling trajectories and are not, in general, the exact real-time quantum dynamics of the original system. Configurational equilibrium averages are unchanged by a consistent choice of positive fictitious masses for the ring-polymer beads.

Ring-polymer molecular dynamics (RPMD) makes an additional dynamical approximation. It assigns the physical particle masses to the ring-polymer beads and propagates the ring-polymer Hamiltonian dynamics without a thermostat in order to approximate Kubo-transformed quantum time-correlation functions.

Combination with other simulation techniques

The simulations done my PIMD can broadly characterize the biomolecular systems, covering the entire structure and organization of the membrane, including the permeability, protein-lipid interactions, along with "lipid-drug interactions, protein–ligand interactions, and protein structure and dynamics."

Applications

PIMD is "widely used to describe nuclear quantum effects in chemistry and physics".

Path Integral Molecular Dynamics can be applied to polymer physics, both field theories, quantum and not, string theory, stochastic dynamics, quantum mechanics, and quantum gravity. PIMD can also be used to calculate time correlation functions

External links

. SMAC-wiki. Archived from (computer code) on May 1, 2016.
John Shumway; Matthew Gilbert (2008). . doi:.{{cite journal}}:Cite journal requires |journal= (help)CS1 maint: multiple names: authors list (link)

Path integral molecular dynamics

Ring-polymer representation of the partition function

Equilibrium estimators and molecular-dynamics sampling

Path integral molecular dynamics

Normal-mode representation

PILE thermostat

Relation to real-time dynamics

Combination with other simulation techniques

Applications

Further reading

External links