Franz J. Vesely > CompPhys Tutorial > Selected Applications > Stochastic Dynamics

 < >
 Part III: Ch. 7

7.2 Path Integral Monte Carlo (PIMC)

So far: only ground state, i.e. $T=0$. $\Longrightarrow$How to treat quantum systems at finite temperature?

Let a quantum system be in contact with a heat bath of temperature $kT > 0$; it is then in a mixed state made up of the energy eigenstates:

$\Psi=\sum \limits_{n} c_{n} \Psi_{n} \, , \;\;\;\; {\rm where} \;\; {\rm H} \Psi_{n} = E_{n} \, \Psi_{n}$     (7.17)

In classical statistical mechanics the probability density of states is given by the Boltzmann factor. In quantum mechanics the density matrix serves the same purpose:

$\begin{eqnarray} \rho(r ,r ';kT) & \equiv & \sum \limits_{n}\Psi_{n}^{*}(r )\, e^{\textstyle - {\rm H}/kT}\,\Psi_{n}(r' ) \\ &=& \sum \limits_{n}\Psi_{n}^{*}(r )\,e^{\textstyle -E_{n}/kT}\,\Psi_{n}(r' ) \end{eqnarray}$     (7.18)

The average of some observable $a(r)$ is given by

$\langle a \rangle =\int a(r )\, \rho(r ,r ;\beta)\, dr \, \left/ \int \rho(r ,r ;\beta)\, dr \equiv Sp[a \rho] \right/ Sp[\rho]$     (7.19)

where $\beta \equiv 1/kT$.

Problem: The explicit form of $\rho(r ,r ;\beta)$ is usually complex or unknown.
Solution: Reduce the above expression to one that contains only known and tractable density matrices.

Back to school: What are the explicit forms of the density matrix for very simple models?

Strategy of PIMC: Express the density matrix of any given system in terms of the free particle density
7.24. This can be done because:

$\begin{eqnarray} \rho(x,x';\beta) &=&\sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta {\rm H}} \Psi_{n}(x')= \\ &=&\sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta {\rm H}/2}\, e^{-\beta {\rm H}/2} \Psi_{n}(x')= \\ &=&\sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta {\rm H}/2} \int dx'' \delta(x'-x'')\, e^{-\beta {\rm H}/2} \Psi_{n}(x'')= \\ &=&\sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta {\rm H}/2} \int dx'' \sum_{m} \Psi_{m}(x') \, \Psi_{m}^{*}(x'')\,e^{-\beta {\rm H}/2} \Psi_{n}(x'')= \\ &=&\int dx'' \left[ \sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta {\rm H}/2} \, \Psi_{n}(x'')\right] \left[ \sum_{m}\Psi_{m}^{*}(x'') \, e^{-\beta {\rm H}/2} \Psi_{m}(x') \right] \end{eqnarray}$

The density matrix may therefore be written as

$\rho(x,x';\beta)=\int dx'' \rho(x,x''; \frac{\textstyle \beta}{2})\, \rho(x'',x';\frac{\textstyle \beta}{2})$     (7.29)

The integrand in the path integral on the r.h.s. consists of two density matrices at $\beta/2$, i.e. double the original temperature.

Now we iterate, writing $x_{0},x_{1},x_{2}\dots$ instead of $x,x',x'',\dots$:

The number $P$ of intermediate steps is called the Trotter number. If $P$ is large ($5-100$) the high temperature $PT$ will act to damp the effect of the potential $U(x)$:

with diagonal element
$\rho(x_{0},x_{0};\beta)= A^{P/2} \int \dots \int dx_{1} \dots dx_{P-1}\, e^{\textstyle -\beta (U_{int}+U_{ext})}$     (7.32)

where

$A \equiv \frac{\textstyle mP}{\textstyle 2\pi \beta \hbar^{2}}\,,\;\;\;\; U_{ext} \equiv \sum_{p=0}^{P-1} U(x_{p})/P$     (7.33)

and

$U_{int} \equiv \frac{\textstyle A \pi}{\textstyle \beta} \sum_{p=0}^{P-1} \left(x_{p}-x_{p+1} \right)^{2}$     (7.34)

Three dimensions:

$\rho(r_{0},r_{0};\beta)= A^{3P/2} \int \dots \int dr_{1} \dots dr_{P-1}\, e^{\textstyle -\beta (U_{int}+U_{ext})}$     (7.35)

with the same $A$ as above, and

$U_{ext} \equiv \sum_{p=0}^{P-1} U(r_{p})/P\,,\;\;\;\; U_{int} \equiv \frac{\textstyle A \pi}{\textstyle \beta} \sum_{p=0}^{P-1} \left|r_{p}-r_{p+1} \right|^{2}$     (7.36)

Visualization / Classical isomorphism: The expression 7.35 looks like a classical Boltzmann factor for a ring polymer of $P$ links [CHANDLER 81]:
Each link is under the influence of an external potential

$u_{ext}(x_{p})=U(x_{p})/P$     (7.37)

while successive links are coupled by a harmonic bond potential

$u_{int}(x_{p},x_{p+1})=\frac{\textstyle A \pi}{\textstyle \beta}\, \left(x_{p}-x_{p+1} \right)^{2}$     (7.38)

with $x_{P}=x_{0}$.

Figure 7.2: Classical isomorphism for one particle

