Monte Carlo normalization of a wave function

Question

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.

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