変分モンテカルロ(VMC)をやってみる~調和振動子~

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

モデル

まずは一番簡単な例として1次元調和振動子についてVMCをやってみることにしよう。
この系のHamiltonianは、

$H=\displaystyle -\frac{1}{2}\frac{d^2}{dx^2}+\frac{1}{2}x^2$

である。よく知られているようにこのHamiltonianの基底状態の解析解は $\exp(-x^2/2)$ であるから、ここでは試行関数として $\psi(x)=\exp(-\alpha x^2)$ を用いることにしよう。ただし $\alpha$ は変分パラメータである。
このように試行波動関数を決めたとしてHamiltonianに対する局所エネルギー $E_L$ を計算してみよう。これは単純に計算すればよく、

$H\psi(x)=\left[\displaystyle -\frac{1}{2}\frac{d^2}{dx^2}+\frac{1}{2}x^2\right]e^{-\alpha x^2}\\ =\left[\displaystyle-\frac{1}{2}\frac{d}{dx}\left(-2\alpha x\right)+\frac{1}{2}\frac{d}{dx}\left(-2\alpha x\right)\right]e^{-\alpha x^2}\\ =\left[\displaystyle \alpha+x^2\left(\frac{1}{2}-2\alpha^2\right)\right]e^{-\alpha x^2}$

であるから、局所エネルギー $E_L$ はその定義から

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

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

このプログラムを書いてみよう。

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

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

int shoki(double data1[]);/*プロトタイプ宣言*/
int metropolis_walk(double data1[], double data2[], double stepsize);/*プロトタイプ宣言*/
int shokika(double data1[], double data2[], double data3[]);/*プロトタイプ宣言*/