Therefore: To solve the quantum problem, apply Monte Carlo method to the classical ring polymer.

 Path integral Monte Carlo for one particle: A particle that is under the influence of an external potential $U(r)$ is represented by a ring polymer consisting of $P$ links. Let $r \equiv \{ r_{0},\dots r_{P-1}\}$ and the according potential energy $U_{pot}(r) \equiv U_{int}(r)+U_{ext}(r)$ be given, with $\begin{eqnarray} U_{int} & \equiv & \frac{\textstyle mP}{\textstyle 2 \beta^{2} \hbar^{2}}\, \sum_{p=0}^{P-1} |r_{p}-r_{p+1}|^{2} \\ U_{ext} & \equiv & \frac{\textstyle 1}{\textstyle P}\, \sum_{p=0}^{P-1}\,U(r_{p}) \end{eqnarray}$ Displace $r$ as a whole by $\Delta r$ (large); also, move each link $r_{p}$ by a small amount $\Delta r_{p}$; the new configuration is called $r'$. Compute $U_{pot}(r')$ and $\Delta U \equiv U_{pot}(r')- U_{pot}(r)$. Metropolis step: Draw $\xi$ from an equidistribution in $[0,1]$;    if $\Delta U \leq 0$, put $r=r'$;    if $\Delta U > 0$ and $\xi \leq e^{-\beta \, \Delta U}$, put $r=r'$ as well;    if $\Delta U > 0$ and $\xi > e^{-\beta \, \Delta U}$, let $r$ remain unchanged. Return to 1.

Figure 7.3: PIMC for one particle

Note:
- The oscillator strength $k=Pm/\beta^{2} \hbar^{2}$ increases with larger Trotter numbers $P$, requiring smaller MC steps
- The external potential $u_{ext}(r_{p})=U(r_{p})/P$ is damped by larger $P$, allowing for larger steps
$\Longrightarrow$ Move the entire ring polymer by a large random step, then displace the individual elements by a small amount.
Or:
Displace the center of mass of the chain by a wide random step, then construct the entire ring polymer anew, sampling the single element positions from the narrow multivariate Gauss distribution (see Chapter "Stochastics"):

$p(r_{0},\dots r_{P-1}) \propto exp \left[ -\frac{\textstyle mP}{\textstyle 2\beta \hbar^{2}} \sum_{p=0}^{P-1}(r_{p}-r_{p+1})^{2} \right]$     (7.39)

$N$ interacting quantum objects:
- Assume a pair potential $u(|r_{j}-r_{i}|)$.
- Each quantum particle is represented by a classical $P$-element ring chain.
- Let $r_{i,p}$ be the position of element $p\,$ in chain (=particle) $i\,$
- The diagonal element of the total density matrix is then

$\begin{eqnarray} \rho(r_{0},r_{0};\beta)&=& A^{3NP/2} \int \dots \int dr_{1,1} \dots dr_{N,P-1}\, exp \left[- A\pi \, \sum_{i=1}^{N} \sum_{p=0}^{P-1} \left( r_{i,p}- r_{i,p+1} \right)^{2} \right] \\ && \cdot exp \left[- \frac{\textstyle \beta}{\textstyle P}\sum_{i=1}^{N}\sum_{j>i}\sum_{p=0}^{P-1} u(|r_{j,p}-r_{i,p}|) \right] \end{eqnarray}$     (7.40)

with $r_{0}\equiv \{ r_{1,0} \dots r_{N,0} \}$. The pair potential acts only between respective links ($p$) of different chains.

EXERCISE: Write a PIMC program treating the case of one particle of mass $m$ in a two-dimensional oscillator potential $U(r)=kr^{2}/2$. Let the Trotter number vary between $2$ and $10$. Determine the positional probability $p(r)$ of the particle from the relative frequency of residence at $r$, averaged over all chain links. Noting that

$p(r) \equiv \rho(r,r;\beta)$     (7.41)

we would expect for the two-dimensional harmonic oscillator (with $\omega_{0}^{2}=k/m$)

$p(r)=2\pi r \left[ \frac{\textstyle A}{\textstyle \pi}\right] \, e^{-Ar^{2}}\,,\;\;\;\; {\rm where}\;\; A \equiv \frac{\textstyle m \omega_{0}^{2}}{\textstyle \hbar}\, tanh\frac{\textstyle \beta \hbar \omega_{0}}{2}$     (7.42)

(For convenience, put $\hbar =1$.) Draw several configurations of the ring polymer that occur in the course of the simulation.

Sample applications:
• Solvated electron in molten KCl [PARRINELLO 84]: localized or not? $\Longrightarrow$ Answer: clear localization.

• Solvation of electrons in simple fluids [COKER 87]
- electron in liquid helium: strongly localized
- electron in liquid xenon smeared out (see Figure)

Reason: the shell of He is rigid and difficult to polarize, repulsing an extra electron; the shell of Xe is easily polarizable, producing a long-ranged, locally flat dipole potential.

Figure 7.4: From Coker et al.: solvated electron a) in liquid helium, b) in liquid xenon

• Solid parahydrogen [ZOPPI 91]:
Kinetic energy in the lattice from PIMC agrees well with results from neutron scattering.
• Recent survey of PIMC applications: [CEPERLEY 95]; see also the web site [CEPERLEY WWW].

vesely 2006

 < >