モンテ・カルロ法による水素原子の電子軌道ホーム>物理・数学のページ>物理シミュレーション>モンテ・カルロと水素原子の電子軌道モンテ・カルロ法による水素原子の電子軌道表示シミュレーションの第一弾はモンテ・カルロ法です。応用といった点でも重要かつ有名な例です。これは一般に動径関数(距離の関数)と球面調和関数(方向の関数)の積の形をしており、なかなかその空間的な分布を把握することは難しいものです。化学分野で重要な混成軌道となると数式からイメージをつかむのはかなり難題になってきます。そこで、第1回は波動関数の分布を視覚的に表現する方法を、DirectXプログラミングの練習も兼ねて作ってみたいと思います。ついてはそのうち自分なりにまとめてみたいと思いますが、ここでは簡単に「波動関数とはなんぞや?」ということについて簡単にお話しておきたいと思います。20世紀初頭の物理学において、光は波動としてだけでなく、粒子としての振る舞いも持つ、という「光量子仮説」から量子力学は始まりました。「波動だと思っていた光が粒子なら、粒子だと思っていた電子が波動として振舞うこともあるのでは?」ということから、電子の波動性、いわゆる「ド・ブロイ波」の概念がでてきます。物質粒子も波動性をもつことが分かりました。その物質の波がどのような形になっているかを決める方程式がシュレーディンガー方程式です。物質が波だとして、「その波は何が振動しているのか?」という問題があります。粒子としての電子は点粒子です(現在の物理学では、と保留をつけたほうがいいかもしれませんが)。電子が雲みたいに広がりをもった連続体(流体)のようなものと考えるのは違います。今回の方法は冒頭のスクリーンショットにもあるように、電子の雲ようにみえますから、強く言っておきます。この振動している何かを量ψ(プサイ)であらわします。これは一般的に複素数で、たとえば干渉は、このψの和によって表すことができるような量です。現在の物理学でもっとも一般的な解釈は、このψの絶対値の2乗、|ψ|2がその点において電子が観測される「確率」を表すというものです。これを「確率解釈」といい場の量子論などでも重要な意味をもつ考え方です。原子が正の電荷を持つ部分と電子で構成されるということは分かっていましたが(電子は19世紀末には発見されていましたが、原子内の正電荷がどのように分布しているかはラザフォードの実験(1911-1913)による)、それが実際どういう形なのか、ということが量子力学以前は大問題でした。それは前期量子論とも言われるボーアの原子モデルである程度説明がつきましたが、その具体的な軌道についてはシュレーディンガー方程式による解が解決を与えました。具体的な解についてはWikipediaの項目を参照してください(数式書くの面倒で・・・すいません)。電子の軌道を決めているのは3つのパラメーター(量子数といいます)です。一つ目は主量子数nといい、電場も磁場もない状態ではエネルギーは主量子数だけできまります。これについてはn=1,2,3と数字で表します。2つ目は方向量子数lといい、感覚的には軌道が楕円形に偏ることに相当し、物理量としては角運動量をもつことに対応します。これについてはl=0から順にs,p,d,fとアルファベットで表します。3つ目は磁気量子数mといい、磁場に対する応答性が変わります。これについては関数の形から座標軸の名前をつけます。これの間には制約があって、|l|<n-1、|m|≦lという条件がつきます。最初の軌道は1s、次が2s、2p、その次が3s、3p、3dとなっていきます。p軌道は3つ(m=-1,0,-1)、d軌道は5つ(m=-2,-1,0,1,2)の異なる形をとるわけです。これによって元素周期の説明とかもできるのですが、スピンの説明もしないといけないのでここでは割愛します。水素に限らず、原子内の電子に関してほぼ同じ波動関数を用いて考えることができます。3.モンテ・カルロ法とは何か?さて、水素原子内の電子について波動関数を求めることができたわけですが、これが具体的にどういう分布なのか、ということは式をみただけではぱっとは分かりにくいです。そこで、今回は「モンテ・カルロ法」という方法をもちいてこれを視覚化してみたいと思います。モンテ・カルロ法は確率的な問題について乱数を用いて計算的に解く方法です。今回の問題では、アルゴリズムを書けば以下のようになります。観測されたところに点を打つことで表示するわけです。4.いざプログラムプログラム自体は簡単です。(最初のソースでa=のところ、不等号が逆になっていました)ご要望がありましたので現状のソースプログラムをダウンロード可能にしておきます。変更したい点があるのですが、とりあえずのバージョンということで。ソースコードのダウンロード:DxHydrogen.zip実行ファイルのダウンロード:ShowProb.zip実行ファイルのほうは、計算を別にして点のデータにしたものを表示するだけになっています。計算にはけっこう時間がかかるので、実際にみてぐるぐる動かしてみたいという場合にはこれで十分だと考えました。ログファイル用ofstream内容は難しくないと思うので説明を省きます。確率が高いところほど点が生成されやすく、点が密にでます。確率が低いところは点がまばらになり、点の密度によって波動関数の大きさをみることができる、というわけです。波動関数の符号によって色を変えています。あとはverticesをDirectXのVertexBufferに登録し、プリミティブトポロジーをD3D10_PRIMITIVE_TOPOLOGY_POINTLISTにして描画するだけです。実際は、カメラ機能を実装するのに苦労したりしましたが、その辺についてはDirectXの記事に反映させるつもりです。機能がかたまったらソースコードを公開する予定です。5.いざ実行!・・・乱数発生の罠さて、実際に動かしてみると挙動が変です。なにか筋が入ってます。最初は波動関数の計算かDirectXの座標変換行列の設定を間違ったかと思ったのですがちがいました。半日ぐらいあれこれいじっても治らず、ふと思い出したことが。「線形合同法の乱数でベクトルを生成するときは気をつけなくてはならない」ということです。C++のrand()関数は線形合同法をつかっています。これは前の乱数から1次漸化式と剰余(合同式)を使って次の乱数を発生させるものです。するとベクトルのように組み合わせたときにある傾向がでてしまうのです。そこで乱数を自前の乱数クラス(足し算法)を作ってやり直してみました。すると・・・解決!こんなにはっきり出てしまうものとは思いませんでした。6.これが波動関数の姿だでは、このプログラムをつかった波動関数の画像をいくつかお見せしましょう。赤い部分(位相が正)がz軸で、それを取り囲むように負の位相部分があります。真上(y軸方向)から見た図。xy平面とyz平面で波動関数が0になっていることが確認できる。y方向の広がりをみることができる。思ったより広がっているというか、きのこみたいです。実際にマウスでぐるぐる回して眺めてみるともっとよく分かります。混成軌道についての注意点とか、そもそも波動関数を視覚化するのはどういう方法がほかに考えられるとか、落ちている情報とか、いろいろと書きたいことはあるのですが、とりあえず実行ファイルをダウンロードして眺めてみてください。複素数としての波動関数には位相の部分があり、それはこの方法でみることはできません。位相を微分したときにエネルギーや運動量といった物理量が出てくるのでこれは今後の課題です。このモンテ・カルロ法は時間がかかります。これこそCUDAの出番だと思いますので、そのうちCUDA版も作ってみたいと思います。