processingの備忘録 -電子軌道-

こんにちは。
今回は「電子軌道（原子軌道）」を描画してみようというお話です。

電子軌道とは

大学化学の始まりと言えば、おそらく電子軌道の話かと思います。
高校化学までではK殻、L殻、M殻、N殻…という「電子殻」があり、そこにそれぞれ2個、8個、18個、32個…と電子が満たされていくという、いわゆる「ボーアモデル」的な考え方でした。しかし大学化学では、軌道は「量子数」なるものをいくつかもち、それの組み合わせで「波動関数」が決まり、それによってs軌道、p軌道、d軌道、f軌道…など種類があらわれるという考え方です。
さらに、s軌道は1つ、p軌道は3つ、d軌道は5つ、f軌道は7つ…とセットになっていて、1つの軌道に電子が2個ずつ入ることができます。高校化学での「殻」に入る電子数はこれにより説明でき、K殻はs軌道のみ（1s軌道）のため2個で閉殻、L殻はs軌道とp軌道（2s軌道、2p軌道）があるため$2+2 \times 3$の8個で閉殻、M殻はs軌道、p軌道、d軌道（3s軌道、3p軌道、3d軌道）があるため$2+2 \times 3+2 \times 5$の18個で閉殻…という形になります。
（ちなみにアイキャッチにあるのは本ブログの方法で描画したd軌道5つです）

例えば2p軌道であれば3種類あり、その1つは以下のような波動関数で表されます。

これは「極座標」で表されているため、変数として$r, \theta, \phi$をもちます。
極座標は図のようになっています。

極座標と直交座標

直交座標は$x, y, z$で表される座標系で、以下のように三角関数を使って相互に変換が可能です。
$$ \begin{align} x &= r ~ \mathrm{sin} ~ \theta ~ \mathrm{cos} ~ \phi \\ y &= r ~ \mathrm{sin} ~ \theta ~ \mathrm{sin} ~ \phi \\ z &= r ~ \mathrm{cos} ~ \theta \end{align} $$

今回のブログでは、簡単のため角度に関する$\theta, \phi$の部分だけを取り出して描画しようと考えているので、描画の結果は電子軌道のおおまかな形のみを表していて、$r$に関する部分は省いているので大きさはわかりません。
（正確には電子軌道に形はなく（「電子雲」というやつ）、「電子が見いだされる確率がこのあたりで高いね」という部分が面として表されていると思ってください。もっというと、波動関数の2乗が存在確率になるので今回の方法で表されるのは存在確率でもありませんが…。また実際は1sより2sの方が大きいなど、大きさも異なります）

ここに書いたことは大雑把な説明なので、詳しく知りたければ本を探してみてください。

電子軌道の描画

それでは描画していきます。
以下、3Dでの描画は以前つくったグリグリ動かせるコードを使っているので、ゴチャゴチャいろいろと書いてあります。

とりあえず、軸を描いてみましょう。（都合上、右手系になるようにz座標に「-」をつけてあります。）

float rotX, rotY;
float protX, protY;
float xc, yc, zc;
float pxc, pyc, pzc;
float s;

void setup(){
    size(500, 500, P3D);
    
    xc = 0;
    yc = 0;
    zc = 500;
    s = 1;
}

void draw(){
    background(255);
    
    camera(0, 0, 1000, xc, yc, zc, 1, 0, 0);
    rotateX(rotX);
    rotateY(rotY);
    scale(s);
    
    //x軸
    stroke(255, 0, 0);
    line(0, 0, 0, 500, 0, 0);
    //y軸
    stroke(0, 255, 0);
    line(0, 0, 0, 0, 500, 0);
    //z軸
    stroke(0, 0, 255);
    line(0, 0, 0, 0, 0, -500);
}

void mouseDragged(){
    if(mouseButton == LEFT){
        if(mouseX-pmouseX >= 0){
            rotX = map(mouseX-pmouseX, 0, width, 0, 2*PI);
        }
        if(mouseX-pmouseX <= 0){
            rotX = map(mouseX-pmouseX, -width, 0, -2*PI, 0);
        }
        rotX += protX;
        protX = rotX;
        
        if(mouseY-pmouseY >= 0){
            rotY = map(mouseY-pmouseY, 0, height, 0, 2*PI);
        }
        if(mouseY-pmouseY <= 0){
            rotY = map(mouseY-pmouseY, -height, 0, -2*PI, 0);
        }
        rotY += protY;
        protY = rotY;
    }
    
    if(mouseButton == CENTER){
        xc += pmouseY-mouseY;
        pxc = xc;
        yc += mouseX-pmouseX;
        pyc = yc;
    }
}

void mouseWheel(MouseEvent e){
    float mw = e.getCount();
    if(mw == 1){
        s *= 0.9;
    }
    if(mw == -1){
        s *= 1.1;
    }   
}

例として、先に挙げた$\mathrm p_x$軌道について描画してみます。

float rotX, rotY;
float protX, protY;
float xc, yc, zc;
float pxc, pyc, pzc;
float s;

float r;
float[][] x = new float[360][180];
float[][] y = new float[360][180];
float[][] z = new float[360][180];

