11.2.3: Path integral molecular dynamics (optional reading)
 Page ID
 5279
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{\!\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)
( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\ #1 \}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\id}{\mathrm{id}}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\kernel}{\mathrm{null}\,}\)
\( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\)
\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\)
\( \newcommand{\norm}[1]{\ #1 \}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)
\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)
\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)
\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vectorC}[1]{\textbf{#1}} \)
\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)
\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)
\( \newcommand{\vectE}[1]{\overset{\!\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{\!\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
Consider once again the path integral expression for the onedimensional canonical partition function (for a finite but large value of \(P\)):
\[Q(\beta) =\left({mP \over 2\pi\beta\hbar^2}\right)^{P/2}\int dx_1 \cdots dx_P exp \left [  \beta \sum _{i=1}^P \left ({1 \over 2}m\omega_P^2(x_{i+1}x_i)^2+ {1 \over P}U(x_i)\right)\right] \nonumber \]  (1) 
(the condition \( {x_{P+1} = x_1} \) is understood). Recall that, according to the classical isomorphism, the path integral expression for the canonical partition function is isomorphic to the classical configuration integral for a certain \(P\)particle system. We can carry this analogy one step further by introducing into the above expression a set of \(P\) momentum integrations:
\[ Q(\beta) =\int dp_1\cdots dp_P\;dx_1\cdots dx_P exp\left [  \beta \sum _{i=1}^P \left ( {P^2_i \over 2m'} + {1 \over 2}m\omega_P^2(x_{i+1}x_i)^2+ {1 \over P}U(x_i)\right)\right] \nonumber \]  (2) 
Note that these momentum integrations are completely uncoupled from the position integrations, and if we were to carry out these momentum integrations, we would reproduce Eq. (1) apart from trivial constants. Written in the form Eq. (2), however, the path integral looks exactly like a phase space integral for a \(P\)particle system. We know from our work in classical statistical mechanics that dynamical equations of motion can be constructed that will generate this partition function. In principle, one would start with the classical Hamiltonian
\[ H = \sum_{i=1}^P \left[{p_i^2 \over 2m'} + {1 \over 2}m\omega_P^2(x_{i+1}x_i)^2 + {1 \over P}U(x_i)\right] \nonumber \]
derive the corresponding classical equations of motion and then couple in thermostats. Such an approach has certainly been attempted with only limited success. The difficulty with this straightforward approach is that the more "quantum'' a system is, the large the paramester \(P\) must be chosen in order to converge the path integral. However, if \(P\) is large, the above Hamiltonian describes a system with extremely stiff nearestneighbor harmonic bonds interacting with a very weak potential \(U/P \). It is, therefore, almost impossible for the system to deviate far harmonic oscillator solutions and explore the entire available phase space. The use of thermostats can help this problem, however, it is also exacerbated by the fact that all the harmonic interactions are coupled, leading to a wide variety of time scales associated with the motion of each variable in the Hamiltonian. In order to separate out all these time scales, one must somehow diagonalize this harmonic interaction. One way to do this is to use normal mode variables, and this is a perfectly valid approach. However, we will explore another, simpler approach here. It involves the use of a variable transformation of the formed used in previous lectures to do the path integral for the freeparticle density matrix.
Consider a change of variables:
\( {u_1} \)  =  \( {x_1} \)  
\( {u_k} \)  =  \(x_k  \tilde {x}_k k = 2, \cdots , P \) 
where
\[ \tilde{x}_i = {(k1)x_{k+1} + x_1 \over k} \nonumber \]
The inverse of this transformation can be worked out in closed form:
\( {x_1}\)  =  \( {u_1} \)  
\( {x_k} \)  =  \( u_1 + \sum_{l=k}^P {(k1) \over (l1)}u_l \) 
and can also be expressed as a recursive inverse:
\( {x_1} \)  =  \( {u_1} \)  
\( {x_k}\)  =  \( u_1 + \sum_{l=k}^P {(k1) \over (l1)}u_l \) 
The term \(k = P \) here can be used to start the recursion. We have already seen that this transformation diagonalized the harmonic interaction. Thus, substituting the transformation into the path integral gives:
\[ Q(\beta) =\int dp_1\cdots dp_P\;du_1\cdots du_P exp\left [\beta \sum _{i=1}^P \left ( {P^2_i \over 2m'_i} + {1 \over 2} m_i \omega_P^2u_i^2+ {1 \over P}U(x_i(u_1,...,u_P))\right)\right] \nonumber \]
The parameters \( {m_i} \) are given by
\( {m_1} \)  =  \( {0} \)  
\( {m_i} \)  =  \( { {i \over i1}m } \) 
Note also that the momentum integrations have been changed slightly to involve a set of parameters \(m'_i\). Introducing these parameters, again, only changes the partition function by trivial constant factors. How these should be chosen will become clear later in the discussion. The notation \(x_i (u_1, \cdots , u_P) \) indicates that each variable \( {x_i} \) is a generally a function of all the new variables \( {u_1, \cdots , u_P }\).
A dynamics scheme can now be derived using as an effective Hamiltonian:
\[ H = \sum_{i=1}^P \left[{p_i^2 \over 2m_i'} + {1 \over 2}m_i\omega_P^2u_i^2 + {1 \over P}U(x_i(u_1,...,u_P))\right] \nonumber \]
which, when coupled to thermostats, yields a set of equations of motion
\( {m_i'\ddot{u}_i}\)  = 
\( m_i\omega_P^2 u_i  {1 \over P}{\partial U \over \partial u_i}\dot{\eta}_i\dot{u}_i \) 

\( Q\ddot{\eta}_i \)  =  \( m_i\dot{u}_i^2  {1 \over \beta} \)  (3) 
These equations have a conserved energy (which is not a Hamiltonian):
\[ H' = \sum_{i=1}^P\left[{1 \over 2}m_i'\dot{u}_i^2 +{1 \over 2}m_i \omega^2_P u^2_i + {1 \over P} U (x_i (u_1, \cdots , u_P)) +{1 \over 2}Q\dot{\eta}_i^2 + {1 \over \beta}\eta_i\right] \nonumber \]
Notice that each variable is given its own thermostat. This is done to produce maximum ergodicity in the trajectories. In fact, in practice, the chain thermostats you have used in the computer labs are employed. Notice also that the time scale of each variable is now clear. It is just determined by the parameters \(\{m_i \} \). Since the object of using such dynamical equations is not to produce real dynamics but to sample the phase space, we would really like each variable to move on the same time scale, so that there are no slow beads trailing behind the fast ones. This effect can be produced by choosing each parameter \( {m'_i}\) to be proportional to \( {m_i} : m'_i = cm_i \) Finally, the forces on the \(u\) variables can be determined easily from the chain rule and the recursive inverse given above. The result is
\( {1 \over P}{\partial U \over \partial u_1} \)  =  \( {1 \over P}\sum_{i=1}^P{\partial U \over \partial x_i} \)  
\( {1 \over P}{\partial U \over \partial u_i} \)  =  \( {1 \over P}\left[{(k2) \over (k1)}{\partial U \over \partial u_{k1} }+ {\partial U \over \partial x_k}\right] \) 
where the first \( ( {i = 1}) \) of these expressions starts the recursion in the second equation.
Later on, when we discuss applications of path integrals, we will see why a formulation such as this for evaluating path integrals is advantageous.