変分モンテカルロ(VMC)をやってみる~水素原子~

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

www.maruzen-publishing.co.jp

水素原子
- モデル
- VMCの方法とソースコード

水素原子

モデル

調和振動子に続いて水素原子についても同じくVMCをやってみよう。水素原子の問題は1次元の問題に帰着できるが、ここでは3次元の問題として扱うことにする。
この系のHamiltonianはBorn-Oppenheimer近似の下で、

$H=\displaystyle-\frac{1}{2}\nabla^2-\frac{1}{r}$

と書ける。水素原子の基底状態の厳密解は $e^{-r}$ であるから、ここでは試行関数として $\psi(r)=e^{-\alpha r}$ を選ぶことにしよう。ここで $\alpha$ は変分パラメータである。このように試行関数を決めた時に局所エネルギーはどうなるかというと、

$H\psi(r) =\displaystyle\left[-\frac{1}{2}\nabla^2-\frac{1}{r}\right]e^{-\alpha r}\\ =\displaystyle\left[-\frac{1}{2}\left(\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}\right)-\frac{1}{r}\right]e^{-\alpha r}\\ =\displaystyle\left[-\frac{1}{r}-\frac{1}{2}\alpha\left(\alpha-\frac{2}{r}\right)\right]e^{-\alpha r}$

より、

$\displaystyle E_L=-\frac{1}{r}-\frac{1}{2}\alpha\left(\alpha-\frac{2}{r}\right)$

