3
$\begingroup$

Im doing a Metropolis Monte Carlo simulation with particles on a sphere and have a question concerning the random movement in a given time step.

I understand that to obtain a uniform distribution of random points on the sphere to begin with it is not enough to use the naive simplest way (use spherical coordinates with constant R and pick random angles theta and phi), but one has to for example use one of these methods: http://mathworld.wolfram.com/SpherePointPicking.html

Looking at another code for a Monte Carlo on a sphere I see a fairly complicated way to generate random moves: pick a random particle, calculate the rotation matrix moving it to the north pole, find a random cartesian vector less than a certain length and move it to the north pole, normalize the cartesian vector and then rotate it back to the vicinity of the original particle position.

This is all to get an unbiased new random position. I don`t understand the rationale completely although I suspect it is connected to detailed balance. But my feeling is that there should be an easier (i.e. faster) way to do this. Actually, intuitively I feel that in this case it is ok to find two small random angles theta and phi and add them to the position of the particle - or would this be a mistake?

  • 0
    Sorry, didn`t know that was frowned upon. Actually someone in the programming forum suggested I try posting the question here.2012-10-15

2 Answers 2

2

$\newcommand{\+}{^{\dagger}}% \newcommand{\angles}[1]{\left\langle #1 \right\rangle}% \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}% \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}% \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}% \newcommand{\dd}{{\rm d}}% \newcommand{\ds}[1]{\displaystyle{#1}}% \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}% \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}% \newcommand{\fermi}{\,{\rm f}}% \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}% \newcommand{\half}{{1 \over 2}}% \newcommand{\ic}{{\rm i}}% \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow}% \newcommand{\isdiv}{\,\left.\right\vert\,}% \newcommand{\ket}[1]{\left\vert #1\right\rangle}% \newcommand{\ol}[1]{\overline{#1}}% \newcommand{\pars}[1]{\left( #1 \right)}% \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}}% \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}% \newcommand{\sech}{\,{\rm sech}}% \newcommand{\sgn}{\,{\rm sgn}}% \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}}% \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ Let's assume we have a mechanism to generate random numbers in $\left[0,1\right)$ and ${\rm P}\pars{\Omega_{\vec{r}}}$ is a distribution function for random points in a sphere of radius $a > 0$. $\Omega_{\vec{r}}$ is the solid angle. In this case, ${\rm P}\pars{\Omega_{\vec{r}}}$ is, indeed, $\Omega_{\vec{r}}\,\,$-independent: $ 1 = \int{\rm P}\pars{\Omega_{\vec{r}}}\,\dd\Omega_{\vec{r}} = {\rm P}\pars{\Omega_{\vec{r}}}\int\dd\Omega_{\vec{r}} = {\rm P}\pars{\vec{r}}\pars{4\pi} \quad\imp\quad{\rm P}\pars{\vec{r}} = {1 \over 4\pi} $ Then, $ 1 = \int_{0}^{\pi}\half\,\sin\pars{\theta}\,\dd\theta\int_{0}^{2\pi} \,{\dd\phi \over 2\pi} $ We can generate random numbers $\xi_{\theta}$ and $\xi_{\phi}$ such that: $ \bracks{~\half\,\sin\pars{\theta}\,\dd\theta = \dd\xi_{\theta}\,, \quad\xi_{0} = 0 \imp \theta = 0~}\ \mbox{and}\ \bracks{~{\dd\phi \over 2\pi} = \dd\xi_{\phi}\,,\quad\xi_{0} = 0 \imp \phi = 0~} $ Those relations yield: $\ds{\half\bracks{-\cos\pars{\theta} + 1} = \xi_{\theta}}$ $\ds{\pars{~\mbox{or/and}\ \sin\pars{\theta/2} = \root{\xi_{\theta}}~}}$ and $\ds{\phi = 2\pi\,\xi_{\phi}}$: $\color{#0000ff}{\large% \theta = 2\arcsin\pars{\root{\xi_{\theta}}}\,,\qquad \phi = 2\pi\xi_{\phi}} $ As we mentioned above, $\xi_{\theta}$ and $\xi_{\phi}$ are uniformly distributed in $\left[0, 1\right)$.

For a sphere of radio $a$ the random points are given by: $ \left\lbrace% \begin{array}{rcl} \color{#0000ff}{\large x} & = & a\sin\pars{\theta}\cos\pars{\phi} = \color{#0000ff}{\large 2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}\cos\pars{2\pi\xi_{\phi}}} \\ \color{#0000ff}{\large y} & = & a\sin\pars{\theta}\sin\pars{\phi} = \color{#0000ff}{\large 2a\root{\xi_{\theta}\pars{1 - \xi_{\theta}}}\sin\pars{2\pi\xi_{\phi}}} \\ \color{#0000ff}{\large z} & = & a\cos\pars{\theta} = \color{#0000ff}{\large a\pars{1 - 2\xi_{\theta}}} \end{array}\right. $
1

I don't clearly see what the second paragraph (uniform distribution) has to do with the rest (unbiased random steps). Anyway: in classical Montecarlo the choosing of displacement is not critical, typically one is only concerned with:

  1. sampling the full configuration space (this typicaly means that the displacement is random enough in "direction", so that any configuration is a priori reachable)

  2. having a "nice" probability (not almost one, nor almost zero) of move acceptance (this typically means that the average displacement is not very large, nor very small)

In your case, I'd try the simplest alternative: for a selected point choose a (random or fixed) displacement, and choose a random direction (a uniform angle). The "fairly complicated" procedure you describe seem to be just a implementation of this approach.

  • 0
    @StevenStadnicki Thanks for the input! I actually changed method after the postings above; the new one is detailed in an answer here: http://stackoverflow.com/questions/12900861/monte-carlo-simulation-on-sphere-unbiased-random-steps (I didn't know double posting was bad at the time). Using this I get the same results as the method described in the question above.2013-09-12