2021/08/28 00:48 更新
拡散量子モンテカルロ法
目次
$$\gdef\pd#1{\frac{\partial}{\partial #1}} \gdef\pdv#1#2{\frac{\partial #2}{\partial #1}} \gdef\pndv#1#2#3{\frac{\partial^{#3} #2}{\partial #1^{#3}}} $$

この記事について

この記事はシュレディンガー方程式を数値的に解く手法のうちの一つである「拡散量子モンテカルロ法」の原理と実際の計算例についてまとめたものです。記事の作成に当たって、https://arxiv.org/pdf/physics/9702023.pdf を大いに参考にしました。このような数値計算手法をまとめた日本語の記事があまり見当たらなかったので、実際に遊んでみるがてら執筆することにした次第です。とはいえ私は趣味で数値計算をやってみただけの素人ですので、より詳しく正確に知りたい方は上記論文を読むのが早いでしょう。
 よく知られているように、シュレディンガー方程式は原子や分子といったミクロな系を記述するとても強力な式ですが、厳密に解けるものは水素原子といった特殊な場合に限られます。特に、多体系は解析的に解くのはほぼ不可能で、これまでに多体系の性質を調べるための様々な近似法や数値計算手法が考案されてきました。今回紹介する拡散量子モンテカルロ法は数値計算手法の一つです。この手法は、方程式そのものを数値的に解くのではなく、波動関数に見立てた仮想的な粒子を大量に設置し、それらをランダムに移動させて定常状態に落ち着かせることにより、基底状態の波動関数及びエネルギー固有値を得るというものです。ただし、波動関数が実数で表される系にしか応用できないという制限があります。
 この記事では、1章で1次元1粒子系の虚時間のシュレディンガー方程式を導入し、2章で経路積分による波動関数の(虚)時間発展を数値的に計算する方法を見ます。これに基づいて3章で拡散量子モンテカルロ法を解説し、4, 5章で多次元多粒子系への拡張と数値計算をするための無次元化について話します。最後に6章では、実際の数値計算の例を、1次元調和振動子、水素原子、ヘリウム原子について載せ、計算の結果が理論値/実験値によく合うことを見ます(図表はあまりないですが・・・)。

1. 1次元虚時間シュレディンガー方程式

1次元シュレディンガー方程式

まず、1次元1粒子のシュレディンガー方程式から出発することにします。そのあとに$d$次元、多粒子系に拡張します。質量$m$の粒子が従うシュレディンガー方程式はハミルトニアンを$\hat{H}$, 波動関数を$\Psi(x,t)$として

$$\begin{aligned} \tag{1} &i\hbar \pdv{t}{\Psi} = \hat{H}\Psi \\ &\hat{H} =-\frac{\hbar^2}{2m}\pndv{x}{}{2}+V(x) \end{aligned}$$

と書けます。一般に$\Psi(x,t)$は、$\hat{H}$のエネルギー固有関数$\phi_n(x)$と固有値 $E_n$を用いて、次のように展開することができます。

$$\tag{2} \Psi(x,t)=\sum_{n=0}^{\infty} c_n\phi(x)e^{-\frac{i}{\hbar}E_nt}$$

ここで$c_n$は展開係数です。$n=0, 1, 2, 3\cdots$に対して、エネルギー固有値$E_n$は

$$\tag{3} E_0 < E_1 \leq E_2\cdots$$

と順序付けされているものとします。


エネルギーのシフト

数値計算でシュレディンガー方程式を安定的に解く際のテクニックとして、エネルギーをシフトしておきます。具体的には、$V(x)\rightarrow V(x)-E_R, ~ E_n\rightarrow E_n-E_R$とエネルギーを$E_R$分ずらしておきます。これに伴ってシュレディンガー方程式(1)および式(2)は