となる。

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=\exp(-2\alpha(r'-r))$ であれば提案を受理して $X^{(k+1)}=X'$ とする。そうでなければ候補を棄却して $X^{(k+1)}=X^{(k)}$ とする。
3. サンプリング点での局所エネルギー $E_L$ を計算する。

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

ここで $X$ は座標であり、今の場合 $X=(x,y,z)$ である。
このプログラムを書いてみよう。

▶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.4;/*metropoliswalkでの正規乱数の分散*/

int main(void)
{
    int sample=30000;
    int i, j;
    int acc[400]={0};

    double alpha_0=0.8;/*変分を始める初期値*/
    double alpha_ini;/*Newton法の始まり値*/
    double alpha_fin;/*Newton法の更新値*/
    double x[400], y[400], z[400];/*walkerの位置配列*/
    double x_backup[400], y_backup[400], z_backup[400];
    double r_ini[400];
    double r_fin[400];
    double r[400];
    double psi_ini[400], psi_fin[400];
    double Sum_E_w[400]={0};
    double E_w[400]={0};
    double E=0;
    double Sum_Edlnpsi_w[400]={0};
    double Edlnpsi_w[400]={0};
    double Edlnpsi=0;
    double Sum_dlnpsi_w[400]={0};
    double dlnpsi_w[400]={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_iniでのエネルギーの微分dE/dalphaを計算する.*/
        /************************************************/

        shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*sumの配列を初期化する.*/

        shokihaichi(x, y, z);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, y, z, x_backup, y_backup, z_backup);/*提案分布に従ってwalkerを移動させる.*/

            /*Metropolis Testを行う*/
            
            for ( j = 0; j < 400; j++)
            {
                r_ini[j]=sqrt(x_backup[j]*x_backup[j]+y_backup[j]*y_backup[j]+z_backup[j]*z_backup[j]);
                psi_ini[j]=exp(-2*alpha_ini*r_ini[j]);/*元の確率分布*/

                r_fin[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                psi_fin[j]=exp(-2*alpha_ini*r_fin[j]);/*更新値での確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                    y[j]=y_backup[j];
                    z[j]=z_backup[j];
                }
                
                /*Metropolis Testの終了*/
                /*変数を更新し終えた*/

                /*walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.*/
                if (i>=4000)
                {   
                    r[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                    Sum_E_w[j]=Sum_E_w[j]-1/r[j]-0.5*alpha_ini*(alpha_ini-2/r[j]);
                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-r[j])*(-1/r[j]-0.5*alpha_ini*(alpha_ini-2/r[j]));
                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-r[j];

                    E_w[j]=Sum_E_w[j]/(i-3999);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);
                }
            }
        }
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 400; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }
        
        E=E/400;
        Edlnpsi=Edlnpsi/400;
        dlnpsi=dlnpsi/400;/*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(x, y, z);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, y, z, x_backup, y_backup, z_backup);/*提案分布に従ってwalkerを移動させる.*/

            /*Metropolis Testを行う*/
            
            for ( j = 0; j < 400; j++)
            {
                r_ini[j]=sqrt(x_backup[j]*x_backup[j]+y_backup[j]*y_backup[j]+z_backup[j]*z_backup[j]);
                psi_ini[j]=exp(-2*alpha_plus*r_ini[j]);/*元の確率分布*/

                r_fin[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                psi_fin[j]=exp(-2*alpha_plus*r_fin[j]);/*更新値での確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                    y[j]=y_backup[j];
                    z[j]=z_backup[j];
                }
                
                /*Metropolis Testの終了*/
                /*変数を更新し終えた*/

                /*walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.*/
                if (i>=4000)
                {   
                    r[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                    Sum_E_w[j]=Sum_E_w[j]-1/r[j]-0.5*alpha_plus*(alpha_plus-2/r[j]);
                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-r[j])*(-1/r[j]-0.5*alpha_plus*(alpha_plus-2/r[j]));
                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-r[j];

                    E_w[j]=Sum_E_w[j]/(i-3999);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);
                }
            }
        }
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 400; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }
        
        E=E/400;
        Edlnpsi=Edlnpsi/400;
        dlnpsi=dlnpsi/400;/*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(x, y, z);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, y, z, x_backup, y_backup, z_backup);/*提案分布に従ってwalkerを移動させる.*/

            /*Metropolis Testを行う*/
            
            for ( j = 0; j < 400; j++)
            {
                r_ini[j]=sqrt(x_backup[j]*x_backup[j]+y_backup[j]*y_backup[j]+z_backup[j]*z_backup[j]);
                psi_ini[j]=exp(-2*alpha_minus*r_ini[j]);/*元の確率分布*/

                r_fin[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                psi_fin[j]=exp(-2*alpha_minus*r_fin[j]);/*更新値での確率分布*/

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                    y[j]=y_backup[j];
                    z[j]=z_backup[j];
                }
                
                /*Metropolis Testの終了*/
                /*変数を更新し終えた*/

                /*walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.*/
                if (i>=4000)
                {   
                    r[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                    Sum_E_w[j]=Sum_E_w[j]-1/r[j]-0.5*alpha_minus*(alpha_minus-2/r[j]);
                    Sum_Edlnpsi_w[j]=Sum_Edlnpsi_w[j]+(-r[j])*(-1/r[j]-0.5*alpha_minus*(alpha_minus-2/r[j]));
                    Sum_dlnpsi_w[j]=Sum_dlnpsi_w[j]-r[j];

                    E_w[j]=Sum_E_w[j]/(i-3999);
                    Edlnpsi_w[j]=Sum_Edlnpsi_w[j]/(i-3999);
                    dlnpsi_w[j]=Sum_dlnpsi_w[j]/(i-3999);
                }
            }
        }
        E=0;/*変数の初期化*/
        Edlnpsi=0;
        dlnpsi=0;
        /*walker全ての平均を取ることにする*/
        for ( j = 0; j < 400; j++)
        {
            E=E+E_w[j];
            Edlnpsi=Edlnpsi+Edlnpsi_w[j];
            dlnpsi=dlnpsi+dlnpsi_w[j];
        }
        
        E=E/400;
        Edlnpsi=Edlnpsi/400;
        dlnpsi=dlnpsi/400;/*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\n", alpha_fin, alpha_ini);

    } while (fabs(alpha_fin-alpha_ini)>1e-6);

    double alpha_va=alpha_fin;/*変分で最適化した変分パラメータ*/

    printf("エネルギーが最小になる変分パラメータは%.10f\n", alpha_va);

    /******************************************************/
    /*alpha_vaでの変分で求めた基底状態のエネルギーを計算しよう*/
    /******************************************************/
    
    shokika(Sum_E_w, Sum_Edlnpsi_w, Sum_dlnpsi_w);/*sumの配列を初期化する.*/

    shokihaichi(x, y, z);

    for ( i = 0; i < sample; i++)
    {
        metropolis_walk(x, y, z, x_backup, y_backup, z_backup);/*提案分布に従ってwalkerを移動させる.*/

        /*Metropolis Testを行う*/
        
        for ( j = 0; j < 400; j++)
        {
            r_ini[j]=sqrt(x_backup[j]*x_backup[j]+y_backup[j]*y_backup[j]+z_backup[j]*z_backup[j]);
            psi_ini[j]=exp(-2*alpha_va*r_ini[j]);/*元の確率分布*/

            r_fin[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
            psi_fin[j]=exp(-2*alpha_va*r_fin[j]);/*更新値での確率分布*/

            double metropolis=genrand_real1();

            if (psi_fin[j]/psi_ini[j]>metropolis)
            {
                acc[j]=acc[j]+1;
            }
            else
            {
                x[j]=x_backup[j];
                y[j]=y_backup[j];
                z[j]=z_backup[j];
            }
            
            /*Metropolis Testの終了*/
            /*変数を更新し終えた*/

            /*walkerについて<E>や<Edlnpsi>,<dlnpsi>を計算する.*/
            if (i>=4000)
            {   
                r[j]=sqrt(x[j]*x[j]+y[j]*y[j]+z[j]*z[j]);
                Sum_E_w[j]=Sum_E_w[j]-1/r[j]-0.5*alpha_va*(alpha_va-2/r[j]);

                E_w[j]=Sum_E_w[j]/(i-3999);
            }
        }
    }
    E=0;/*変数の初期化*/
    /*walker全ての平均を取ることにする*/
    for ( j = 0; j < 400; j++)
    {
        E=E+E_w[j];
    }
    
    E=E/400;

    printf("基底状態のエネルギーは%.10f\n", E);


    return 0;
}

int shokihaichi(double data1[], double data2[], double data3[])
{
    int i;
    for ( i = 0; i < 400; 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[400];
    double q[400];
    double pi=asin(1e0)*2e0;
    for ( i = 0; i < 400; 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 < 400; i++)
    {
        data1[i]=0;
        data2[i]=0;
        data3[i]=0;
    }
    return 0;
}

このプログラムを実行すると、 $\alpha=1$ の時にエネルギーが最小になり、その時のエネルギーは $\langle E\rangle=-1/2$ になることがわかる。この結果は厳密解と一致している。

MiyanTarumi’s blog

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

変分モンテカルロ(VMC)をやってみる~水素原子~

水素原子

モデル

VMCの方法とソースコード