void setup(){
    size(500, 500, P3D);
    
    xc = 0;
    yc = 0;
    zc = 500;
    s = 1;
    
    r = 100;
}

void draw(){
    background(255);
    
    camera(0, 0, 1000, xc, yc, zc, 1, 0, 0);
    rotateX(rotX);
    rotateY(rotY);
    scale(s);
    
    stroke(255, 0, 0);
    line(0, 0, 0, 500, 0, 0);
    stroke(0, 255, 0);
    line(0, 0, 0, 0, 500, 0);
    stroke(0, 0, 255);
    line(0, 0, 0, 0, 0, -500);
    
    for(int theta = 0; theta < 360; theta ++){ //修正箇所(tをthetaに)
    for(int phi = 0; phi < 180; phi ++){ //修正箇所(pをphiに)
        float t = radians(theta); //修正箇所(radians()の処理を追加)
        float p = radians(phi); //修正箇所(radians()の処理を追加)
        
        float Y;
        Y = sin(t)*cos(p);
      
        x[theta][phi] = abs(r*Y) * sin(t) * cos(p); //修正箇所([]内のt、pをそれぞれtheta、phiに)
        y[theta][phi] = abs(r*Y) * sin(t) * sin(p); //修正箇所([]内のt、pをそれぞれtheta、phiに)
        z[theta][phi] = abs(r*Y) * cos(t); //修正箇所([]内のt、pをそれぞれtheta、phiに)
        
        if(r*Y > 0) stroke(0, 0, 255);
        if(r*Y < 0) stroke(255, 0, 0);
        point(x[theta][phi], y[theta][phi], -z[theta][phi]); //修正箇所(t、pをそれぞれtheta、phiに)
    }
    }
}

void mouseDragged(){
    if(mouseButton == LEFT){
        if(mouseX-pmouseX >= 0){
            rotX = map(mouseX-pmouseX, 0, width, 0, 2*PI);
        }
        if(mouseX-pmouseX <= 0){
            rotX = map(mouseX-pmouseX, -width, 0, -2*PI, 0);
        }
        rotX += protX;
        protX = rotX;
        
        if(mouseY-pmouseY >= 0){
            rotY = map(mouseY-pmouseY, 0, height, 0, 2*PI);
        }
        if(mouseY-pmouseY <= 0){
            rotY = map(mouseY-pmouseY, -height, 0, -2*PI, 0);
        }
        rotY += protY;
        protY = rotY;
    }
    
    if(mouseButton == CENTER){
        xc += pmouseY-mouseY;
        pxc = xc;
        yc += mouseX-pmouseX;
        pyc = yc;
    }
}

void mouseWheel(MouseEvent e){
    float mw = e.getCount();
    if(mw == 1){
        s *= 0.9;
    }
    if(mw == -1){
        s *= 1.1;
    }   
}

座標については、二次元の配列（x[t][p]の形、tはtheta$~\theta$、pはphi$~\phi$）を使うことで$\theta, \phi$についてに表せるようにしていて、$\theta$は0~360、$\phi$は0~180の値をとるので配列の数はそれぞれ360と180にしています。
極座標から直交座標に変換するためにabs(rY)の後に三角関数がついています。（abs()は絶対値を返す関数）
大きさ$r$に関する部分は考えていないので100に固定しています。角度$\theta, \phi$に関する部分は$\mathrm{sin} ~ \theta ~ \mathrm{cos} ~ \phi$となっています。 何故かわかりませんが$\theta, \phi$に関する部分はまとめて$Y(\theta, \phi)$（球面調和関数）と書かれることが多いのでそうしていて、rYはある角度での原点からの距離を表すものなので、実際に使う際は絶対値をとっています。
波動関数（ここではrY）が正となるときには青で、負となるときには赤で表すことが多いので、if()で場合分けして塗り分けています。

追記（2021/11/17）
角度を弧度法に直すのを忘れていました！（通りで図が何か帯みたくなるなと思った…）
for文が始まってすぐにradians()でtheta、phiを弧度法に変換するようにし、配列の引数についてtをthetaに、pをphiに変更しました。
（それに伴って、アイキャッチも更新しました。）
以上

実行してみると、
x軸方向に広がりを持っていることがわかります。

ちなみに、波動関数は2乗すると電子の存在確率になるということなのでやってみると、
少し引き伸ばされたような形になりました。

おわりに

波動関数は訳のわからない形をしていますが計算によって求められたものなので、やろうと思えばコードに組み込んで「この条件のときはこんな波動関数でこんな形」という風に機械に解いてもらって描画することも可能です。いわゆる「Schrödinger方程式を解く」ということになるわけですが、自分はそこまで理解できているわけではないので触れないことにします。
ですが、文献を漁れば書いてあることではあるので、rやYを対応するものに置き換えれば大きさも反映した形になるでしょうし、2乗してみればより実際の形に近いものが描画できるかと思います。

というわけで、今回は電子軌道を描画してみました。
ここからいろいろと付け足して遊んでみるのも良いかもしれません。
それでは、また次回。