$$\tag{4} i\hbar \pdv{t}{\Psi}= \left[ -\frac{\hbar^2}{2m}\pndv{x}{}{2} +\left( V(x)-E_R \right) \right] \Psi$$$$\tag{5} \Psi(x,t)=\sum_{n=0}^{\infty}c_n\phi_n(x)e^{-i\frac{E_n-E_R}{\hbar}t}$$

のように変更されます。


実時間から虚時間への変換

虚時間と聞くと何か仰々しいものと思われるかもしれませんが、特に難しいことはなく、単に$\tau = it$という変換を施すだけです。この変換で式(4)および式(5)は

$$\tag{6} \hbar \pdv{\tau}{\Psi}= \left[ \frac{\hbar^2}{2m}\pndv{x}{}{2} - \left( V(x)-E_R \right) \right] \Psi$$$$\tag{7} \Psi(x,\tau)=\sum_{n=0}^{\infty}c_n\phi_n(x)e^{-\frac{E_n-E_R}{\hbar}\tau}$$

となります。

 $\tau \rightarrow \infty$のときの波動関数の挙動を調べてみます。式(3)からエネルギー固有値は順序付けされているので、$E_R$と$E_0$の大小で振る舞いが異なってきます。式(7)より、

$$ \lim_{\tau\rightarrow \infty}\Psi(x,\tau)= \left\{ \begin{matrix} 0 &(E_R>E_0) \\ \infty & (E_R < E_0) \\ c_0\phi_0(x) & (E_R=E_0) \end{matrix} \right.$$

となり、$E_R=E_0$のときにのみ、波動関数は収束して基底状態に落ち着くことがわかります。拡散量子モンテカルロ法では、$E_R$を適切に選ぶことによって最終的に基底状態の波動関数及びエネルギー固有値を求めます。

2. 経路積分による定式化

拡散量子モンテカルロ法の方法論を議論するときは、拡散方程式とのアナロジーから出発するやり方と、ファインマンの経路積分から出発するやり方があるそうです。ここでは、私が親しみ深かった経路積分による方法で議論したいと思います。

経路積分

 経路積分法によると、$\Psi(x, \tau)$は時間$\tau = 0$における波動関数$\Psi(x_0,0)$を用いて

$$ \tag{8} \Psi(x,\tau) = \int_{-\infty}^{\infty} dx_0 K(x,\tau|x_0,0)\Psi(x_0,0)$$

と書き表すことができます。ここで$K(x,\tau|x_0,0)$はファインマン核と呼ばれるもので、次のように書き表せます。

$$ \tag{9} K(x,\tau|x_0,0)= \lim_{N \rightarrow \infty} \int_{-\infty}^{\infty}dx_1 \cdots \int_{-\infty}^{\infty}dx_{N-1} \left( \frac{m}{2\pi \hbar \Delta \tau} \right)^{N/2} \exp{ \left[ -\frac{\Delta \tau}{\hbar} \sum_{j=0}^{N} \left\{ \frac{m}{2\Delta \tau^2}(x_j-x_{j-1})^2 + V(x_j)-E_R \right\} \right] }$$

ただし、$\Delta \tau=\tau/N$です。ここで、次のように$P(x_n, x_{n-1})$と$W(x_n)$を導入します。

$$ \tag{10} P(x_n,x_{n-1})=\left( \frac{m}{2\pi \hbar \Delta \tau} \right)^{1/2} \exp{\left\{ -\frac{m(x_n-x_{n-1})^2}{2\hbar \Delta \tau} \right\}}$$$$ \tag{11} W(x_n)=\exp{ \left\{ -\frac{\Delta \tau(V(x_n)-E_R)}{\hbar} \right\} }$$

これらの式を使って式(8)を書き直すと、

$$\tag{12} \Psi(x, \tau)=\lim_{N\rightarrow \infty} \left( \prod_{j=0}^{N-1} \int_{-\infty}^{\infty} dx_j \right) \prod_{n=1}^{N}W(x_n)P(x_n,x_{n-1}) \Psi(x_0, 0)$$

とります。式(10)をよく見てみると、これは分散$\sigma^2 = \hbar \Delta \tau /m$の正規分布になっています。また、次の関係式を満たすことがわかります。

$$\tag{13} \int_{-\infty}^{\infty}dyP(x,y)=1$$$$\tag{14} \int_{-\infty}^{\infty} dy P(x,y)P(y,z)=P(x,z)$$

モンテカルロ積分

さて、ここで$\mathcal{P}(x_0,\cdots, x_{N-1})$を次のように定義します。

$$\tag{15} \mathcal{P}(x_0,\cdots,x_{N-1})= \prod_{n=1}^{N} P(x_n,x_{n-1})$$

この$\mathcal{P}(x_0,\cdots,x_{N-1})$は式(13)と式(14)から

$$\tag{16} \left( \prod_{j=0}^{N-1} \int_{-\infty}^{\infty} dx_j \right) \mathcal{P}(x_0,\cdots x_{N-1}) =1$$

を満たし、多変数の確率密度関数とみなすことができます。この性質より、十分大きな$N$に対して$\Psi(x, \tau)$は次のように近似できます。

$$\begin{aligned} \tag{17} \Psi(x,\tau) &\sim \left( \prod_{j=0}^{N-1} \int_{-\infty}^{\infty} dx_j \right) \mathcal{P}(x_0,\cdots x_{N-1}) f(x_0,\cdots,x_{N-1}) \\ &\sim \frac{1}{\mathcal{N}} \sum_{\substack{i=1\\x^{(i)} \in \mathcal{P}}}^{\mathcal{N}} f(x_0^{(i)},\cdots, x_{N-1}^{(i)}) \end{aligned}$$

ここで、

$$ \tag{18} f(x_0^{(i)},\cdots, x_{N-1}^{(i)}) = \Psi(x_0, 0) \prod_{n=1}^NW(x_n)$$

としました。式(17)の $x^{(i)}\in \mathcal{P}$は、$\mathcal{P}(x_0,\cdots x_{N-1})$を確率密度関数とする確率分布からサンプリングしたもので、いわゆるモンテカルロ積分によって積分を評価しています。したがって得られる積分の値は厳密なものではありませんが、誤差は$\mathcal{O}(1/\sqrt{\mathcal{N}})$に比例して小さくなります。
 $\mathcal{P}$に従う確率分布から$(x_0,\cdots x_{N-1})$をサンプリングするには次のようにします。そもそも、$\mathcal{P}(x_0,\cdots x_{N-1})$の定義は式(15)でした。まず初期位置$x_0$から初めて、平均$x_0$、分散$\sigma^2 = \hbar \Delta \tau /m$に従う正規分布から $x_1$を発生させます。次に、平均$x_1$,分散$\sigma$に従う正規分布から$x_2$を発生させます。・・・というふうにこれを$x_N$まで繰り返します。通常は、毎回平均値を変えながらサンプリングするのではなく、標準正規分布(平均$0$, 分散$1$)からサンプリングした数$\rho_n$を使って次のように順々に発生させます。

$$ \tag{19} x_n=x_{n-1}+\sigma \rho_n$$

3. 拡散量子モンテカルロ

上記の方法によって$\Psi(x, \tau)$を計算することができますが、毎回与えられた$x,\tau$に対して計算することになり、基底状態の波動関数とエネルギー固有値を求めるにはあまり実際的ではありません。しかし、上記の経路積分の方法を改善することによって、基底状態の波動関数とエネルギーを同時に求めることができます。それが今回紹介する拡散量子モンテカルロ法というものです。

考え方

 拡散量子モンテカルロ法の考え方は大雑把に言うと次のようなものです:空間に「粒子」を大量に生成した後、それらをランダムウォークで十分な時間拡散させて、最終的な粒子の分布を波動関数(確率分布)とみなす。基底状態のエネルギーの求め方については後ほど説明します。ここでいう「粒子」というのは、実際に物理的に興味があるミクロな粒子のことではなく、シミュレーション上において仮想的に作る粒子です。この拡散モンテカルロ法では先の説明の通り、この仮想的な粒子の位置の分布そのものを波動関数とみなすために、基底状態の波動関数が常に正となる系にしか応用できないという制約があります。
 拡散モンテカルロ法では式(12)の被積分関数部分を次のように解釈し、次のように手順$0\sim2N$を順次実行していきます。

$$ \tag{20} \underbrace{\Psi(x_0, 0)}_{手順~0} \underbrace{P(x_1, x_0)}_{手順~1} \underbrace{W(x_1)}_{手順~2} \cdots \underbrace{P(x_N,x_{N-1})}_{手順~2N-1} \underbrace{W(x_N)}_{手順~2N} $$

初期状態の決定

手順0では、ランダムウォークさせる「粒子」を空間に$\mathcal{N}$個配置します。この初期状態は波動関数$\Psi(x_0,0)$に対応し、通常$\Psi(x_0,0)=\delta(x-x_0)$とします。この状態はある位置$x=x_0$に全粒子を配置することで実現されます。


ランダムウォーク

手順1以降では、配置した$\mathcal{N}$個の粒子を式(19)に従ってそれぞれ移動させます。すなわち手順1では、各$i$番目の粒子に対して、$x_1^{(i)}=x_0^{(i)}+\sigma\rho_0^{(i)}$ のように移動させます。


仮想的な粒子の生成・消滅過程

手順2では、粒子の位置$x_1^{(i)}$ごとに$W(x_1^{(i)})$を計算するのではなく、その代わりに$W(x_1^{(i)})$に応じて粒子を生成させたり消滅させたりするという少しトリッキーなことを行います(実際、私もなぜこのようなトリックを施すかは完全に理解できていませんが、非常にうまく行きます)。$W(x_1^{(i)})$(式(14))を計算するにあたって、$E_R$をどうするかという問題がありますが、先に$W(x_1^{(i)})$が計算できたとして方法論を説明します。この過程では、$W(x_1^{(i)})$の値を使って、次で決定される整数値$m_1$を計算します。

$$\tag{21} m_1 = \min\left( \mathrm{int}{ \left[W(x_1^{(i)})+u\right],3 } \right)$$

ここで$\mathrm{int}[x]$は$x$の整数部分、$u$は区間$[0,1]$からランダムにサンプリングした実数です。そして、得られた$m_1$に対して次のように粒子を生成もしくは消滅させます。

  • $x_1^{(i)}=0$ ならば、粒子$i$を消滅させる。
  • $x_1^{(i)}=1$ ならば、粒子$i$をそのままにする(何もしない)。
  • $x_1^{(i)}=2$ ならば、粒子$i$と同じ位置に粒子を$1$つ生成(コピー)する。
  • $x_1^{(i)}=3$ ならば、粒子$i$と同じ位置に粒子を$2$つ生成(コピー)する。

このようにする直感的な理由としては、ポテンシャルエネルギーの低いところでは粒子の存在確率が大きく、逆に高いところでは小さいということを反映するための処方だと考えられます。また、$m_1\leq3$に制限しているのは、数値計算を安定させるためだそうです。


$E_R$ の決定

$E_R$が基底状態のエネルギー$E_0$に一致しないと、波動関数が基底状態に収束しないことはすでに話しました。$E_R$が$E_0$より大きいと波動関数は$0$に収束し、逆に$E_R$が$E_0$より小さいと波動関数は発散します。今の拡散モンテカルロ法に置き換えると、生成・消滅の過程を経て粒子の数が$0$または$\infty$になることを意味します。うまいこと粒子数を一定に保つために$E_R$をどのように定めたらよいでしょうか。
 粒子の生成・消滅の基準となる$m_1$は、式(21)の$W(x_1^{(i)})$に大きく左右されます。粒子数を一定にするには、$W$が$1$に近い方がよいので、式(11)から$V(x_1^{(i)})\approx E_R$となるのが望ましいです。$V(x_1^{(i)})$の値は各粒子の位置で異なるので、ここでは平均をとります。すなわち、

$$ \tag{22} E_R =\langle V\rangle_1 = \frac{1}{\mathcal{N}_1} \sum_{j=1}^{\mathcal{N}_1} V(x_1^{(i)})$$

のように$E_R$を取ります。各生成・消滅過程で粒子数は変化するので、そのたびに$E_R$の値を計算します。
 しかし、これだけだと粒子の分布が基底状態から遠い場合に粒子数が急激に変化する場合があります。これを防ぐために、式(22)で得た$E_R$に補正を加えます。具体的には正の実数$\alpha$を適切に選んで、

$$ \tag{23} E_R'=E_R+\alpha \left( 1-\frac{\mathcal{N}_1}{\mathcal{N}_0} \right)$$

と補正します。実際の$W(x_1^{(i)})$の計算では、この$E_R'$を使います。$\alpha$としては、

$$ \tag{24} \alpha = \frac{\hbar}{\Delta \tau}$$

と選ぶことが多いようです。


手順$2N$まで繰り返す。

上記の議論では、手順$0$から$2$までの方法を例に説明してきました。あとは手順$1$と手順$2$を手順$2N$まで繰り返すだけです。


波動関数

十分時間がたった後の粒子の位置分布は波動関数そのものとみなすことができます。空間を$x_\mathrm{min}$から$x_\mathrm{max}$までに限定し、$n_b$等分することで、各区間の粒子の個数を度数分布で表します。各区間の粒子の個数を$N_i$ $(i=1,\cdots,n_b)$とすると、波動関数$\phi(x)$は

$$ \tag{25} \phi(x_i)=\frac{ N_i }{ \sqrt{w\sum_{i=1}^{n_b}N_i^2} }$$

と評価できます。ここで、$w=(x_\mathrm{max}-x_\mathrm{min})/n_b$であり、度数分布のビンの幅に対応します。

4. $d$次元$S$粒子系への拡張

これまでは、1次元1粒子系で議論を進めてきました。これを$d$次元$S$粒子系に拡張します。

$d$次元への拡張

$d$次元$1$粒子系におけるハミルトニアンは、

$$ \tag{26} \hat{H}=-\frac{\hbar^2}{2m} \sum_{\alpha=1}^{d}\pndv{x_\alpha}{}{2} +V(x_1,\cdots x_d)$$

となります。拡散モンテカルロ法を実行するときには、$1$次元系のときと同じように、粒子を$d$次元それぞれの方向にランダムウォークさせます。$E_R$を計算する際は、式(21)においてポテンシャルエネルギーをそのまま$V(x_1,\cdots , x_d)$に置き換えれば大丈夫です。


S粒子系への拡張

$d$次元$S$粒子系のハミルトニアンを書くと、

$$ \tag{27} \hat{H}= \sum_{j=1}^{S} \left[ -\frac{\hbar^2}{2m_j} \sum_{\alpha=1}^{d} \pndv{x_{j\alpha}}{}{2} + V\left( \left\{ x_{j\alpha} \right\} \right) \right]$$

となります。ここで、$m_j$は$j$番目の粒子(言うまでもないですが、この文脈では物理的な粒子です)の質量で、$x_{j\alpha}$は$j$番目の粒子の$\alpha$方向における位置です。この場合、次のように座標をスケールします。

$$ \tag{28} x_{j\alpha}=\sqrt{ \frac{m}{m_j} } x'_{j\alpha}$$

ここで、$m>0$は任意に選ぶことができますです。このように選ぶことによって、あたかも$d'=S\times d$次元空間を移動する質量$m$の$1$粒子系とみなすことができ、式(26)の場合と同じように扱うことができます。

5. 無次元化

通常シミュレーションを実行する際には、数値的に扱いやすくなるように方程式を無次元化します。式(4)のシュレディンガー方程式において、次のように時間$\tau$、空間$x$を無次元化します。

$$\tag{29} \tau\rightarrow \tau T,~~~~x\rightarrow x L$$

この置き換えによって式(6)は

$$\tag{30} \pdv{\tau}{\Psi}= \frac{\hbar T}{2mL^2} \pndv{x}{\Psi}{2}-\frac{T \mathcal{E}}{\hbar} \left[V(x)-E_R \right]\Psi$$

となります。注意すべきは、新しい変数$"x,t"$が無次元量となり、代わりに$L,T$にそれぞれ長さと時間の次元を持たせているということです。ここで、$x\rightarrow xL$の置き換えによって$V(x)$が$V(x)\rightarrow \mathcal{E}V(x)$と変換されるとしました。また、$L$および$T,\mathcal{E}$を次の関係を満たすように選ぶと便利です。

$$\tag{31} \begin{aligned} \frac{\hbar T}{2mL^2}&=\frac{1}{2}\\ \frac{T\mathcal{E}}{\hbar}&=1 \end{aligned}$$

したがって式(30)は

$$\tag{32} \pdv{\tau}{\Psi}= \frac{1}{2} \pndv{x}{\Psi}{2}- \left[V(x)-E_R \right]\Psi$$

となります。式(31)より、$L,T,\mathcal{E}$のうち1つを決めると残りの2つは自動的に決まります。この量は次元解析によって決めるとよいです。
 この無次元化によって、式(10)と(11)は

$$ \tag{33} P(x_n,x_{n-1})=\left( \frac{1}{2\pi \Delta \tau} \right)^{1/2} \exp{\left\{ -\frac{(x_n-x_{n-1})^2}{2\Delta \tau} \right\}}$$$$ \tag{34} W(x_n)=\exp{ \left\{ -\Delta \tau(V(x_n)-E_R) \right\} }$$

と変更されます(論文の方では、式(34)には$\Delta \tau$が入っていませんが、おそらく誤植だと思います)。また、式(33)から

$$ \tag{35} \sigma = \sqrt{\Delta \tau}$$

となることがわかります。

6. 数値計算の例

それでは実際に数値計算を行っていきます。ソースコードは以下に置きました。
https://github.com/RoofOZI/DiffusionQuantumMonteCarlo.git
以下の数値計算では、"粒子"の数を$2000$個、$\Delta \tau=0.05$とし、ランダムウォークを$2000$回行いました($\tau=100$に相当)。また、式(23)における$\alpha$の値は$1$としました。エネルギー固有値は、$\tau=50\sim100$ に渡る $E_\mathrm{R}$の平均値で評価しています。

1次元調和振動子

1次元調和振動子の系のポテンシャルエネルギーは

$$ \tag{36} V(x)=\frac{1}{2}m\omega^2 x^2$$

と書き表されます。$m,\omega$はそれぞれ粒子の質量、系を特徴付ける振動数です。次元解析から、$L=\sqrt{\hbar/m\omega}$とします。このとき、$T=1/\omega,~\mathcal{E}=\hbar \omega$です。そうすると、式(36)のポテンシャルエネルギーは単純に

$$ \tag{37} V(x)=\frac{1}{2}x^2$$

となります。エネルギー固有値の理論値は$E_0=0.5$ です。一方、数値計算の結果では、 $E_0=0.50085(37)$という理論値に近い値が得られました。数値計算の際、"粒子"の初期位置は$x=0$としました。(また、上記リンクに得られた波動関数の図を載せています。)


水素原子

水素原子の場合におけるポテンシャルエネルギーはクーロンポテンシャルであり、

$$\tag{38} V(x,y,z)=-\alpha\frac{\hbar c}{\sqrt{x^2+y^2+z^2}}$$

と表されます。$\alpha\approx1/137$は微細構造定数です。この時次元解析から、$L=\alpha^{-1}\frac{\hbar}{mc}(\approx 0.53\times 10^{-10} [\mathrm{m}]),$ $T=\alpha^{-2}\frac{\hbar}{mc^2},$ $E=\alpha^2mc^2(\approx27.21~[\mathrm{eV}])$となります。無次元化により、

$$\tag{39} V(x,y,z)=-\frac{1}{\sqrt{x^2+y^2+z^2}}$$

となります。粒子の初期位置を$(x,y,z)=(0,0,1)$として計算した結果、エネルギー固有値として$E_0=-0.50082(54)$が得られました。理論値は$E_0=-0.5$なので、非常に良く一致していることが分かります。


ヘリウム原子

ヘリウム原子の場合のハミルトニアンを書き下してみます。

$$\tag{40} \hat{H}= -\frac{\hbar}{2m} \left( \nabla_1^2 + \nabla_2^2 \right) +\alpha\hbar c \left\{ -\frac{2}{r_1^2} - \frac{2}{r_2^2} + \frac{1}{r_{12}^2} \right\}$$

ヘリウム原子は$2$電子系であり、それぞれの位置座標を$\boldsymbol{r}_1,$ $\boldsymbol{r}_2$とし、$r_{1}=\|\boldsymbol{r}_1\|,$ $r_{2}=\|\boldsymbol{r}_2\|,$ $r_{12}=\|\boldsymbol{r}_1 - \boldsymbol{r}_2\|$としました。無次元化は水素原子の場合と同様にでき、

$$\tag{41} \hat{H}= -\frac{1}{2} \left( \nabla_1^2 + \nabla_2^2 \right) + \left\{ -\frac{2}{r_1^2} - \frac{2}{r_2^2} + \frac{1}{r_{12}^2} \right\}$$

となります。
 ヘリウム原子の基底状態は、スピン部分が$1$重項状態で反対称、空間部分が対称となっており、全体で反対称となっています。ハミルトニアンにスピンに関わる項がないので、空間部分の波動関数$\phi(\boldsymbol{r}_1,\boldsymbol{r}_2)$ はそのまま$\hat{H}\phi(\boldsymbol{r}_1,\boldsymbol{r}_2)=E \phi(\boldsymbol{r}_1,\boldsymbol{r}_2)$ を満たします。これを踏まえて、ポテンシャルエネルギーはそのまま、

$$ \tag{42} V(\boldsymbol{r}_1,\boldsymbol{r}_2)=-\frac{2}{r_1^2} - \frac{2}{r_2^2} + \frac{1}{r_{12}^2}$$

となります。
 粒子の初期位置を$\boldsymbol{r}_1=(0,0,1),$ $\boldsymbol{r}_2=(1,0,0)$とし、今回は$\Delta \tau=0.01$で数値計算を行った結果、$E_0=-2.8932(27)$ $(\approx -78.72 ~[\mathrm{eV}])$ という値が得られました。実験値は$E_0=-2.90338583(13)$ (下記URLより) なので、良く一致しています。
(https://en.wikipedia.org/wiki/Helium_atom)

課題

エネルギー固有値は非常によく理論値/実験値と一致していましたが、もちろん"粒子の個数"、初期位置や$\Delta \tau$など、パラメータの与え方によって得られる値は少なからず変動します。このような誤差の評価をどうするかは難しい問題です。
 また、$1$次元のときの"粒子"の分布は波動関数とよく一致していたのですが、$3$次元の場合、例えば動径方向の分布は形こそ似ているものの、全体的に引き延ばされたような分布になっている、というような結果になりました。多次元のときの粒子の分布の取り扱いは見直す必要がありそうです。
 モンテカルロ法で数値計算を行う際、符号問題という問題によく直面するそうです。今回はこれには触れませんでしたが、数値計算の分野では非常に有名らしいので、何か今後数値計算をする人は(私も含めて)知っておくべき問題だと思います。