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?
- Density matrix for the free particle:
One-dimensional box of length $L$; particle mass $m$;
energy operator
$
{\rm H}=-\frac{\textstyle \hbar^{2}}{\textstyle 2m} \,
\frac{\textstyle \partial^{2}}{\textstyle \partial x^{2}}
$
(7.20)
The eigenfunctions of ${\rm H}$ normalized over $[-L/2,L/2]$ are
$
\Psi_{n}=\frac{\textstyle 1}{\textstyle \sqrt{L}}\,
e^{\textstyle ik_{n}x}\,,\;\;\;\;{\rm with}\;\;\;
k_{n}\equiv\frac{\textstyle 2\pi n}{\textstyle L}\;\;\;\;{\rm and}\;\;\;
E_{n}=\frac{\textstyle \hbar^{2}}{\textstyle 2m}\,k_{n}^{2}
$
(7.21)
yielding the density matrix
$
\rho_{0}(x,x';\beta)=\sum_{n}\frac{\textstyle 1}{\textstyle L}
e^{\textstyle -ik_{n}(x-x')}
\,e^{\textstyle -\beta \hbar^{2} k_{n}^{2}/2m}
$
(7.22)
Now let $L \rightarrow \infty$, treating the discrete wave number
$k_{n}$ as a continuous variable $k$ with differential
$dk \approx \Delta k =k_{n+1}-k_{n}=2 \pi /L$. Then
$
\rho_{0}(x,x';\beta)=\frac{\textstyle 1}{\textstyle 2\pi}
\int\limits_{-\infty}^{\infty}
e^{\textstyle -ik(x-x')}\,e^{\textstyle -\beta \hbar^{2} k^{2}/2m}\, dk
$
(7.23)
Thus we find for the density matrix of the free particle
$
\rho_{0}(x,x';\beta)=
\left[ \frac{\textstyle m}{\textstyle 2\pi \beta \hbar^{2}}\right]^{1/2}
e^{\textstyle -m(x-x')^{2}/2\beta \hbar^{2}}
$
(7.24)
For the diagonal element we have
$
\rho_{0}(x,x;\beta)=
\left[ \frac{\textstyle m}{\textstyle 2\pi \beta \hbar^{2}}\right]^{1/2}
$
(7.25)
(independent of $x$).
- Density matrix for the harmonic oscillator:
Particle mass $m$; harmonic potential
$
U(x)=\frac{\textstyle m \omega_{0}^{2}}{\textstyle 2}\, x^{2}
$
(7.26)
Density matrix
[KUBO 71]:
$
\begin{eqnarray}
\rho(x,x';\beta)&=&
\left[ \frac{\textstyle m \omega_{0}}{\pi \hbar}\, tanh \frac{\textstyle \beta \hbar \omega_{0}}{2}
\right]^{1/2}
\\
&&
\cdot exp \left\{ -\frac{\textstyle m \omega_{0}}{4 \hbar}
\left[ (x+x')^{2} \, tanh \frac{\textstyle 1}{2} \beta \hbar \omega_{0}
+(x-x')^{2} \, coth \frac{\textstyle \beta \hbar \omega_{0}}{2} \right] \right\}
\end{eqnarray} $
(7.27)
with diagonal elements
$
\displaystyle
\rho(x,x;\beta)=
\left[ \frac{\textstyle m \omega_{0}}{\pi \hbar}\, tanh \frac{\textstyle \beta \hbar \omega_{0}}{2}
\right]^{1/2}\, exp \left\{ -\frac{\textstyle m \omega_{0}}{\hbar}
x^{2} \, tanh \frac{\textstyle \beta \hbar \omega_{0}}{2} \right\}
$
(7.28)
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$:
$
\rho(x_{0},x_{P};\beta)=\int \dots \int dx_{1}\, dx_{2}\dots dx_{P-1}
\;\rho(x_{0},x_{1};\frac{\textstyle \beta}{\textstyle P})\dots
\rho(x_{P-1},x_{P};\frac{\textstyle \beta}{\textstyle P})
$
(7.30)
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)$:
$
\begin{eqnarray}
\rho(x_{p}, x_{p+1}; \frac{\textstyle \beta}{\textstyle P})
& \longrightarrow & \sum_{n}\Psi_{n}^{*}(x_{p})\,
e^{\textstyle -(\beta/P)\left[{\rm H}_{free}+U(x) \right]}\,
\Psi_{n}(x_{p+1})
\\
&&
\approx \sum_{n}\Psi_{n}^{*}(x_{p})\,
e^{\textstyle -(\beta/P){\rm H}_{free}} \Psi_{n}(x_{p+1})
e^{\textstyle -(\beta/2P)\left[U(x_{p})+U(x_{p+1}) \right]}
\\
&&
= \rho_{0}(x_{p},x_{p+1};\frac{\textstyle \beta}{\textstyle P})\,
e^{\textstyle -(\beta/2P)\left[U(x_{p})+U(x_{p+1}) \right]}
\end{eqnarray} $
(7.31)
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
|