変分モンテカルロ(VMC)をやってみる~ヘリウム原子~
Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。
ヘリウム原子
モデル
最後にヘリウム原子についても同じくVMCをやってみよう。
この系のHamiltonianはBorn-Oppenheimer近似の下でを電子1の座標、を電子2の座標とすると、
と書ける。これまでの時と同様に試行関数を与えたいのだが、この場合は良い試行関数を作ること自体が問題になってしまうので、ここでは天下り的であるが試行関数として次のようなものを使用することにしよう。
ここでとし、これまでと同様には変分パラメータである。
次にこの試行関数を用いた時に局所エネルギーがどのようになるのかを計算してみる*1。
微分演算子の入っていない部分はそのまま係数としてかかるだけなので微分演算子の部分だけを計算すればいいだろう。そこでまずについて計算することにしよう。合成関数のラプラシアンに対する作用が、
であることに注意すると、に関係する部分の計算は、
を計算すればよいことになる。まず初項について計算しよう。この項の計算は簡単で1体の時のラプラシアンの計算と全く同じであるから、
より、
となる。
続いて第二項について計算する。ここでこれ以降の計算の表記を簡単にするためにが入っている部分の指数の肩を
と書くことにし、その微分を
と書くことにする。すると第二項の計算は、
となる。ここでは電子1の座標の成分を、は電子2の座標の成分をそれぞれ表すとした。
最後に第三項の計算を行おう。この項ではの入った項の微分を二回行う必要があるのでの二回微分を表すものとして
とを定義しよう。すると第三項は、
となる。以上から
となることがわかった。同様にしてについて計算すると、
を得ることができる*2。
以上から
と計算でき、よって局所エネルギーは、
と求めることができた。
VMCの方法とソースコード
よってプログラミングで行う操作は次の通りである。
**REPAT**
始めの変分パラメータを決める。
ランダムな場所に個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**
1. 番目のサンプリング点をとする。walkerの移動後の座標として番目の候補をを中心とする正規分布に従ってとして提案する。
2. メトロポリステスト:
の間の一様乱数を生成し、であれば提案を受理してとする。そうでなければ候補を棄却してとする。
3. サンプリング点での局所エネルギーを計算する。
**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値を求める。
Newton法に従って変分パラメータを更新する。
**END**
ここでは座標であり、今の場合である。
このプログラムを書いてみよう。
▶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; }
このプログラムを実行すると、となり、その時の基底状態のエネルギーとしてが得られる。これはHartree-Fockの値やDFTの値、厳密値と比肩しうる数値になっている。
変分モンテカルロ(VMC)をやってみる~水素原子~
Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。
水素原子
モデル
調和振動子に続いて水素原子についても同じくVMCをやってみよう。水素原子の問題は1次元の問題に帰着できるが、ここでは3次元の問題として扱うことにする。
この系のHamiltonianはBorn-Oppenheimer近似の下で、
と書ける。水素原子の基底状態の厳密解はであるから、ここでは試行関数としてを選ぶことにしよう。ここでは変分パラメータである。このように試行関数を決めた時に局所エネルギーはどうなるかというと、
より、
となる。
VMCの方法とソースコード
よってプログラミングで行う操作は次の通りである。
**REPAT**
始めの変分パラメータを決める。
ランダムな場所に個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**
1. 番目のサンプリング点をとする。walkerの移動後の座標として番目の候補をを中心とする正規分布に従ってとして提案する。
2. メトロポリステスト:
の間の一様乱数を生成し、であれば提案を受理してとする。そうでなければ候補を棄却してとする。
3. サンプリング点での局所エネルギーを計算する。
**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値を求める。
Newton法に従って変分パラメータを更新する。
**END**
ここでは座標であり、今の場合である。
このプログラムを書いてみよう。
▶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; }
このプログラムを実行すると、の時にエネルギーが最小になり、その時のエネルギーはになることがわかる。この結果は厳密解と一致している。
変分モンテカルロ(VMC)をやってみる~調和振動子~
Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。
調和振動子
モデル
まずは一番簡単な例として1次元調和振動子についてVMCをやってみることにしよう。
この系のHamiltonianは、
である。よく知られているようにこのHamiltonianの基底状態の解析解はであるから、ここでは試行関数としてを用いることにしよう。ただしは変分パラメータである。
このように試行波動関数を決めたとしてHamiltonianに対する局所エネルギーを計算してみよう。これは単純に計算すればよく、
であるから、局所エネルギーはその定義から
と求まる。
VMCの方法とソースコード
よってプログラミングで行う操作は次の通りである。
**REPAT**
始めの変分パラメータを決める。
ランダムな場所に個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**
1. 番目のサンプリング点をとする。walkerの移動後の座標として番目の候補をを中心とする正規分布に従ってとして提案する。
2. メトロポリステスト:
の間の一様乱数を生成し、であれば提案を受理してとする。そうでなければ候補を棄却してとする。
3. サンプリング点での局所エネルギーを計算する。
**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値を求める。
Newton法に従って変分パラメータを更新する。
**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; }
このプログラムを実行すると、が得られ、その時にエネルギーは、になる。これは厳密解と一致する。
変分モンテカルロ(VMC)をやってみる~準備~
Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。
変分モンテカルロ法(VMC)
量子力学の多くの問題においてSchrödinger方程式を厳密に解くことは難しく、何らかの形で近似的に解く必要がある。そこで変分による方法で基底状態について近似的に解くことを考えてみよう。
VMCのための前準備
量子力学での変分原理
Hamiltonianをとし、そのエネルギー固有値をとしたときにそれに対応した固有関数を(は座標)と書くことにしよう*1。
ここで唐突だが、ある波動関数(はパラメータ)を考えてみよう。この波動関数をエネルギー固有関数でモード展開すると次のように展開できる:
ここではで展開するときの展開係数であり、それは
で与えられる。
量子力学においてエネルギーの期待値は、
で計算できるのであった。ここにモード展開された表式での波動関数を代入してあげると、
となる。ここでが規格化されているとして規格化条件を用いた。
ここから波動関数の期待値は基底状態のエネルギーを下回ることはないことがわかる。なのでのパラメータを期待値が最小になるように調整すれば波動関数は基底状態の波動関数に近づくことになり、基底状態のエネルギーの近似値を得ることができる。
さて、ここまでは試行波動関数が規格化されているものとして話を進めてきたが、一般に試行波動関数として複雑なものを採用した場合は規格化を行うことは面倒であるし、そもそも積分を解析的に行うことが可能かどうかすら怪しいであろう。そこでこの先では試行波動関数があらかじめには規格化されていないとし、試行波動関数をと書くことにしよう。すると、エネルギー期待値は、
と書き直される。よって変分原理は、
のように書き直せばよい。
原理的にはこの変分原理に従ってエネルギーが最小になるように変分パラメータを調整すればよいのだが、先に述べたように試行関数やHamiltonianが複雑になれば解析的に実行することは困難であり、積分を区分求積的に数値計算で行うにしても1次元系ならともかく3次元などになれば問題が複雑になればなるほど積分の次元が飛躍的に大きくなり次元の呪いに嵌ってしまう。例えば3次元の2粒子系(Helium原子の電子など)では積分に必要な次元は6次元になる。
そこで高次元でも数値計算を行うためにマルコフ連鎖モンテカルロ法を用いて変分における積分計算を行うことにする(だから変分"モンテカルロ"という。)。
マルコフ連鎖モンテカルロ法(MCMC)
前節での要請に従ってMCMCを導入する。
モンテカルロ法
先にも述べたように高次元の数値積分は難しい。このことを実感するために領域として長さ1の次元超立方体を積分することを考えてみよう。仮にの間を1000分割して区分求積的に積分を行うことを考えると、積分領域は個の領域に分けられることになり、考えなければならない領域の数が指数関数的に増加してしまい高次元では計算不可能になってしまうことがわかるだろう。
この計算コストの指数的増加による次元の呪いから逃れるために積分領域からランダムに値を選んできてそれらの平均をとることで積分を求めるということを考えてみる。同様にの領域での積分を考えよう。ランダムに点を選んで持ってきて平均を取るというアイデアに従えば積分は、
と表される。この式自体はをランダムに選ぶことを除けば区分求積とほぼ同じ形なのでが大きくなれば厳密な積分に近づくことは容易にわかる。そこでこの近づき方がどのように振る舞うかを確認してみよう。そのためにはこの式の分散の振る舞いを見ればよい。の分散は次の式で表せる。
この式において初項のの時は後ろの項と相殺して消え、のときに初項からの期待値が個でることを考えて計算すると、
となり、のように次元によらない形で振る舞うことがわかる。
以上からモンテカルロ法を用いて積分計算をすれば次元の呪いを回避しながら数値計算できそうだと予想できるだろう。
マルコフ連鎖モンテカルロ法
さて、上の例ではランダムに点を選ぶときに点が選ばれる確率は一様にサンプリングされるとした。しかし、高次元の積分を行うときには一様にサンプリングすることは計算資源の無駄遣いになってしまう。というのも多くの被積分関数は領域によって積分への寄与が全く異なり、次元が上がれば積分への寄与が大きい適切なサンプリング点がサンプリングされる確率が下がってしまうからである。例えばGauss積分を考え、サンプリングはから行うとする。よく知られているようにGauss関数はに近いところでのみ大きな値を持ち原点から離れるにしたがって急速にに近づく。そこで有意な値を持つ範囲をであるとすると、のときの有意なサンプリング領域はである。しかし、では有意な領域は全積分領域のうちしかない。こうしたことから高次元で安易に一様サンプリングをしてしまうと大きな統計誤差を含んだ結果が出てしまう。つまり一見モンテカルロ法を用いて次元の呪いを回避できたように見えてサンプリングの問題として残ってしまっている。
このようなサンプリングの問題を踏まえると、サンプリング数を高次元でも膨大なものにしないために一様にサンプリングを行うのではなく、積分領域のうち有意な領域を重点的にサンプリングすればいいだろう。この有意な領域への重点的なサンプリングを行うために望んでいる確率分布を構成する方法の一つとしてマルコフ連鎖モンテカルロ法(MCMC)がある。
構成したい確率分布をとしよう。MCMCではの値を各点が現れる滞在時間をに比例するように変化させるようにし、これによって十分長い時間で統計平均の良い近似を得ることができる。そこでこのようなMCMCの方法として今回用いるメトロポリス法によるMCMCを導入することにしよう*2。
メトロポリス法
MCMCの一つの方法としてメトロポリス法の手順を確認する。
先と同様に構成したい確率分布をとする。メトロポリス法の手続きは次のようなものである。
1. 番目のサンプリング点をとする。実数をランダムに選び、番目の候補をとして提案する。ただし、MCMCの条件を満たすために:とが同じ確率で現れるように候補の選び方を提案する。
2. メトロポリステスト:
の間の一様乱数を生成し、であれば提案を受理してとする。そうでなければ候補を棄却してとする(この条件はで提案を受理するということと同じである。)。
1において述べたようにとが出現する確率が同じであればよいのでその提案確率はの間の一様乱数でもを中心としたガウス乱数でも何でもよい。なお、記事でVMCをする際にはガウス乱数を用いて計算することにする。
VMC
このようにして得られるメトロポリス法によるMCMCを用いてどのように変分計算を行うかを考えてみよう。
量子力学のおける変分原理によれば、試行関数を最適化するためには、
を最小化すればよいのであった。唐突だがこの式を少し変形してみよう。形を確率分布を挟んだ積分らしく変形すると、は局所エネルギーを用いて
と書くことができる。ここでとについて
とした。
このように書き換えることでエネルギーの期待値は確率分布がである局所エネルギーの期待値と思うことができる。そこでMCMCを用いてサンプリング点をから選び出して平均を取ってやれば大数の法則からエネルギーの期待値を求めることができるであろう。これがVMCにおける多次元積分の基本的な方針である。
最後にここまでの話を踏まえてVMCでのメトロポリス法を整理しておこう。
**REPAT**
始めの変分パラメータを決める。
ランダムな場所に個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**
1. 番目のサンプリング点をとする。walkerの移動後の座標として番目の候補を提案分布に従ってとして提案する。ただし、MCMCの条件を満たすために:とが同じ確率で現れるように候補の選び方を提案する。
2. メトロポリステスト:
の間の一様乱数を生成し、であれば提案を受理してとする。そうでなければ候補を棄却してとする。
3. サンプリング点での局所エネルギーを計算する。
**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値などを求める。
Newton法に従って変分パラメータを更新する。
**END**
Newton法
最後に変分パラメータの調整の仕方を確認する。
エネルギーが最小になるではのによる微分がになるのでのゼロ点を求めればよい。自体を数値微分したものを使ってもよいが二階導関数を数値微分の数値微分で求めることになるので大きな誤差が出てしまう。そこで一階導関数は解析的に計算することにしよう。の表式に戻って微分を行うと、
となる。ここで試行関数として実関数のみを用いることを仮定すると、計算結果として現れるのは実数であることに注意すると、
と右側に落とした微分を左側に持ってくることができる。よっては、
と計算できる。これをMCMCを用いて数値計算し、あとはNewton法の更新式
に従って変分パラメータを更新すればよい。ただしはの数値微分である。
さて、これで準備は終わりなので以上を踏まえて具体的な系で計算をしてみよう。
参考文献
[1]. J. M. Thijssen, 「計算物理学」,丸善出版, 2012.
[2]. 花田 政範, 松浦 壮, 「ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門」, 講談社, 2020.
2022年度 大阪大学 大学院理学研究科 物理 解答
阪大の理学研究科 物理学・宇宙地球科学専攻の2022年度院試問題の解答です.間違いがあれば教えてください.
問題は,ここにあります.(9/10現在はまだ上がってませんが,そのうちアップロードされると思います.)
問題I
I
(1)
密度分布はデルタ関数と次で定義される関数
を用いると,
と書ける.よって棒の慣性モーメントは,回転軸を軸とすると,
と求められる.
(2)
棒に働く力のモーメントがなことから回転の運動方程式は,が十分小さいのでと近似できることに注意すると,
となる.よって棒の振り子運動の角振動数は,
とわかる.
(3)
(2)で得た運動方程式の一般解は,
と表せる.ここで,は積分定数である.
まずでであることからであることが従う.さらに角運動量と力積の関係から
が成立するので,
よりは,
と求められる.以上からは,
と決まる.
よって振れ幅の最大値は,
と求められる.
II
(4)
重心の座標についてこの時間微分は
であるから,重心の運動エネルギーは
となる.よって台Aと棒からなる系の運動エネルギーはこれと棒の重心周りの回転エネルギー,そして台の運動エネルギーを足して
となる.
一方で位置エネルギーについて基準面をにとることと,軸を鉛直上向きに取ることに注意すると,
となる.
(5)
について十分小さいとし,と近似することにする.,に関する3次以上の項及び定数項を無視することにすれば,Lagrangianは,
となる.よって,,,,とわかる.
(6)
についてのLagrange方程式を求めると
となる.一方でについてのLagrange方程式は,
となる.
(7)
の時を考える.この時についてのLagrange方程式において効くのはがかかっている項のみであるから,この極限で運動方程式は
となる.において台がにあり,そして速度がであることを踏まえると台は動かない.
をに関する運動方程式に代入すれば
となる.よって振り子運動の角振動数は,
とわかる.
(8)
の時を考えよう.このときについての運動方程式は
となる.棒の重心に注目する.の二回微分について計算してみると,
となる.重心の初速度がであることを考えると,台Aと棒は固定された重心のの周りで運動する.
最後に棒の振り子運動の角振動数を求める.についての運動方程式からであるから,これをに関する運動方程式に代入すると,
となる.よって角振動数は,
である.
(9)
まず,についてのの運動方程式から
となる.これをについての運動方程式に代入すると,
と変形できる.よっては角振動数
の単振動をする.ゆえに,この方程式の一般解は積分定数,を用いて
と書ける.初期条件についてにおいて,であることを踏まえれば,,と決まる.以上からまずのついて解が
となることがわかる.
一方についての運動方程式について両辺を積分すると,,を積分定数として
となる.において,であることを考えると,,であることがわかる.よっては,
と求まる.
問題2
I
(1)
円筒座標系としてを取る.この時に電荷分布は,
と表される.基底ベクトルを用いてとを,と表すことにすれば,
となる.よって静電ポテンシャルは,
と求められる.
(2)
であるから,円筒座標におけるナブラ演算子が
であることと,がにのみ依存することに注意すると,
と電場が求められる.原点において電場は明らかにであるから電荷が受ける力はである.
(3)
であるからである.
なので,
を計算すればよい.ここで,,の積分はになることを踏まえてnon zeroな項のみを抜き出せば,
と求められる.念のために(2)の答えとの整合性を確認しておくとの下で近似を行い最小次の項のみ残せば(3)での解と一致する.
II
(4)
同じく円筒座標を用いて考える.,とすると,
と近似できるので,ベクトルポテンシャルは,であることから,
となる.再びここでもnon zeroな項のみ抜き出してやることにすれば,
とベクトルポテンシャルが計算できる.
(5)
であるから,これを計算すればよい.磁束密度の,成分についてはがに依存しないことやがであることが容易にわかるのでのみ求めればよいであろう.実際に計算してみると,
となる.以上から磁束密度は,
と計算できた.
(6)
円周上の電荷は秒間に回,回転することになる.ここでは角速度である.また角運動量 の大きさについては,慣性モーメントを用いてと書けるのであった.よって電流は,
と表せる.よって磁束密度についてこれを代入すれば
ととの関係がわかる最後に磁気モーメントと磁束密度の相互作用エネルギーについて,その定義より,
となる.
III
(7)
であるから,磁束密度は,
と求められる.さて,
であるから
とも書ける.よって最終的に
となる.
(8)
わからない.
問題3
I
(1)
であるから
(2)
となる.ここでデルタ関数がかかった項の積分についてでもでもどちらでも結果は変わらないが簡単な方を選んだ.ここでのlimitをとると右辺の積分はになることから,上の式は
と変形される.を用いると,接続条件は
となる.
(3)
(1)と(2)の結果をとについて解くと
となる.よって,
となる.
(4)
であることと,なことからのときであり,であるときに対応する.
まず,のときであったから
となる.これは入射粒子のエネルギーがであるときには全く透過しないであろうと予想されることから考えれば当然の結果である.
一方でのときはであったから,
となる.この結果もエネルギーがであるときにはすべて壁を乗り越えうることを思えば妥当な結果になっている.
最後にであるときにとするとどうなるかを考える.この操作によっての変化は(2)の式を踏まえるととの虚部の符号が反転することのみである.つまり,これは複素共役をとったことに対応するが,複素共役をとったことによって絶対値は変化しないのでとはこの操作では変化しない.
II
(5)
とにおけるそれぞれのポテンシャル障壁でのやはこの場合でも変化しないと考えられる.
まずについてはが反射したときの寄与とが透過したときの寄与が考えられ,それらを足し合わせたものがになる.ゆえに,
となる.
続いてについては,が透過したときの寄与とが反射したときの寄与が考えることができる.よって,
が成立する.
(6)
まずについて
とまとめる.これを先の二式目に代入すると,
となる.以上から行列は,
と求められる.
(7)
に対する2つの表式について等置すれば
となり,
と求まる.よって,
と書くことができる.以上から
とわかった.
(8)
(7)における表式のとの関係についてこれらはどちらもが中心の波であるから(5)と同様に考えよう.
まずへの寄与としては,からの反射とからの透過の寄与が考えられるので
となる.一方でについては,からの透過とからの反射が考えられるので
となる.これをとについて解けばよいのだが,これは(5)の式と全く同一なので結果をそのまま流用すればよい.なので,
と書ける.これと(6),(7)での結果を踏まえれば
となる.以上からは,
とわかる.
(9)
であるからとなるのでは単位行列に他ならない.よってを求めるためにはを求めればよい.そこでまずをを用いて表すことにしよう.(3)の結果からそれぞれ
となる.よっては,
と表せる.よって,は,
と計算できる.以上から(8)において求めた式で,,,を代入すると
となる.ゆえにまず,
と求まる.そしてについても
となる.よって,,は,
と計算できる.最後にこれを(3)における結果と比較すると,
となる.つまり,II(9)における反射率の方がI(3)の反射率よりも大きいことがわかる.これはポテンシャル障壁が1つから2つになった結果により入射粒子がより反射されやすくなったからであると考えられる.
問題4
I
(1)
Gibbs free energyの全微分を計算し,Helmholtz free energyの全微分を代入すると,
となる.よって,,,である.
(2)
(1)で求めたGibbs free energyの全微分から
である.は今によらないのでで積分すると,
と,をとの2つの熱力学変数の積で表せる.そしてこの式の全微分を取り,それと(1)において求めたGibbs free energyの全微分を等置すると,
となり,Gibbs-Duhemの関係式を得る.ただしここで,,とした.
II
(3)
分配関数をまず計算しよう:
さて,Helmholtz free energyについてであるから,
とHelmholtz free energyが求まった.
(4)
Helmholtz free energyの全微分から化学ポテンシャルについて
となる.ここで今考えているものは理想気体であるから状態方程式を用いて表式からを消去すると,
と求まる.
(5)
一粒子状態密度(以下DOS)を求めるために,まずエネルギーが以下であるような状態の数を求めよう.
まず自由粒子のSchrödinger方程式
の解は
となる.ここでと置いた.今,周期境界条件などを課すことにすれば(は箱の一辺の長さ),
が得られる.
エネルギーが以下の状態の総数を求めるためには半径がである3次元球の中にある状態の数を数えればよいのだが,-空間において1つの状態が占める体積は周期境界条件より得られた条件式からであるので,球の体積をこれで割ればよい.つまり,は,
と求められる.よってDOSは,
となる.以上より,
とわかる.
(6)
Maxwellの関係式
より両辺をで積分すれば圧力が求められる.
さて,粒子数の平均について
で与えられたから,これをで微分し上式の右辺に代入すると次のようになる.
この両辺をで積分しよう.積分範囲はにとるが,でであるから積分の下限については省略することにする.
ここでは以上の定数であり,である.
上の積分においてと置換して積分公式を用いると,
となり,圧力の表式を得ることができる.
(7)
Helmholtz free energyの全微分から圧力は,
で求められるのであった.まず,風船の中の圧力について求めよう.風船内のHelmholtz free energyは,
であったからは,
と計算できる.一方,風船外部の圧力については,風船内部の気体と風船外部の気体は同一のものであるからそのHelmholtz free energyは風船内の気体のHelmholtz free energyと同一であると考えられる.ゆえには,
となる.よって圧力差は,
である.よって,とは,
とわかる.
(8)
であるから,が成立する.ここにGibbs-Duhemの関係式から得た式
を代入すると,
という式を得ることができる.
さらにの左辺に(7)の結果を代入し,そのうえで両辺の全微分を取ると
という式を得ることができるので,これを先に得た式のに代入すると,
となり,が満たす微分方程式が得られた.
最後に微分方程式をの一次のオーダーで近似する.のときと近似できることを用いると,
を得られる.
(9)
を(8)で得られた近似された微分方程式に代入すると,
となる.両辺を積分すると,
を得る.ここで(7)より,であるからこれを代入すると,
と解を求めることができた.最後に積分定数を境界条件から定めよう.与えられた境界条件はであるから,
となり,とわかる.以上よりは,
と求められた.
二重ポテンシャル障壁による共鳴透過
一次元量子系のポテンシャル障壁による散乱において、のとき、を障壁の長さとすると
を満たす場合、透過率がとなり障壁を反射されることなく完全に通り抜けることが知られていて、これは共鳴透過と呼ばれている。
ポテンシャル障壁が1つしかないときは、の場合にしか共鳴透過は起こらないのだが、ポテンシャル障壁が2つある場合にはの場合にもこの共鳴透過が起こる*1。このことを実際に散乱問題を解くことで確認してみよう。
問題設定
上のようなポテンシャルの散乱問題を考えよう。そして、における波動関数を、の波動関数を、の波動関数を、の波動関数を、の波動関数をと置くことにする。
Schrödinger方程式とその接続条件
各区間におけるSchrödinger方程式を並べてやると次のようになる。
これらを解けばを積分定数として、
ここで、
である。
接続条件は、ポテンシャルが有限なことを踏まえると各波動関数とその微分が各区間端において連続であればよい。ゆえに、
であり、これらの方程式から積分定数を決めればよい。しかし、これを強引に解くのはかなり骨が折れるので行列を用いて解くことにする。そのために、それぞれを整理して次のように書いてみよう。面倒なのでと書くことにして、
と書き直せる。積分定数を求めるにはこれらの連立方程式を解けばよい。
積分定数の決定1
さて、上で求めた連立方程式を解くためにまず準備をいくつか行っておく。
行列について
ととを定義しよう。
また、連立方程式を解くためにの逆行列を求めると、その逆行列は、
となっている。
以上の下で上の連立方程式を改めて書き直すと、
となる。ここで、共通して出てくるのはという行列なのでこれをと書くことにしよう。そしてこれを具体的に計算してみると、
となる。これを用いて連立方程式を解いてみる。上の連立方程式を順に代入していくと、
となり、行列の計算を行えばこの連立方程式を解けたことになる。
具体的に行列を計算していけばよいのだが、少々面倒なのでかけていく行列の規則性に注目してみよう。上の式を見ると、はとのペアの繰り返しで出てきている。そこでこの二つをまとめてと表すことにし、これを転送行列と呼ぶことにする。なお、連立方程式を見直せばこの転送行列一つ一つ一つがポテンシャル障壁を一つ通り抜けることに対応していることがわかるだろう。
転送行列の成分を具体的に計算してみよう。定義から
とわかる。
今の場合の転送行列2つにおいては共通しているので、の成分のうちのみによる部分を文字で置くことにしよう。
このように置いたうえで連立方程式の解をこれらを用いて表すと次のようになる。
ここで最後の表式についてなどの引数について省略した。
以上から、連立方程式を完全に解くことができては、
とわかる。
ここで物理的な要請から変数を減らす。粒子はから入射してくるとしよう。するとの領域では反射波がないのでと置ける。また、簡単のためとしよう。これらから残りの変数について完全に決定することができる。
まず、について2つ目の式から
となる。そして、についても
と計算できる。
積分定数の決定2
あとはを代入すればよいのだが、先にについてある程度式変形を行っておこう:
これらを踏まえたうえで計算を行おう。
まず、の分子の部分について計算を行う。
ここで
である。
次にの分子の部分の計算を行う。まず、
である。なお最後の式変形においてを用いた。
よって、
とわかった。
最後にとに共通する分母の部分の計算を行おう。
ここで、がその定義から純虚数であったことを思い出そう。このことに注意したうえで実部と虚部に上の式を分けると、実部について
と整理できる。
続いて虚部については
となる。よって、分母の部分は、
積分定数の決定3
以上からとについて
と求められた。
反射率と透過率
今入射波の絶対値をと置いているので反射率と透過率はそれぞれ係数の絶対値二乗で与えられることに注意しよう。
面倒なのは分母の部分なので分母のみを先に計算しておく。
となる。
よって反射率と透過率は、が純虚数なことに注意して
であることがわかった。これは明らかにを満たす。
こうして得られた透過率について、であるときに、つまりエネルギーがポテンシャルよりも小さいのにも関わらず粒子が完全に透過してしまうことがわかるだろう。これが二重ポテンシャル障壁における共鳴透過と呼ばれる現象であり、江崎玲於奈によるトンネル電流の最も簡単なモデルである。
共鳴透過が起こる条件を調べよう。の定義より、は次のように書ける。
そして、左辺について、、を代入して計算すると、
となるので、結局、共鳴透過の条件は、
となる。
そこで、の単位系において、横軸に入射エネルギーを取り、、、、としてプロットしたのが以下の図である。
上の図において、赤線は共鳴透過の条件の左辺、青線は条件式の右辺であり、これらの交点が共鳴透過が起こるようなエネルギーに対応している。そして、緑線は透過率の対数を取ったものであり、である点で共鳴透過が起こっている。
共鳴透過の条件と井戸型ポテンシャルの固有状態
なぜこのような不思議な現象が起きているかを考えてみよう。
そのヒントは、二重ポテンシャル障壁において障壁と障壁の間が有限井戸型ポテンシャルのようになっていることにある。
有限井戸型ポテンシャルの固有値を決める方程式は、
であった。なお、ここで用いた井戸型ポテンシャルのモデルは一般に使われるような軸対象の問題設定とは少し違うことに注意せよ。具体的には、井戸の始まりはであり、そこから]まで井戸の底が続き、井戸の終わりがにあるような問題設定である。
さて、これと先の共鳴透過の条件式を見比べると、が右辺にかかっているかどうかの違いしかないことがわかる。がでになることを考えれば、有限井戸型ポテンシャルは二重ポテンシャル障壁においてとしたものであり、そしてこれは数値計算すればわかることなのであるが、共鳴透過が起こるエネルギーと有限井戸型ポテンシャルのエネルギー固有値は、右辺にがかかるかどうかの違いしかないのでほぼ一致する。つまり、入射粒子のエネルギーがポテンシャル障壁のはざまによって作られた疑似的な井戸型ポテンシャルが作っているであろうエネルギー固有状態のエネルギーと一致しているときに共鳴が起こり、入射粒子が減衰されることなく透過するのである。
終わりに
二重ポテンシャル障壁において共鳴透過が起こることの原因は、ポテンシャル障壁の間にできた疑似的な井戸型ポテンシャルの中でエネルギー固有状態のようなものが形成されていることであったことを考えると、ポテンシャル障壁の数を3つ、4つと増やしていっても同じような現象が起きそうである。けど、転送行列の計算がめっちゃ面倒そうなので、こういう行列が簡単にできる方法があれば誰か教えてください。
量子力学の演義でやらされるポテンシャル障壁の散乱にもう一つ障壁を付け加えるだけで変な現象が起きたりするのは結構面白いとおもう。その分計算も結構面倒になるけど。。。
*1:直観的にはかなり気持ち悪い現象だと思う。
統計力学における形式的なアンサンブルの作り方
統計力学においてカノニカルアンサンブルやグランドカノニカルアンサンブル、あるいはアンサンブルなどを構成する方法を紹介する。
なお、この記事は統計力学のはじめの一歩 2020年版4章などの内容をまとめ直したものである。
統計力学の復習
Boltzmann原理によれば系のマクロなエントロピーとエネルギーがの間にあるような状態の数の間には
の関係があるのであった。
またこのBoltzman原理は状態密度を用いれば
と書くこともできる。
ミクロカノニカルアンサンブルから話を始めよう。
ミクロカノニカルアンサンブルにおいて、エネルギーがの間にある状態である確率は、
となるのであった。
そして、導出は省略する*1が温度の熱浴と対象系を全系としてミクロカノニカルアンサンブルを適用してやれば、対象系のエネルギーがであるような確率は、
で与えられる。ここには分配関数と呼ばれ、規格化定数のようなものである。そしてそれは、
で計算される。
この他にもグランドカノニカルアンサンブルなども同様に化学ポテンシャルと温度の熱浴を考え、それと対象系を含めたものにミクロカノニカルアンサンブルを適用すれば得ることができる。
各アンサンブルについてこのように物理的な状況を考え、それに対してミクロカノニカルアンサンブルを適用し新たなアンサンブルを得てもよいが、もっと形式的に新たなアンサンブルを得る方法がないか考えてみよう。
分配関数について
カノニカルアンサンブルにおいて規格化定数のような役割を果たす分配関数についてもう少し考えてみる。分配関数とは、
であった。この積分を計算することを考える。まず、Boltzmann原理を用いて状態密度をエントロピーを用いて書き換えると、
となる。そして、被積分関数が指数関数であることから鞍点法を用いて積分を評価しよう。
であるような点をと書くことにすると、被積分関数はこのの近くでのみ大きな値を持つのでこの付近でのみ積分を考えればよい。すると2次までの近似で被積分関数は、
と近似される。積分されるのは2次の項であり、それはGauss積分で計算でき、結局
となる。
さて、この分配関数の対数をとれば係数の部分はなので無視できることを考えると、
となり、Mathieu関数あるいは、Helmholtz自由エネルギーが得られる。
これらは、ミクロカノニカルアンサンブルにおけるBoltzman原理のカノニカルアンサンブルにおける対応物と考えられる。
つまり、ミクロカノニカルアンサンブルにおいて規格化定数のような役割を果たしていた状態密度がBoltzman原理によってマクロ系のエントロピーと結びついていたのと同様に、カノニカルアンサンブルでは規格化定数にあたる分配関数がマクロ系のMathieu関数(あるいはHelmholtz自由エネルギー)と結びついていることがわかる。
では、Mathieu関数(あるいはHelmholtz自由エネルギー)とは熱力学においてどのような量であっただろうか。
エントロピーの微分を思い出すと、
なのであった。そして、これをについてLugendre変換したものがMathieu関数(あるいはHelmholtz自由エネルギー)である。
ここで、分配関数の定義を見直せば、これは状態密度のを変数とするLaplace変換になっている。即ち、熱力学においてLugendre変換によって各熱力学関数が結びついていたのと対応して、統計力学では各アンサンブルがLpalace変換によって結びつくことになる。
形式的にミクロカノニカルアンサンブルからカノニカルアンサンブルを得る方法
先の過程を逆にたどってみよう。
まず第一にエントロピーからLugendre変換で結びつく何らかの熱力学関数(今の場合はMathieu関数)がミクロな系の量とBoltzman原理のように結びついてほしい。つまり、
という関係にあってほしい。ここではまだわからないミクロカノニカルアンサンブルにおけるにあたる量であり、はMathieu関数である。
次にこのようなを構成することを考える。上の式を変形すると、
となる。
さて、Mathieu関数においてはLugendre変換からによって定まっている*2ので、上の式におけるもこの関係を満たしていなければならない。すると、
というならば鞍点法を用いればLugendre変換から決まるに関する条件を満たしてなおかつ対数を取った時に熱力学関数が出るような規格化定数になる。そして、最後にBoltzman原理を用いれば
となり、これによってエネルギーがの時の相対確率がと決まる。
つまり、分配関数を構成したければ、エントロピーからLugendre変換したい変数についてLugendre変換後の自然な変数をBoltzman定数で割ったもの(つまり逆温度)を変数とするLaplace変換すればよいことがわかる:
こうしてMathieu関数とBoltzman原理からカノニカルアンサンブルを構成できた。
以上の手続きを形式的にまとめると次のようにまとめられる。
- エントロピーから使いたい、即ちBoltzman原理のようにそれによってミクロ系とマクロ系の対応を与える熱力学関数をLugendre変換によって求める。
- 状態密度をLugendre変換によって変えた後の自然な変数(をBoltzman定数で割ったもの)を用いてLaplace変換し、それを新たなアンサンブルにおける分配関数として定義する。
- 等重率の仮定より分配関数からある状態の確率を決める。
グランドカノニカルアンサンブル
形式的にアンサンブルを構成する方法の応用としてミクロカノニカルアンサンブルからグランドカノニカルアンサンブルを求めてみる。
エントロピーについて、その微分は、
であった。
まず、エントロピーの変数のとについてLugendre変換し、新たな熱力学関数を作る。
ここで、とについては、Lugendre変換からとによって決まっている。
この熱力学関数 のLugendre変換後の自然な変数について、がに、がになっているので、変換後の自然な変数をBoltzman定数 で割ったもの、つまりとを変数とするLaplace変換を行えばよい。つまり、
であり、このがグランドカノニカルアンサンブルにおける大分配関数である。
そして、を大分配関数としたことから、エネルギーがであり、粒子数がであるような状態の出現確率が、
と決まる。
こうしてBoltzman原理とミクロカノニカルアンサンブルからグランドカノニカルアンサンブルが構成することができた。
なお、構成方法からわかるようにグランドカノニカルアンサンブルにおけるBoltzman原理の対応物は、
である。ここで、はグランドポテンシャル。
T-pアンサンブル
アンサンブルについても同様に導出してみよう。
まず、エントロピーについて、とについてLugendre変換し、新たな熱力学関数*3を構成する:
ここで、とについては、Lugendre変換からとによって決まっている。
この熱力学関数のLugendre変換後の自然な変数について、がに、がになっているので、変換後の自然な変数をBoltzman定数で割ったもの、つまりとを変数とするLaplace変換を行えばよい。つまり、
であり、このがアンサンブルにおける分配関数である。
そして、を分配関数としたことから、エネルギーがであり、体積がであるような状態の出現確率が、
と決まる。
こうしてBoltzman原理とミクロカノニカルアンサンブルからアンサンブルを構成することができた。
なお、構成方法からわかるようにアンサンブルにおけるBoltzman原理の対応物は、
である。
その他のアンサンブル
同様のことを行えば物理的に使えるかどうかはともかくとして様々なアンサンブルを作ることが可能である。
また、ゴムの熱力学などで使えるアンサンブルについてもこの手法を用いて、それに対応するアンサンブルを構成することができる。この場合エントロピーの微分は、を張力、を変位として
なので、とについてLugendre変換した熱力学関数を考え、それに対応する分配関数を構成してやれば、
となる。そして、エネルギーが、変位がの状態の出現確率は、
となる。