int main(void)
{
    int i, j;
    int sample=30000;

    double alpha_ini;/*Newton法のためのalphaの元値*/
    double alpha_fin;/*alphaの更新値*/
    double alpha_plus;
    double alpha_minus;
    double alpha_0=0.4;/*alphaの初期値*/
    double stepsize=0.7;
    double x[400];
    double x_backup[400];
    double psi_ini[400];
    double psi_fin[400];
    double ene_walkersum[400]={0};
    double ene_walker[400]={0};
    double E=0;
    double EDpsi_sum[400]={0};
    double EDpsi[400]={0};
    double Edlnpsidalpha=0;
    double Dpsi_sum[400]={0};
    double Dpsi[400]={0};
    double dlnpsidalph=0;
    double acc[400]={0};
    double dEdalpha;
    double h=1e-3;/*数値微分の幅*/

    init_genrand((unsigned)time(NULL));

    alpha_fin=alpha_0;
    alpha_ini=alpha_fin+1;/*while文に引っかかるようにずらしておく*/

    while (fabs(alpha_fin-alpha_ini)>1e-6)
    {
        /*alphaを変化させながら最小になるようなalphaを探しに行く*/
        alpha_ini=alpha_fin;/*前ループの結果を代入*/

        shokika(ene_walkersum, EDpsi_sum, Dpsi_sum);/*sumを収納する変数の初期化を行う.*/
        
        /*メトロポリス法を行ってenergyの微分を求める*/
        shoki(x);

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, x_backup, stepsize);/*walkさせる*/

            /*MetropolisTestを行う*/
            for ( j = 0; j < 400; j++)
            {
                psi_ini[j]=exp(-2*alpha_ini*x_backup[j]*x_backup[j]);
                psi_fin[j]=exp(-2*alpha_ini*x[j]*x[j]);

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc[j]=acc[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                }
                /*MetropolisTestの終了*/

                /*各walkerでの<E>, <E dln(psi)/d(alpha)>, <dln(psi)/d(alpha)>の計算 */
                if (i>=4000)
                {
                    ene_walkersum[j]=ene_walkersum[j]+(alpha_ini+x[j]*x[j]*(0.5-2*alpha_ini*alpha_ini));
                    ene_walker[j]=ene_walkersum[j]/(i-3999);

                    EDpsi_sum[j]=EDpsi_sum[j]+(-x[j]*x[j])*(alpha_ini+x[j]*x[j]*(0.5-2*alpha_ini*alpha_ini));
                    EDpsi[j]=EDpsi_sum[j]/(i-3999);

                    Dpsi_sum[j]=Dpsi_sum[j]-x[j]*x[j];
                    Dpsi[j]=Dpsi_sum[j]/(i-3999);
                }
                 
            }
        
        }
        /*walker全ての平均を取る*/
        /*初期化を行う*/
        E=0;
        Edlnpsidalpha=0;
        dlnpsidalph=0;
        for ( j = 0; j < 400; j++)
        {
            E=E+ene_walker[j];
            Edlnpsidalpha=Edlnpsidalpha+EDpsi[j];
            dlnpsidalph=dlnpsidalph+Dpsi[j];
        }
        E=E/400;/*400個の平均だから割っておく*/
        Edlnpsidalpha=Edlnpsidalpha/400;
        dlnpsidalph=dlnpsidalph/400;
        
        dEdalpha=2*(Edlnpsidalpha-E*dlnpsidalph);/*dE/dalphaが計算できた*/


        alpha_plus=alpha_ini+h;
        alpha_minus=alpha_ini-h;

        /***********************************************************************/
        /*dE/dalphaの数値微分を求めるためにalpha+とalpha-でのdE/dalphaの値を求める.*/
        /***********************************************************************/


        /****************************************/
        /****************************************/
        /*alpha+でのdE/dalphaを計算することにする.*/
        /****************************************/
        /****************************************/



        shoki(x);
        shokika(ene_walkersum, EDpsi_sum, Dpsi_sum);/*再利用するためにsumを収納する変数の初期化を行う.*/

        double acc_plus[400]={0};
        double dEdalpha_plus;

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, x_backup, stepsize);/*walkさせる*/

            /*MetropolisTestを行う*/
            for ( j = 0; j < 400; j++)
            {
                psi_ini[j]=exp(-2*alpha_plus*x_backup[j]*x_backup[j]);
                psi_fin[j]=exp(-2*alpha_plus*x[j]*x[j]);

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc_plus[j]=acc_plus[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                }
                /*MetropolisTestの終了*/

                /*各walkerでの<E>, <E dln(psi)/d(alpha)>, <dln(psi)/d(alpha)>の計算 */
                if (i>=4000)
                {
                    ene_walkersum[j]=ene_walkersum[j]+(alpha_plus+x[j]*x[j]*(0.5-2*alpha_plus*alpha_plus));
                    ene_walker[j]=ene_walkersum[j]/(i-3999);

                    EDpsi_sum[j]=EDpsi_sum[j]+(-x[j]*x[j])*(alpha_plus+x[j]*x[j]*(0.5-2*alpha_plus*alpha_plus));
                    EDpsi[j]=EDpsi_sum[j]/(i-3999);

                    Dpsi_sum[j]=Dpsi_sum[j]-x[j]*x[j];
                    Dpsi[j]=Dpsi_sum[j]/(i-3999);
                }
                 
            }
        
        }
        /*walker全ての平均を取る*/
        /*初期化を行う*/
        E=0;
        Edlnpsidalpha=0;
        dlnpsidalph=0;
        for ( j = 0; j < 400; j++)
        {
            E=E+ene_walker[j];
            Edlnpsidalpha=Edlnpsidalpha+EDpsi[j];
            dlnpsidalph=dlnpsidalph+Dpsi[j];
        }
        E=E/400;/*400個の平均だから割っておく*/
        Edlnpsidalpha=Edlnpsidalpha/400;
        dlnpsidalph=dlnpsidalph/400;
        
        dEdalpha_plus=2*(Edlnpsidalpha-E*dlnpsidalph);/*dE/dalpha+が計算できた*/



        /****************************************/
        /****************************************/
        /*alpha-でのdE/dalphaを計算することにする.*/
        /****************************************/
        /****************************************/


        shoki(x);
        shokika(ene_walkersum, EDpsi_sum, Dpsi_sum);/*再利用するためにsumを収納する変数の初期化を行う.*/

        double acc_minus[400]={0};
        double dEdalpha_minus;

        for ( i = 0; i < sample; i++)
        {
            metropolis_walk(x, x_backup, stepsize);/*walkさせる*/

            /*MetropolisTestを行う*/
            for ( j = 0; j < 400; j++)
            {
                psi_ini[j]=exp(-2*alpha_minus*x_backup[j]*x_backup[j]);
                psi_fin[j]=exp(-2*alpha_minus*x[j]*x[j]);

                double metropolis=genrand_real1();

                if (psi_fin[j]/psi_ini[j]>metropolis)
                {
                    acc_plus[j]=acc_plus[j]+1;
                }
                else
                {
                    x[j]=x_backup[j];
                }
                /*MetropolisTestの終了*/

                /*各walkerでの<E>, <E dln(psi)/d(alpha)>, <dln(psi)/d(alpha)>の計算 */
                if (i>=4000)
                {
                    ene_walkersum[j]=ene_walkersum[j]+(alpha_minus+x[j]*x[j]*(0.5-2*alpha_minus*alpha_minus));
                    ene_walker[j]=ene_walkersum[j]/(i-3999);

                    EDpsi_sum[j]=EDpsi_sum[j]+(-x[j]*x[j])*(alpha_minus+x[j]*x[j]*(0.5-2*alpha_minus*alpha_minus));
                    EDpsi[j]=EDpsi_sum[j]/(i-3999);

                    Dpsi_sum[j]=Dpsi_sum[j]-x[j]*x[j];
                    Dpsi[j]=Dpsi_sum[j]/(i-3999);
                }
                 
            }
        
        }
        /*walker全ての平均を取る*/
        /*初期化を行う*/
        E=0;
        Edlnpsidalpha=0;
        dlnpsidalph=0;
        for ( j = 0; j < 400; j++)
        {
            E=E+ene_walker[j];
            Edlnpsidalpha=Edlnpsidalpha+EDpsi[j];
            dlnpsidalph=dlnpsidalph+Dpsi[j];
        }
        E=E/400;/*400個の平均だから割っておく*/
        Edlnpsidalpha=Edlnpsidalpha/400;
        dlnpsidalph=dlnpsidalph/400;
        
        dEdalpha_minus=2*(Edlnpsidalpha-E*dlnpsidalph);/*dE/dalpha-が計算できた*/

        double ddE=(dEdalpha_plus-dEdalpha_minus)/(2*h);/*dE/dalphaのalpha_iniでの数値微分*/

        alpha_fin=alpha_ini-dEdalpha/ddE;/*alphaの更新*/

        printf("%.10f %.10f %.10f\n", alpha_fin, alpha_ini, dEdalpha);

    }

    double alpha_va=alpha_fin;/*変分で最小化された変分パラメータ*/
    double E_val;/*変分で求めたエネルギー*/

    /*alpha_valでのエネルギーを計算*/

    shoki(x);
    shokika(ene_walkersum, EDpsi_sum, Dpsi_sum);/*再利用するためにsumを収納する変数の初期化を行う.*/

    double acc_va[400]={0};


    for ( i = 0; i < sample; i++)
    {
        metropolis_walk(x, x_backup, stepsize);/*walkさせる*/

        /*MetropolisTestを行う*/
        for ( j = 0; j < 400; j++)
        {
            psi_ini[j]=exp(-2*alpha_va*x_backup[j]*x_backup[j]);
            psi_fin[j]=exp(-2*alpha_va*x[j]*x[j]);

            double metropolis=genrand_real1();

            if (psi_fin[j]/psi_ini[j]>metropolis)
            {
                acc_va[j]=acc_va[j]+1;
            }
            else
            {
                x[j]=x_backup[j];
            }
            /*MetropolisTestの終了*/

            /*各walkerでの<E>, <E dln(psi)/d(alpha)>, <dln(psi)/d(alpha)>の計算 */
            if (i>=4000)
            {
                ene_walkersum[j]=ene_walkersum[j]+(alpha_va+x[j]*x[j]*(0.5-2*alpha_va*alpha_va));
                ene_walker[j]=ene_walkersum[j]/(i-3999);
            }
             
        }
    
    }
    /*walker全ての平均を取る*/
    /*初期化を行う*/
    E=0;
    for ( j = 0; j < 400; j++)
    {
        E=E+ene_walker[j];
    }
    E=E/400;/*400個の平均だから割っておく*/
    E_val=E;

    printf("変分で求めたエネルギーは%.10f\n", E_val);
    return 0;
}

int shoki(double data1[])
{
    int i;
    for ( i = 0; i < 400; i++)
    {
        data1[i]=(genrand_real1()-0.5)*2*2;/*[-2,2]の間にランダムに初期*/
    }
    return 0;
}

int metropolis_walk(double data1[], double data2[], double stepsize)/*data1がx, data2がx_backup*/
{
    int i;
    double p[400];
    double q[400];
    double pi=asin(1e0)*2e0;
    for ( i = 0; i < 400; i++)
    {
        data2[i]=data1[i];/*data1の情報をdata2に保存*/
        p[i]=genrand_real1();
        q[i]=genrand_real1();
        data1[i]=data1[i]+stepsize*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=0.5$ が得られ、その時にエネルギーは、 $\langle E\rangle=1/2$ になる。これは厳密解と一致する。

MiyanTarumi’s blog

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

変分モンテカルロ(VMC)をやってみる~調和振動子~

調和振動子

モデル

VMCの方法とソースコード