MiyanTarumi’s blog

MiyanTarumi’s blog

感想とか色々。数物系は書くか微妙。。。

変分モンテカルロ(VMC)をやってみる~ヘリウム原子~

Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。

www.maruzen-publishing.co.jp

ヘリウム原子

モデル

最後にヘリウム原子についても同じくVMCをやってみよう。
この系のHamiltonianはBorn-Oppenheimer近似の下で\textbf{r}_1を電子1の座標、\textbf{r}_2を電子2の座標とすると、

H=\displaystyle-\frac{1}{2}\nabla_1^2-\frac{1}{2}\nabla_2^2-\frac{2}{r_1}-\frac{2}{r_2}+\frac{1}{|\textbf{r}_1-\textbf{r}_2|}

と書ける。これまでの時と同様に試行関数を与えたいのだが、この場合は良い試行関数を作ること自体が問題になってしまうので、ここでは天下り的であるが試行関数として次のようなものを使用することにしよう。

\displaystyle \psi(\textbf{r}_1,\textbf{r}_2)=e^{-2r_1}e^{-2r_2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

ここでr_{12}=|\textbf{r}_1-\textbf{r}_2|とし、これまでと同様に\alphaは変分パラメータである。
次にこの試行関数を用いた時に局所エネルギーE_Lがどのようになるのかを計算してみる*1。 微分演算子の入っていない部分はそのまま係数としてかかるだけなので微分演算子の部分だけを計算すればいいだろう。そこでまず\nabla_1^2について計算することにしよう。合成関数のラプラシアンに対する作用が、

\nabla^2\left(fg\right)=\nabla\cdot\nabla\left(fg\right))=\nabla\cdot\left(\left(\nabla f\right)g+f\nabla g\right)\\ =\left(\nabla\cdot\nabla f\right)g+\nabla f\cdot\nabla g+\nabla f\cdot\nabla g+f\nabla\cdot\nabla g\\ =\nabla^2 f+2\nabla f\cdot\nabla g+f\nabla^2 g

であることに注意すると、\textbf{r}_1に関係する部分の計算は、

\nabla_1^2\left( e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right) =\left(\nabla_1^2 e^{-2r_{1}}\right)e^{\frac{r_{12}}{2(1+\alpha r_{12})}}+\nabla_1 e^{-2r_1}\cdot\nabla_1e^{\frac{r_{12}}{2(1+\alpha r_{12})}}+e^{-2r_1}\nabla_1^2e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\

を計算すればよいことになる。まず初項について計算しよう。この項の計算は簡単で1体の時のラプラシアンの計算と全く同じであるから、

\nabla_1^2 e^{-2r_{1}} =\displaystyle\left[\frac{\partial^2}{\partial r_1^2}+\frac{2}{r_1}\frac{\partial}{\partial r_1}\right]e^{-2r_1}\\ =\displaystyle\left[4-\frac{4}{r_1}\right]e^{-2r_1}

より、

