2

I would like to normalize a quantum mechanical multi-particle wave function numerically, and since the result is a multidimensional integral I thought Monte Carlo methods might be appropriate. So, I'm looking for $\mathcal{N}$ where

$\psi=\frac{1}{\sqrt{\mathcal{N}}}\tilde{\psi}$

and $\psi$ and $\tilde{\psi}$ are the normalized and unnormalized wave function respectively. So, introducing $\langle f(\textbf{x})\rangle$ as the expectation value of $f$ in the normalized state $\psi$, with position coordinates $\textbf{x}$ (a vector containing all components of position vectors for all particles), I find the Monte Carlo estimate

$\frac{1}{\mathcal{N}}=\frac{1}{\mathcal{N}A^N_d}\int d^d\textbf{x}=\frac{1}{A^N_d}\int d^d\textbf{x}\big|\frac{\psi(\textbf{x})}{\tilde{\psi}(\textbf{x})}\big|^2\approx\frac{1}{A^N_dn}\sum_{j=1}^n\frac{1}{|\tilde{\psi}(\textbf{x}_j)|^2}$

where $A_d$ is the (finite) $d$-dimensional area of the system and $\textbf{x}_j$ are $n$ position configurations drawn from the (normalized) probability density $|\psi|^2$. I'm using the Metropolis algorithm, which automatically uses a normalized probability, so that I have no choice in that matter.

But when I try to do this I get an estimate for $\frac{1}{\mathcal{N}}$ which increases with number of Monte Carlo samples $n$. And the strange thing is, so does the standard deviation, which one expects to have a behaviour $\sim n^{-1/2}$ in a Monte Carlo. I've checked this up to millions of configurations, for a wave function $\psi$ that is expected to converge relatively fast.

Can anyone find some mistake in my logic? Or does someone have a tip for a better way to calculate the normalization?

EDIT: In reply to a comment by Isidore Seville I'll explain why I need the normalization of $\tilde{\psi}$ in the first place.

My wave function is itself given as an integral over auxiliary coordinates $\textbf{y}$:

$\tilde{\psi}(\textbf{x})=\int d^d\textbf{y}\phi_1(\textbf{y})\phi_2(\textbf{y},\textbf{x})$

When performing a Metropolis Monte Carlo using $|\tilde{\psi}(\textbf{x})|^2$ as the probability density I basically have two possibilities:

  1. For every random move in the particle coordinates $\textbf{x}$, perform a MC integration in $\textbf{y}$ to get the resulting $|\tilde{\psi}(\textbf{x})|^2$
  2. Introduce another set of auxiliary coordinates $\textbf{z}$ and use $|\tilde{\psi}(\textbf{x})|^2=\int d^d\textbf{y}d^d\textbf{z}\ \phi_1(\textbf{y})^*\phi_1(\textbf{z})\phi_2(\textbf{y},\textbf{x})^*\phi_2(\textbf{z},\textbf{x})$. Then treat the MC as an integral over all three sets of coordinates, randomly moving $\textbf{x}$, $\textbf{y}$ and $\textbf{z}$ for every iteration.

For now I decided to go for the second possibility, which for e.g. the expectation value of a function $f(\textbf{x})$ thus entails calculating

$\langle f(\textbf{x})\rangle=\int d^d\textbf{x}f(\textbf{x})|\tilde{\psi}(\textbf{x})|^2=\int d^d\textbf{x}d^d\textbf{y}d^d\textbf{z}\ f(\textbf{x})\phi_1(\textbf{y})^*\phi_1(\textbf{z})\phi_2(\textbf{y},\textbf{x})^*\phi_2(\textbf{z},\textbf{x})$

But this introduces a new problem: although of course $|\tilde{\psi}(\textbf{x})|^2$ is real, the integrand above isn't necessarily until the integrals over $\textbf{y}$ and $\textbf{z}$ converge. Therefore I have to use another probability density $P(\textbf{x},\textbf{y},\textbf{z})$ that is real and nonnegative to find the MC estimate (using some shorthand below):

$\langle f(\textbf{x})\rangle=\int d\Omega\ f\phi_1^*\phi_1\phi_2^*\phi_2\frac{P}{P}\approx\frac{1}{n}\sum_{i=1}^N\frac{f\phi_1^*\phi_1\phi_2^*\phi_2}{P}$

where the summand above is evaluated on configurations drawn from the probability $P$ ($n$ in total).

So this is the situation: To evaluate the summand I need the normalized $\phi_1$, $\phi_2$ and $P$ (the latter since there is no way to do the Metropolis MC with an unnormalized $P$). Also, there is no way to just do the integral using $|\tilde{\psi}|^2$ in the normal way, since this itself includes an integral.

