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

 < >
 Part III: Ch. 7

## 7.1 Diffusion Monte Carlo (DMC)

First suggested in the Forties (see [MEYERS 56]); rediscovered in the eighties ([CEPERLEY 80]).

Originally applied to the ground state of a bosonic system such as $\scriptstyle 4$He [KALOS 74, WHITLOCK 79]. Later: extended to fermions and excited states ( [BARNETT 86, CEPERLEY 88]); see also [CEPERLEY 96].

Time-dependent Schroedinger equation for a particle of mass $m$ in a potential $U(r)$

$i \hbar \frac{\textstyle \partial \Psi(r,t)}{\textstyle \partial t} = {\rm H} \, \Psi(r,t)$    (7.1)

with the energy operator
${\rm H} \equiv - \frac{\textstyle \hbar^{2}}{\textstyle 2m} \nabla^{2} + \left[ U(r) - E_{T} \right]$    (7.2)

and a trial energy $E_{T}$

Define an imaginary time variable $s \equiv it/\hbar$; then
$\frac{\textstyle \partial \Psi(r,s)}{\textstyle \partial s} = D \nabla^{2} \Psi(r,s) -\left[ U(r) - E_{T} \right] \Psi(r,s)$    (7.3)

with $D \equiv \hbar^{2}/2m$

$\Longrightarrow$ Diffusion with autocatalysis! Visualize $\Psi$ as describing the density of bacteria diffusing in a fluid with locally varying nutrient concentration.

Expand $\Psi$ in eigenfunctions $\Psi_{n}$ of ${\rm H}$:
• If $E_{T}=E_{0}$ (ground state energy), then all $\Psi_{n}$ except $\Psi_{0}$ will fade out for large "times" $s$:

$\lim_{s \rightarrow \infty} \Psi(r,s)=\Psi_{0}(r)$    (7.4)

• If $E_{T} > E_{0}$, the total momentary weight $I(s) \equiv \int\Psi(r,s) \, d r$ will grow exponentially in time.
• If $E_{T} < E_{0}$, the integral $I(s)$ decreases exponentially in time.
$\Longrightarrow$ Solve 7.3 for various values of $E_{T}$: find that $E_{T}$ for which $I(s)$ remains stationary. Then $E_{T}=E_{0}$ and $\Psi = \Psi_{0}$.

How to solve 7.3? Consider the diffusion and the autocatalysis parts of the equation separately.

Diffusion part:

$\frac{\textstyle \partial n(r,t)}{\textstyle \partial t} = D\, \nabla^{2} n(r,t)$    (7.5)

May be solved by a random walk procedure: Start $N$ Brownian walkers