\left(\nabla_1^2 e^{-2r_{1}}\right)e^{\frac{r_{12}}{2(1+\alpha r_{12})}} =\displaystyle\left[4-\frac{4}{r_1}\right]e^{-2r_1}e^{-2r_2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

となる。
続いて第二項について計算する。ここでこれ以降の計算の表記を簡単にするためにr_{12}が入っている部分の指数の肩を

U_{12}=\displaystyle\frac{r_{12}}{2(1+\alpha r_{12})}

と書くことにし、その微分を

DU_{12}=\displaystyle\frac{\partial}{\partial r_{12}}U_{12}=-\frac{1}{2}\frac{\alpha r_{12}}{(1+\alpha r_{12})^2}+\frac{1}{2}\frac{1}{1+\alpha r_{12}}

と書くことにする。すると第二項の計算は、

\nabla_1 e^{-2r_1}\cdot\nabla_1e^{\frac{r_{12}}{2(1+\alpha r_{12})}} =\displaystyle\sum_{i}\left(\frac{\partial}{\partial x_{1i}}e^{-2r_1}\right)\left(\frac{\partial }{\partial x_{1i}}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\sum_{i}\left(\frac{\partial r_{1}}{\partial x_{1i}}\frac{\partial}{\partial r_1}e^{-2r_1}\right)\left(\frac{\partial r_{12}}{\partial x_{1i}}\frac{\partial }{\partial r_{12}}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\sum_{i}\left[(-2)\cdot DU_{12}\right]\frac{\partial r_1}{\partial x_{1i}}\frac{\partial r_{12}}{\partial x_{1i}}e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle \sum_{i}\left[(-2)\cdot DU_{12}\right]\frac{x_{1i}}{r_1}\frac{x_{1i}-x_{2i}}{r_{12}}e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle\frac{-2DU_{12}}{r_{12}}\frac{r_1^2-\textbf{r}_1\cdot\textbf{r}_2}{r_1}e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

となる。ここでx_{1i}は電子1の座標の成分を、x_{2i}は電子2の座標の成分をそれぞれ表すとした。
最後に第三項の計算を行おう。この項ではr_{12}の入った項の微分を二回行う必要があるのでU_{12}の二回微分を表すものとして

DDU_{12}=\displaystyle\frac{\partial^2}{\partial r_{12}^2}U_{12}=\frac{\partial}{\partial r_{12}}\left(-\frac{1}{2}\frac{\alpha r_{12}}{(1+\alpha r_{12})^2}+\frac{1}{2}\frac{1}{1+\alpha r_{12}}\right)\\ =\displaystyle\frac{\alpha^2 r_{12}}{(1+\alpha r_{12})^3}-\frac{\alpha }{(1+\alpha r_{12})^2}

とDDU_{12}を定義しよう。すると第三項は、

\nabla_1^2e^{\frac{r_{12}}{2(1+\alpha r_{12})}} =\displaystyle \sum_{i}\frac{\partial^2}{\partial x_{1i}^2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}} =\displaystyle\sum_{i}\frac{\partial }{\partial x_{1i}}\left(DU_{12}\frac{x_{1i}-x_{2i}}{r_{12}}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\sum_{i}\left[\frac{DU_{12}}{r_{12}}+\left(\frac{DDU_{12}}{r_{12}}-\frac{DU_{12}}{r_{12}^2}+\frac{(DU_{12})^2}{r_{12}}\right)(x_{1i}-x_{2i})\frac{\partial r_{12}}{\partial x_{1i}}\right]e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle\sum_{i}\left[\frac{DU_{12}}{r_{12}}+\left(\frac{DDU_{12}}{r_{12}}-\frac{DU_{12}}{r_{12}^2}+\frac{(DU_{12})^2}{r_{12}}\right)\frac{(x_{1i}-x_{2i})^2}{r_{12}}\right]e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle\left[\frac{3DU_{12}}{r_{12}}+\left(\frac{DDU_{12}}{r_{12}}-\frac{DU_{12}}{r_{12}^2}+\frac{(DU_{12})^2}{r_{12}}\right)r_{12}\right]e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\ =\displaystyle\left[\frac{2DU_{12}}{r_{12}}+DDU_{12}+(DU_{12})^2\right]e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\\

となる。以上から

\nabla_1^2\left( e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\left[4-\frac{4}{r_1}-\frac{4DU_{12}}{r_{12}}\frac{r_1^2-\textbf{r}_1\cdot\textbf{r}_2}{r_1}+\frac{2DU_{12}}{r_{12}}+DDU_{12}+(DU_{12})^2\right]e^{-2r_1}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

となることがわかった。同様にして\nabla_2^2について計算すると、

\nabla_2^2\left( e^{-2r_2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}\right)\\ =\displaystyle\left[4-\frac{4}{r_2}+\frac{4DU_{12}}{r_{12}}\frac{\textbf{r}_1\cdot\textbf{r}_2-r_2^2}{r_2}+\frac{2DU_{12}}{r_{12}}+DDU_{12}+(DU_{12})^2\right]e^{-2r_2}e^{\frac{r_{12}}{2(1+\alpha r_{12})}}

を得ることができる*2。
以上から

H\psi(\textbf{r}_1,\textbf{r}_2)\\ =\displaystyle\left[-4-\frac{1}{2}\left(-\frac{4DU_{12}}{r_{12}}\frac{r_1^2-\textbf{r}_1\cdot\textbf{r}_2}{r_1}+\frac{4DU_{12}}{r_{12}}\frac{\textbf{r}_1\cdot\textbf{r}_2-r_2^2}{r_2}\right.\\ \left.+\frac{4DU_{12}}{r_{12}}+2DDU_{12}+2(DU_{12})^2\right)+\frac{1}{r_{12}}\right]\psi(\textbf{r}_1,\textbf{r}_2)

と計算でき、よって局所エネルギーE_Lは、

E_L=\displaystyle-4-\frac{1}{2}\left(-\frac{4DU_{12}}{r_{12}}\frac{r_1^2-\textbf{r}_1\cdot\textbf{r}_2}{r_1}+\frac{4DU_{12}}{r_{12}}\frac{\textbf{r}_1\cdot\textbf{r}_2-r_2^2}{r_2}\right.\\ \left.\displaystyle+\frac{4DU_{12}}{r_{12}}+2DDU_{12}+2(DU_{12})^2\right)+\displaystyle\frac{1}{r_{12}}

と求めることができた。

VMCの方法とソースコード

よってプログラミングで行う操作は次の通りである。

調和振動子のVMC

**REPAT**
始めの変分パラメータ\alphaを決める。 ランダムな場所にN個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**

1. k番目のサンプリング点をX^{(k)}とする。walkerの移動後の座標としてk+1番目の候補をX^{k}を中心とする正規分布に従ってX'=X^{(k)}+\Delta Xとして提案する。
2. メトロポリステスト:
[0,1]の間の一様乱数rを生成し、r< |\psi(X')/\psi(X)|^2であれば提案を受理してX^{(k+1)}=X'とする。そうでなければ候補を棄却してX^{(k+1)}=X^{(k)}とする。
3. サンプリング点での局所エネルギーE_Lを計算する。

**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値を求める。
Newton法に従って変分パラメータ\alphaを更新する。
**END**

ここでXは座標であり、今の場合X=(x_1,y_1,z_1, x_2, y_2, z_2)である。
このプログラムを書いてみよう。

▶C言語のプログラム（クリックすると展開されます）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include "MT.h"

int shokihaichi(double data1[], double data2[], double data3[]);/*プロトタイプ宣言*/
int metropolis_walk(double data1[], double data2[], double data3[], double data4[], double data5[], double data6[]);/*プロトタイプ宣言*/
int shokika(double data1[], double data2[], double data3[]);/*プロトタイプ宣言*/

double sigma=0.15;

int main(void)
{
    int sample=50000;
    int i, j;
    int acc[500];

    double alpha_0=0.05;
    double alpha_ini;
    double alpha_fin;
    double x1[500], y1[500], z1[500];
    double x2[500], y2[500], z2[500];
    double x1_backup[500], y1_backup[500], z1_backup[500];
    double x2_backup[500], y2_backup[500], z2_backup[500];
    double r1_ini[500], r2_ini[500], r12_ini[500];
    double r1_fin[500], r2_fin[500], r12_fin[500];
    double r1[500], r2[500], r12[500];
    double psi_ini[500], psi_fin[500];
    double Denom[500];/*エネルギー計算に出てくる分母の関数*/
    double U12[500], DU12[500], DDU12[500];
    double R1DR12[500], R2DR12[500];
    double Sum_E_w[500]={0};
    double E_w[500]={0};
    double E=0;
    double Sum_Edlnpsi_w[500]={0};
    double Edlnpsi_w[500]={0};
    double Edlnpsi=0;
    double Sum_dlnpsi_w[500]={0};
    double dlnpsi_w[500]={0};
    double dlnpsi=0;
    double dE;
    double h=1e-3;/*数値微分の幅*/

    init_genrand((unsigned)time(NULL));

    alpha_fin=alpha_0;

    do
    {
        alpha_ini=alpha_fin;/*前のループでの結果を始まりのalphaに代入*/

        /************************************************/
        /*dEを求める.*************************************/
        /*alpha_iniでの<E>, <Edlnpsi>, <dlnpsi>を計算する.*/
        /************************************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*Sumの配列を初期化*/

        shokihaichi(x1, y1, z1);/*ランダムにwalkrを配置*/
        shokihaichi(x2, y2, z2);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x1, y1, z1, x1_backup, y1_backup, z1_backup);/*提案分布に従ってwalkerを変化させる*/
            metropolis_walk(x2, y2, z2, x2_backup, y2_backup, z2_backup);

            /*Metropolis Testを行う*/

            for ( j = 0; j < 500; j++)
            {
                r1_ini[j]=sqrt(x1_backup[j]*x1_backup[j]+y1_backup[j]*y1_backup[j]+z1_backup[j]*z1_backup[j]);
                r2_ini[j]=sqrt(x2_backup[j]*x2_backup[j]+y2_backup[j]*y2_backup[j]+z2_backup[j]*z2_backup[j]);
                r12_ini[j]=sqrt((x1_backup[j]-x2_backup[j])*(x1_backup[j]-x2_backup[j])+(y1_backup[j]-y2_backup[j])*(y1_backup[j]-y2_backup[j])+(z1_backup[j]-z2_backup[j])*(z1_backup[j]-z2_backup[j]));

                psi_ini[j]=exp(-4*(r1_ini[j]+r2_ini[j])+r12_ini[j]/(1+alpha_ini*r12_ini[j]));/*始めの確率分布*/

                r1_fin[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);
                r2_fin[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                r12_fin[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));

                psi_fin[j]=exp(-4*(r1_fin[j]+r2_fin[j])+r12_fin[j]/(1+alpha_ini*r12_fin[j]));/*候補の確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x1[j]=x1_backup[j];
                    y1[j]=y1_backup[j];
                    z1[j]=z1_backup[j];
                    x2[j]=x2_backup[j];
                    y2[j]=y2_backup[j];
                    z2[j]=z2_backup[j];
                }

                /*Metropolis Testの終了*/
                /*walkerについての<E>, <Edlnpsi>, <dlnpsi>を計算する.*/

                if (i>=4000)
                {
                    r1[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);/*更新し終えた位置でパラメタを計算*/
                    r2[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                    r12[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));
                    R1DR12[j]=x1[j]*(x1[j]-x2[j])+y1[j]*(y1[j]-y2[j])+z1[j]*(z1[j]-z2[j]);
                    R2DR12[j]=x2[j]*(x1[j]-x2[j])+y2[j]*(y1[j]-y2[j])+z2[j]*(z1[j]-z2[j]);
                    Denom[j]=1/(1+alpha_ini*r12[j]);
                    U12[j]=0.5*r12[j]*Denom[j];
                    DU12[j]=-0.5*alpha_ini*r12[j]*Denom[j]*Denom[j]+0.5*Denom[j];
                    DDU12[j]=alpha_ini*alpha_ini*r12[j]*Denom[j]*Denom[j]*Denom[j]-alpha_ini*Denom[j]*Denom[j];

                    Sum_E_w[j]=Sum_E_w[j]-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j];
                    E_w[j]=Sum_E_w[j]/(i-3999);

                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]))*(-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j]);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);

                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);

                }
                
            }
            
        }
        
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 500; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }

        E=E/500;
        Edlnpsi=Edlnpsi/500;
        dlnpsi=dlnpsi/500;/*dE/dalphaの計算をするための部品を計算した.*/

        dE=2*(Edlnpsi-E*dlnpsi);/*dE/dalphaの計算ができた.*/

        /*****************************/
        /*ddE/ddalphaの数値微分を行う.*/
        /****************************/

        double alpha_plus=alpha_ini+h;/*微分のためにalphaを変化させる*/
        double alpha_minus=alpha_ini-h;

        /***********************************/
        /***********************************/
        /*alpha_plusでのdE/dalphaを計算する.*/
        /***********************************/
        /***********************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*Sumの配列を初期化*/

        shokihaichi(x1, y1, z1);/*ランダムにwalkrを配置*/
        shokihaichi(x2, y2, z2);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x1, y1, z1, x1_backup, y1_backup, z1_backup);/*提案分布に従ってwalkerを変化させる*/
            metropolis_walk(x2, y2, z2, x2_backup, y2_backup, z2_backup);

            /*Metropolis Testを行う*/

            for ( j = 0; j < 500; j++)
            {
                r1_ini[j]=sqrt(x1_backup[j]*x1_backup[j]+y1_backup[j]*y1_backup[j]+z1_backup[j]*z1_backup[j]);
                r2_ini[j]=sqrt(x2_backup[j]*x2_backup[j]+y2_backup[j]*y2_backup[j]+z2_backup[j]*z2_backup[j]);
                r12_ini[j]=sqrt((x1_backup[j]-x2_backup[j])*(x1_backup[j]-x2_backup[j])+(y1_backup[j]-y2_backup[j])*(y1_backup[j]-y2_backup[j])+(z1_backup[j]-z2_backup[j])*(z1_backup[j]-z2_backup[j]));

                psi_ini[j]=exp(-4*(r1_ini[j]+r2_ini[j])+r12_ini[j]/(1+alpha_plus*r12_ini[j]));/*始めの確率分布*/

                r1_fin[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);
                r2_fin[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                r12_fin[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));

                psi_fin[j]=exp(-4*(r1_fin[j]+r2_fin[j])+r12_fin[j]/(1+alpha_plus*r12_fin[j]));/*候補の確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x1[j]=x1_backup[j];
                    y1[j]=y1_backup[j];
                    z1[j]=z1_backup[j];
                    x2[j]=x2_backup[j];
                    y2[j]=y2_backup[j];
                    z2[j]=z2_backup[j];
                }

                /*Metropolis Testの終了*/
                /*walkerについての<E>, <Edlnpsi>, <dlnpsi>を計算する.*/

                if (i>=4000)
                {
                    r1[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);/*更新し終えた位置でパラメタを計算*/
                    r2[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                    r12[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));
                    R1DR12[j]=x1[j]*(x1[j]-x2[j])+y1[j]*(y1[j]-y2[j])+z1[j]*(z1[j]-z2[j]);
                    R2DR12[j]=x2[j]*(x1[j]-x2[j])+y2[j]*(y1[j]-y2[j])+z2[j]*(z1[j]-z2[j]);
                    Denom[j]=1/(1+alpha_plus*r12[j]);
                    U12[j]=0.5*r12[j]*Denom[j];
                    DU12[j]=-0.5*alpha_plus*r12[j]*Denom[j]*Denom[j]+0.5*Denom[j];
                    DDU12[j]=alpha_plus*alpha_plus*r12[j]*Denom[j]*Denom[j]*Denom[j]-alpha_plus*Denom[j]*Denom[j];

                    Sum_E_w[j]=Sum_E_w[j]-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j];
                    E_w[j]=Sum_E_w[j]/(i-3999);

                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]))*(-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j]);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);

                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);

                }
                
            }
            
        }
        
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 500; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }

        E=E/500;
        Edlnpsi=Edlnpsi/500;
        dlnpsi=dlnpsi/500;/*dE/dalpha_plusの計算をするための部品を計算した.*/

        double dE_plus=2*(Edlnpsi-E*dlnpsi);/*dE/dalpha_plusの計算ができた.*/

        /***********************************/
        /***********************************/
        /*alpha_minusでのdE/dalphaを計算する.*/
        /***********************************/
        /***********************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*Sumの配列を初期化*/

        shokihaichi(x1, y1, z1);/*ランダムにwalkrを配置*/
        shokihaichi(x2, y2, z2);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x1, y1, z1, x1_backup, y1_backup, z1_backup);/*提案分布に従ってwalkerを変化させる*/
            metropolis_walk(x2, y2, z2, x2_backup, y2_backup, z2_backup);

            /*Metropolis Testを行う*/

            for ( j = 0; j < 500; j++)
            {
                r1_ini[j]=sqrt(x1_backup[j]*x1_backup[j]+y1_backup[j]*y1_backup[j]+z1_backup[j]*z1_backup[j]);
                r2_ini[j]=sqrt(x2_backup[j]*x2_backup[j]+y2_backup[j]*y2_backup[j]+z2_backup[j]*z2_backup[j]);
                r12_ini[j]=sqrt((x1_backup[j]-x2_backup[j])*(x1_backup[j]-x2_backup[j])+(y1_backup[j]-y2_backup[j])*(y1_backup[j]-y2_backup[j])+(z1_backup[j]-z2_backup[j])*(z1_backup[j]-z2_backup[j]));

                psi_ini[j]=exp(-4*(r1_ini[j]+r2_ini[j])+r12_ini[j]/(1+alpha_minus*r12_ini[j]));/*始めの確率分布*/

                r1_fin[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);
                r2_fin[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                r12_fin[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));

                psi_fin[j]=exp(-4*(r1_fin[j]+r2_fin[j])+r12_fin[j]/(1+alpha_minus*r12_fin[j]));/*候補の確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x1[j]=x1_backup[j];
                    y1[j]=y1_backup[j];
                    z1[j]=z1_backup[j];
                    x2[j]=x2_backup[j];
                    y2[j]=y2_backup[j];
                    z2[j]=z2_backup[j];
                }

                /*Metropolis Testの終了*/
                /*walkerについての<E>, <Edlnpsi>, <dlnpsi>を計算する.*/

                if (i>=4000)
                {
                    r1[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);/*更新し終えた位置でパラメタを計算*/
                    r2[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                    r12[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));
                    R1DR12[j]=x1[j]*(x1[j]-x2[j])+y1[j]*(y1[j]-y2[j])+z1[j]*(z1[j]-z2[j]);
                    R2DR12[j]=x2[j]*(x1[j]-x2[j])+y2[j]*(y1[j]-y2[j])+z2[j]*(z1[j]-z2[j]);
                    Denom[j]=1/(1+alpha_minus*r12[j]);
                    U12[j]=0.5*r12[j]*Denom[j];
                    DU12[j]=-0.5*alpha_minus*r12[j]*Denom[j]*Denom[j]+0.5*Denom[j];
                    DDU12[j]=alpha_minus*alpha_minus*r12[j]*Denom[j]*Denom[j]*Denom[j]-alpha_minus*Denom[j]*Denom[j];

                    Sum_E_w[j]=Sum_E_w[j]-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j];
                    E_w[j]=Sum_E_w[j]/(i-3999);

                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]))*(-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j]);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);

                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-0.5*(r12[j]*Denom[j])*(r12[j]*Denom[j]);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);

                }
                
            }
            
        }
        
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 500; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }

        E=E/500;
        Edlnpsi=Edlnpsi/500;
        dlnpsi=dlnpsi/500;/*dE/dalpha_minusの計算をするための部品を計算した.*/

        double dE_minus=2*(Edlnpsi-E*dlnpsi);/*dE/dalpha_minusの計算ができた.*/

        /***************************************************/
        /***************************************************/
        /*alphaをNewton法で求めるための部品をすべて計算し終えた*/
        /***************************************************/
        /***************************************************/

        /*alphaの更新を行う*/

        double ddE=(dE_plus-dE_minus)/(2*h);/*dE/dalphaの数値微分ddEが計算できた.*/
        
        alpha_fin=alpha_ini-dE/ddE;/*alphaを更新する.*/
        printf("alpha_fin=%.10f, alpha_ini=%.10f dE=%.10f\n", alpha_fin, alpha_ini, dE);
    } while (fabs(alpha_fin-alpha_ini)>1e-4);

    double alpha_val=alpha_fin;/*最適化した変分パラメータ*/
    printf("変分パラメータはalpha=%.10f", alpha_val);

    /*****************************/
    /*****************************/
    /*alphaでのエネルギーを計算する*/
    /*****************************/
    /*****************************/

    shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*Sumの配列を初期化*/

    shokihaichi(x1, y1, z1);/*ランダムにwalkrを配置*/
    shokihaichi(x2, y2, z2);

    for ( i = 0; i < sample; i++)
    {
        metropolis_walk(x1, y1, z1, x1_backup, y1_backup, z1_backup);/*提案分布に従ってwalkerを変化させる*/
        metropolis_walk(x2, y2, z2, x2_backup, y2_backup, z2_backup);

        /*Metropolis Testを行う*/

        for ( j = 0; j < 500; j++)
        {
            r1_ini[j]=sqrt(x1_backup[j]*x1_backup[j]+y1_backup[j]*y1_backup[j]+z1_backup[j]*z1_backup[j]);
            r2_ini[j]=sqrt(x2_backup[j]*x2_backup[j]+y2_backup[j]*y2_backup[j]+z2_backup[j]*z2_backup[j]);
            r12_ini[j]=sqrt((x1_backup[j]-x2_backup[j])*(x1_backup[j]-x2_backup[j])+(y1_backup[j]-y2_backup[j])*(y1_backup[j]-y2_backup[j])+(z1_backup[j]-z2_backup[j])*(z1_backup[j]-z2_backup[j]));

            psi_ini[j]=exp(-4*(r1_ini[j]+r2_ini[j])+r12_ini[j]/(1+alpha_val*r12_ini[j]));/*始めの確率分布*/

            r1_fin[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);
            r2_fin[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
            r12_fin[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));

            psi_fin[j]=exp(-4*(r1_fin[j]+r2_fin[j])+r12_fin[j]/(1+alpha_val*r12_fin[j]));/*候補の確率分布*/

            double metropolis=genrand_real1();

            if (psi_fin[j]/psi_ini[j]>metropolis)
            {
                acc[j]=acc[j]+1;
            }
            else
            {
                x1[j]=x1_backup[j];
                y1[j]=y1_backup[j];
                z1[j]=z1_backup[j];
                x2[j]=x2_backup[j];
                y2[j]=y2_backup[j];
                z2[j]=z2_backup[j];
            }

            /*Metropolis Testの終了*/
            /*walkerについての<E>, <Edlnpsi>, <dlnpsi>を計算する.*/

            if (i>=4000)
            {
                r1[j]=sqrt(x1[j]*x1[j]+y1[j]*y1[j]+z1[j]*z1[j]);/*更新し終えた位置でパラメタを計算*/
                r2[j]=sqrt(x2[j]*x2[j]+y2[j]*y2[j]+z2[j]*z2[j]);
                r12[j]=sqrt((x1[j]-x2[j])*(x1[j]-x2[j])+(y1[j]-y2[j])*(y1[j]-y2[j])+(z1[j]-z2[j])*(z1[j]-z2[j]));
                R1DR12[j]=x1[j]*(x1[j]-x2[j])+y1[j]*(y1[j]-y2[j])+z1[j]*(z1[j]-z2[j]);
                R2DR12[j]=x2[j]*(x1[j]-x2[j])+y2[j]*(y1[j]-y2[j])+z2[j]*(z1[j]-z2[j]);
                Denom[j]=1/(1+alpha_val*r12[j]);
                U12[j]=0.5*r12[j]*Denom[j];
                DU12[j]=-0.5*alpha_val*r12[j]*Denom[j]*Denom[j]+0.5*Denom[j];
                DDU12[j]=alpha_val*alpha_val*r12[j]*Denom[j]*Denom[j]*Denom[j]-alpha_val*Denom[j]*Denom[j];

                Sum_E_w[j]=Sum_E_w[j]-4-0.5*(2*(-2)/r1[j]*R1DR12[j]*DU12[j]/r12[j]-2*(-2)/r2[j]*R2DR12[j]*DU12[j]/r12[j]+4*DU12[j]/r12[j]+2*DU12[j]*DU12[j]+2*DDU12[j])+1/r12[j];
                E_w[j]=Sum_E_w[j]/(i-3999);
            }
            
        }
        
    }
    
    E=0;/*変数の初期化*/
    /*walker全ての平均を取ることにする*/
    for ( j = 0; j < 500; j++)
    {
        E=E+E_w[j];
    }

    E=E/500;
    
    printf("基底状態のエネルギーは%.10f", E);
    return 0;
}

