ホーム＞物理・数学のページ＞物理シミュレーション＞モンテ・カルロと水素原子の電子軌道

モンテ・カルロ法による水素原子の電子軌道表示

HYDRO2PY.JPG - 90,012BYTES

シミュレーションの第一弾はモンテ・カルロ法です。

水素原子の電子軌道をあらわす波動関数は、シュレーディンガー方程式を解く練習問題としても、
また量子力学の実際の物理的な現象への適用、量子化学分野への応用といった点でも重要かつ有名な例です。
ですが、これは一般に動径関数（距離の関数）と球面調和関数（方向の関数）の積の形をしており、
なかなかその空間的な分布を把握することは難しいものです。
さらに化学分野で重要な混成軌道となると数式からイメージをつかむのはかなり難題になってきます。

そこで、第1回は波動関数の分布を視覚的に表現する方法を、DirectXプログラミングの練習も兼ねて作ってみたいと思います。

 

1.波動関数と確率密度分布

　量子力学についてはそのうち自分なりにまとめてみたいと思いますが、ここでは簡単に「波動関数とはなんぞや？」ということについて
簡単にお話しておきたいと思います。
　20世紀初頭の物理学において、光は波動としてだけでなく、粒子としての振る舞いも持つ、という「光量子仮説」から量子力学は始まりました。
　「波動だと思っていた光が粒子なら、粒子だと思っていた電子が波動として振舞うこともあるのでは？」ということから、電子の波動性、いわゆる「ド・ブロイ波」の概念がでてきます。
　波動として振舞う、というのは具体的には「回折」や「干渉」という現象を起こすということですが、これはその後実験的に確かめられ、
物質粒子も波動性をもつことが分かりました。これを「波動・粒子二重性」といいます。
　その物質の波がどのような形になっているかを決める方程式がシュレーディンガー方程式です。
　さて、物質が波だとして、「その波は何が振動しているのか？」という問題があります。
　あくまで粒子としての電子は点粒子です（現在の物理学では、と保留をつけたほうがいいかもしれませんが）。
　電子が雲みたいに広がりをもった連続体（流体）のようなものと考えるのは違います。特に今回の方法は冒頭のスクリーンショットにもあるように、電子の雲ようにみえますから、強く言っておきます。
　では何なのか。この振動している何かを量ψ（プサイ）であらわします。これは一般的に複素数で、たとえば干渉は、このψの和によって表すことができるような量です。現在の物理学でもっとも一般的な解釈は、このψの絶対値の２乗、｜ψ｜2がその点において電子が観測される「確率」を表すというものです。これを「確率解釈」といい場の量子論などでも重要な意味をもつ考え方です。

2.電子軌道の波動関数

　原子が正の電荷を持つ部分と電子で構成されるということは分かっていましたが（電子は19世紀末には発見されていましたが、原子内の正電荷がどのように分布しているかはラザフォードの実験(1911-1913)による）、それが実際どういう形なのか、ということが量子力学以前は大問題でした。
　それは前期量子論とも言われるボーアの原子モデルである程度説明がつきましたが、その具体的な軌道についてはシュレーディンガー方程式による解が解決を与えました。具体的な解についてはWikipediaの項目を参照してください（数式書くの面倒で・・・すいません）。
　電子の軌道を決めているのは３つのパラメーター（量子数といいます）です。
　一つ目は主量子数nといい、電場も磁場もない状態ではエネルギーは主量子数だけできまります。これについてはn=1,2,3と数字で表します。
　２つ目は方向量子数lといい、感覚的には軌道が楕円形に偏ることに相当し、物理量としては角運動量をもつことに対応します。これについてはl=0から順にs,p,d,fとアルファベットで表します。
　３つ目は磁気量子数mといい、磁場に対する応答性が変わります。これについては関数の形から座標軸の名前をつけます。
　これの間には制約があって、|l|<n-1、|m|≦lという条件がつきます。
　つまり、最初の軌道は1s、次が2s、2p、その次が3s、3p、3dとなっていきます。そしてｐ軌道は３つ(m=-1,0,-1)、ｄ軌道は５つ(m=-2,-1,0,1,2)の異なる形をとるわけです。
　これによって元素周期の説明とかもできるのですが、スピンの説明もしないといけないのでここでは割愛します。
　水素に限らず、原子内の電子に関してほぼ同じ波動関数を用いて考えることができます。

3.モンテ・カルロ法とは何か？

　さて、水素原子内の電子について波動関数を求めることができたわけですが、これが具体的にどういう分布なのか、ということは式をみただけではぱっとは分かりにくいです。
　そこで、今回は「モンテ・カルロ法」という方法をもちいてこれを視覚化してみたいと思います。
　モンテ・カルロ法は確率的な問題について乱数を用いて計算的に解く方法です。
　今回の問題では、アルゴリズムを書けば以下のようになります。

psi(x,y,z) : 波動関数

1. 乱数で適当な範囲の(x,y,z)と0から1の実数pを生成
2. pp = psi(x,y,z)^2を計算
3. p < pp ならばその点は観測されるので配列に格納
   そうでなければ、その点を捨てる

これを20000個の点が生成されるまで繰り返す

つまり、20000回の観測を状態に影響を与えることなく（物理的には無理なわけですが）繰り返し、
観測されたところに点を打つことで表示するわけです。

4.いざプログラム