$r_{i}(t_{n+1})=r_{i}(t_{n})+\$\bf\xi _{i}\,,\;\;\;i=1,\dots N $(7.6) where$\xi_{x,y,z}$are drawn from a Gauss distribution with$\sigma^{2}=2D\, \Delta t$. Now consider an ensemble of$M$such$N$-particle systems and write the local distribution density at time$t$as$ p(r,t) \equiv \langle \delta \left[ r_{i}(t)-r \right] rangle = \frac{\textstyle 1}{\textstyle M} \frac{\textstyle 1}{\textstyle N} \sum \limits_{l=1}^{M} \sum \limits_{i=1}^{N} \delta \left[ r_{i,l}(t)-r \right] $(7.7) This is an estimate for the solution$n(r,t)$of the diffusion equation 7.5. Autocatalysis part:$ \frac{\textstyle \partial n(r,t)}{\textstyle \partial t} = f(r) \, n(r,t) $(7.8) Instead of writing the formal solution as$ n(r,t)=n(r,0)\,exp \left[f(r)t\right] $(7.9) we use once more a stochastic method. Consider an ensemble of$M$systems of$N$particles each, at fixed positions. Let the number$M$of systems in the ensemble be allowed to vary: - Systems that contain many particles located at positions with high values of$f(r)$are replicated - Systems with unfavorable configurations are deleted Procedure: Let the ensemble be given at step$t_{n}$. • For each of the$M(t_{n})$systems, determine the multiplicity (see equ. 7.9)$ K_{l}=exp\left[ \sum \limits_{i=1}^{N} f(r_{i,l}) \, \Delta t \right]\,, \;\;\;l=1, \dots M(t_{n}) $(7.10) • Replicate the$l$-th system such that on the average$K_{l}$copies are present. To achieve this, produce first$int (K_{l})-1$copies (where$ int(..)$= next smaller integer) and then, with probability$w \equiv K_{l}- int(K_{l})$, one additional copy. (In practice, draw$\xi$equidistributed$\in [0,1]$and check whether$\xi \leq w$.) If$K_{l} < 1$, remove, with probability$1-K_{l}$, the$l$-th system from the ensemble. Again, the distribution density 7.7 is an estimate of the density at position$r$Combining the two stochastic techniques for solving the diffusion and autocatalytic equations we obtain the following procedure. Figure 7.1: Quantum mechanical diffusion Monte Carlo $N$(non-interacting) particles of mass$m$, distributed at random in a given spatial region, are subject to the influence of a potential$U(r)$. Determine the "diffusion constant"$D=\hbar^{2}/2m$; choose a trial energy$E_{T}$, a time step$\Delta s$and an initial ensemble size$M(s_{0})$. For each system$l$($=1, \dots \, M(s_{0})$) in the ensemble and for each particle$i$($=1, \dots \, N$) perform a random displacement step$ r_{i,l}(s_{n+1})=r_{i,l}(s_{n})+\xi_{i,l}\,,\;\;\; i=1,\dots N $where the components of the vector$\xi_{i,l}$are picked from a Gaussian distribution with$\sigma^{2}=2D\,\Delta s$. For each system$l$determine the multiplicity$K_{l}$according to$ K_{l}=exp\left\{ \left[ \sum_{i=1}^{N} U(r_{i,l})-E_{T} \right] \Delta s \right\} $Produce$\mbox{int}(K_{l})-1$copies of each system ($\mbox{int}(...)$denoting the nearest smaller integer;) with probability$w =K_{l}-\mbox{int}(K_{l})$produce one additional copy, such that on the average there are$K_{l}$copies in all. If$ K_{l}<1 $, purge the system with probability$1-K_{l}$from the ensemble. If the number$M$of systems contained in the ensemble increases systematically (i.e. for several successive steps), choose a smaller$E_{T}$; if$M$increases, take a larger$E_{T}$. Repeat until$M$remains constant; then the ground state energy is$E_{0}=E_{T}$and$ \Psi_{0}(r) = \langle \delta (r_{i,l}-r) \rangle $So far:$\Psi$must be real and$\ge 0\Longrightarrow$Bosons ($\scriptstyle 4 $He or similar). Generalization for fermions: fixed node and released node approximation [CEPERLEY 88]. Note: The analogy between the wave function$\Psi(r,t)$and the local density$n(r,t)$is purely formal. It must be distinguished from the physical interpretation of$|\Psi (r)|^{2}=prob\{ quantum\;object\;to\; be\; found\; at\; r\}$Importance sampling DFT: If the potential$U(r)$is highly negative in some region of space, the autocatalytic term may get out of control and must be handled differently: - Introduce an estimate$\Psi_{T}(r)$of the correct solution$\Psi_{0}(r)$; - define an auxiliary function$ f(r,s) \equiv \Psi_{T}(r)\,\Psi(r,s) $(7.11) - insert this in 7.3 to find$ \frac{\textstyle \partial f}{\textstyle \partial s} = D \nabla^{2} f - \left[ \frac{\textstyle {\rm H} \Psi_{T}}{\textstyle \Psi_{T}} - E_{T} \right] \, f -D \,\nabla \cdot \left[ f \, \nabla \ln | \Psi_{T} |^{2} \right] $(7.12) Since$ \frac{\textstyle {\rm H}\,\Psi_{T}}{\textstyle \Psi_{T}} \approx E_{0} \approx E_{T} $(7.13) the autocatalytic term is now well-behaved, and the multiplicity$K_{l}$will remain bounded. Visualisation of equ. 7.12: The new term looks like an advective contribution. In the image of a diffusing and multiplying bacterial strain there is now an additional driving force$ F(r) \equiv \nabla \, \ln |\Psi_{T}(r)|^{2} $(7.14) which creates a flow, or drift. This means that the individual diffusors follow a preferred direction along$F(r)$:$ r_{i,l}(s_{n+1})=r_{i,l}(s_{n})+ \xi_{i,l}+D \Delta s \, F(r_{i,l}(s_{n})) $(7.15) The multiplicity$K_{l}$is now$ K_{l}=exp \left\{ \left[ \frac{\textstyle {\rm H} \Psi_{T}}{\textstyle \Psi_{T}} - E_{T} \right] \, \Delta s \right\} \$    (7.16)

Green's function Monte Carlo (GFMC): Another formulation of the DMC procedure [SKINNER 85].

Recent literature: See [CEPERLEY 96] and web sites [CEPERLEY WWW] or [CAVENDISH WWW].

vesely 2006

 < >