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

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})$.
  1. 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$.

  2. 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\} $

  3. 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.

  4. 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}$.

  5. 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)


$ \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

< >