jorgen
  • 123
  • 4
  • 1
    I don't understand your method. Why aren't you simply integrating the unnormalized wavefunction with MC integration? i.e. $1/N = \int |\tilde\psi|^2 dx$, then MC integrate the RHS. –  Feb 04 '14 at 14:31
  • Without any importance sampling (i.e. just take $x$ randomly, uniform distribution) I see how that might work, but I'm assuming that would converge slowly. I'd like to use $\sim|\psi|^2$ itself as a probability density, and in order to do that I would need to do $\int|\tilde{\psi}|^2dx=\int|\tilde{\psi}/\psi|^2\cdot|\psi|^2dx\approx\sum_{i=1}^n|\tilde{\psi}/\psi|^2(x_i)$. The summand can't be evaluated since then I would need $\mathcal{N}$ itself. And the Metropolis algorithm automatically uses the normalized $|\psi|^2$, so I don't see any way around. Or am I misunderstanding something? –  Feb 04 '14 at 16:46
  • This sounds more like a question for [scicomp.SE]. (We'll migrate it there if people agree.0 – David Z Feb 04 '14 at 16:56
  • Huh, I didn't know about that one, sounds very useful. I can move it there myself –  Feb 04 '14 at 17:08
  • 1
    Hello, jorgen. I don't understand why you would need the normalization factor. From my limited experience, usually only the expectation values matter, and these quantities can be obtained in MC with standard metropolis algorithm. Could you explain a little why you want to compute the normalization factor? –  Feb 05 '14 at 05:09
  • I need the normalization for a different MC where I calculate the energy of a wave function $\phi$ but use $\tilde{\psi}$ as probability: $\langle E\rangle=\int E|\phi|^2dx=\int E|\phi/\psi|^2|\psi|^2dx\approx\sum_{i=1}^nE(x_i)|\phi(x_i)/\psi(x_i)|^2/n$. I need the normalized $\psi$ to calculate the summand since the MC automatically use it instead of $\tilde{\psi}$ even though I use the latter in the code. –  Feb 05 '14 at 08:23

1 Answers1

1

My understanding is that, you have computed the energy of a given wave function by doing the following:

$$ E[\phi] = \langle \phi|H|\phi\rangle = \frac{\int E(x) |\phi|^2 dx}{\int |\phi|^2 dx} = \frac{\int E(x) (|\phi|^2/|\psi|^2) |\psi|^2 dx}{\int |\psi|^2 dx}\times \frac{\int |\psi|^2 dx}{\int |\phi|^2 dx}. $$

Then, the following quantity is measured in MC:

$$ A = \frac{\int E(x) (|\phi|^2/|\psi|^2) |\psi|^2 dx}{\int |\psi|^2 dx}, $$

by treating $|\psi|^2$ as unnormalized probability and $E(x) (|\phi|^2/|\psi|^2)$ as a variable. Now we need the value of the following quantity:

$$ B = \frac{\int |\psi|^2 dx}{\int |\phi|^2 dx} = \frac{\int (|\psi|^2/|\phi|^2) |\phi|^2 dx}{\int |\phi|^2 dx}, $$

We can measure $B$ by treating $|\psi|^2/|\phi|^2$ as a variable. To make it work, the reference state $\psi$ and the state $\phi$ should be "reasonably close".

I hope the above helps.

Isidore Seville
  • 111
  • 1
  • Thanks! But if I understand correctly the above method would require me to calculate $B$ in a MC drawing configurations from $|\phi|^2$ - and this is intractable (the reason for the whole multiply/divide by $|\psi|^2$ in the first place). – jorgen Feb 05 '14 at 14:58
  • @jorgen Could you explain a bit why you cannot sample $\phi$ directly? This is quite unusual from my experience in variational MC. – Isidore Seville Feb 05 '14 at 16:52
  • I edited the OP to answer your question. I'd be interested to hear if you have another suggestion on how to do it! – jorgen Feb 05 '14 at 18:09
  • @jorgen Thanks for updating. It seems that the kernel $\phi_2(x,y)$ can be thought of some operator acting on the reference state $\phi_1(x)$. What is the structure of $\phi_2$? If it is not "sparse" or has simple, analytic expression, then doing MC is indeed a challenging task. – Isidore Seville Feb 05 '14 at 19:33
  • $\phi_2$ is not sparse but it does have a relatively simple analytic expression. That was why I originally thought the MC might have a chance to work, but I'm kinda stuck on this normalization.. – jorgen Feb 05 '14 at 20:02