int shokihaichi(double data1[], double data2[], double data3[])
{
    int i;
    for ( i = 0; i < 500; i++)
    {
        data1[i]=(genrand_real1()-0.5)*2*2;
        data2[i]=(genrand_real1()-0.5)*2*2;
        data3[i]=(genrand_real1()-0.5)*2*2;
    }
    return 0;
}

int metropolis_walk(double data1[], double data2[], double data3[], double data4[], double data5[], double data6[])/*1~3がx,y,z, 4~6はそれらのバックアップ*/
{
    int i;
    double p[500];
    double q[500];
    double pi=asin(1e0)*2e0;
    for ( i = 0; i < 500; i++)
    {
        /*x*/
        data4[i]=data1[i];/*バックアップにxの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data1[i]=data1[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);

        /*y*/
        data5[i]=data2[i];/*バックアップにyの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data2[i]=data2[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);

        /*z*/
        data6[i]=data3[i];/*バックアップにxの値を保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data3[i]=data3[i]+sigma*sqrt(-2*log(p[i]))*cos(2*pi*q[i]);
    }
    return 0;
}

int shokika(double data1[], double data2[], double data3[])
{
    int i;
    for ( i = 0; i < 500; i++)
    {
        data1[i]=0;
        data2[i]=0;
        data3[i]=0;
    }
    return 0;
}

このプログラムを実行すると、\alpha=0.1433となり、その時の基底状態のエネルギーとして-2.8783が得られる。これはHartree-Fockの値-2.8617やDFTの値-2.83、厳密値-2.9037と比肩しうる数値になっている。

*1:Thijssenの教科書に載っている式は和訳も原著(2nd Edition)もこの式が間違っているので教科書の式でプログラムを組むと間違った結果が出てしまうので注意。付録のソースコードでは正しい式になってるのにどうして。。。

*2:要はr_{12}のx_{2i}での微分一回につき負号がでるので第二項のみ係数が変わる。