　プログラム自体は簡単です。
(最初のソースでa=のところ、不等号が逆になっていました）

ご要望がありましたので現状のソースプログラムをダウンロード可能にしておきます。
いろいろと変更したい点があるのですが、とりあえずのバージョンということで。

ソースコードのダウンロード：DxHydrogen.zip

実行ファイルのダウンロード：ShowProb.zip

実行ファイルのほうは、計算を別にして点のデータにしたものを表示するだけになっています。
計算にはけっこう時間がかかるので、実際にみてぐるぐる動かしてみたいという場合にはこれで十分だと考えました。

    srand((unsigned)time(NULL));
    float    x,y,z, p , pp ,a ,sum = 0.0f;
    const float L = 10.0f;
    float maxpp = 0.0f;
    int nn = 0;

    for( int i=0; i<N ; ++i)
    {
        do {

            x = L * (1.0f - 2.0f * (float)rand() / RAND_MAX);
            y = L * (1.0f - 2.0f * (float)rand() / RAND_MAX);
            z = L * (1.0f - 2.0f * (float)rand() / RAND_MAX);
            p = (float)rand() / RAND_MAX;

            pp = psi(x,y,z);
            a = (pp>0) - (pp<0);
            pp *= pp;
            nn++;
            sum += pp;
        } while( p >= pp );

        vertices[i].Pos.x = x;
        vertices[i].Pos.y = y;
        vertices[i].Pos.z = z;

        vertices[i].Col.r = a > 0.0f ? 0.8f : 0.0f;
        vertices[i].Col.b = 0.8f;
        vertices[i].Col.g = a < 0.0f ? 0.8f : 0.0f;
        vertices[i].Col.a = 1.0f;

        ofs << std::setw(6) << i << ": (" << x << "," << y << "," << z << ") : p = " << p << "/pp = " << pp << std::endl;

        if ( maxpp < pp)
            maxpp = pp;


    ofs << "maxpp = " << maxpp << std::endl;
    ofs << "nn = " << nn << ", ave = " << sum / nn << std::endl;
    ofs.close();

変数の説明

変数名 説明
x,y,z 乱数で発生した座標
p 乱数で発生した確率用の数値
pp 波動関数の絶対値の2乗
a 波動関数の符号
vertices[] VertexBuffer用配列
nn 試行回数
maxpp ppの最大値
L 発生させる座標の範囲(-L〜L)
ofs ログファイル用ofstream

内容は難しくないと思うので説明を省きます。
確率が高いところほど点が生成されやすく、点が密にでます。
確率が低いところは点がまばらになり、点の密度によって波動関数の大きさをみることができる、というわけです。
また波動関数の符号によって色を変えています。
あとはverticesをDirectXのVertexBufferに登録し、プリミティブトポロジーをD3D10_PRIMITIVE_TOPOLOGY_POINTLISTにして描画するだけです。

実際は、カメラ機能を実装するのに苦労したりしましたが、その辺についてはDirectXの記事に反映させるつもりです。

もうちょっと機能がかたまったらソースコードを公開する予定です。

5.いざ実行！・・・乱数発生の罠

さて、実際に動かしてみると挙動が変です。

HYDROERR.JPG - 115,744BYTES

なにか筋が入ってます。
最初は波動関数の計算かDirectXの座標変換行列の設定を間違ったかと思ったのですがちがいました。

で、半日ぐらいあれこれいじっても治らず、ふと思い出したことが。

「線形合同法の乱数でベクトルを生成するときは気をつけなくてはならない」ということです。
Visual C++のrand()関数は線形合同法をつかっています。
これは前の乱数から1次漸化式と剰余（合同式）を使って次の乱数を発生させるものです。
そうするとベクトルのように組み合わせたときにある傾向がでてしまうのです。

そこで乱数を自前の乱数クラス（足し算法）を作ってやり直してみました。

すると・・・解決！こんなにはっきり出てしまうものとは思いませんでした。

 

6.これが波動関数の姿だ

では、このプログラムをつかった波動関数の画像をいくつかお見せしましょう。

まずは2py軌道から

HYDRO2PY.JPG - 90,012BYTES

3d3z2-r2軌道
斜め上から見た図。赤い部分（位相が正）がｚ軸で、それを取り囲むように負の位相部分があります。

HYDRO3D3RZ.JPG - 54,748BYTES

真横から見た図。

HYDRO3D3RZ_2.JPG - 79,329BYTES

3dzx軌道
真上（ｙ軸方向）から見た図。ｘｙ平面とｙｚ平面で波動関数が0になっていることが確認できる。

HYDRO3DZX.JPG - 102,498BYTES

斜めから見た図。ｙ方向の広がりをみることができる。

HYDROG3DZX_2.JPG - 117,076BYTES

ｓｐ混成軌道　思ったより広がっているというか、きのこみたいです。

HYDRO2SP.JPG - 62,800BYTES

実際にマウスでぐるぐる回して眺めてみるともっとよく分かります。
混成軌道についての注意点とか、そもそも波動関数を視覚化するのはどういう方法がほかに考えられるとか、落ちている情報とか、いろいろと書きたいことはあるのですが、とりあえず実行ファイルをダウンロードして眺めてみてください。

7. 今後の展開

今回は波動関数の確率密度としての側面をとりあげたわけですが、
複素数としての波動関数には位相の部分があり、それはこの方法でみることはできません。
そして位相を微分したときにエネルギーや運動量といった物理量が出てくるので
これは今後の課題です。

また、このモンテ・カルロ法は時間がかかります。
これこそCUDAの出番だと思いますので、そのうちCUDA版も作ってみたいと思います。

＜前のページへ＞