利用Gnuplot軟體繪製似氫原子的徑向分佈函數（下）

利用Gnuplot軟體繪製似氫原子的徑向分佈函數（下）(Hydrogen-like radial distribution Plots Using Gnuplot (II))
國立臺灣師範大學化學系兼任教授 邱智宏

連結：利用Gnuplot軟體繪製似氫原子的徑向分佈函數（上）

接下來，探究第一個疑問，一般會出現這樣的矛盾，主要是將機率密度和徑向分佈函數 (radial distribution function) 弄混，前者是空間某一特定點電子出現的機率密度，而密度要乘上體積才是距原子核某特定距離電子出現機率的大小。

如果要探究距離原子核從 $$r$$ 到 $$r+dr$$，$$\theta$$ 到 $$\theta+d\theta$$，$$\phi$$ 到 $$\phi+d\phi$$ 之間的單位體積內找到電子的機率有多大？則必需計算其間體積的變化量再乘上該處的電子機率密度 $$(|\varphi|^2)$$。由圖五可看出極座標單位體積 $$(dV)$$ 的變化量為：$$dV=r\sin\theta~d\phi\times rd\theta\times dr=r^2\sin\theta~d\theta~d\phi~dr$$

圖五、極座標中單位體積 $$(dV)$$ 的表示法示意圖（圖片來源：參考資料6）

由於 $$1s$$、$$2s$$ 和 $$3s$$ 軌域和 $$\theta$$、$$\phi$$ 無關，因此可單看徑度變化的影響如圖六所示，探索在從 $$r$$ 到 $$r+dr$$ 間圓殼 (spherical shell) 內找到電子的機率有多少？

圖六、$$1s$$ 軌域徑度從 $$r$$ 到 $$r+dr$$ 間圓殼體積的示意圖（圖片來源：http://images.slideplayer.com/23/6601420/slides/slide_44.jpg）

若將 $$dV$$ 乘上機率密度並將 $$\theta$$、$$\phi$$ 積分，即能檢驗徑度變化對電子出現機率的影響：

$$\begin{array}{ll}\displaystyle |\varphi|^2r^2dr\int^{2\pi}_{0}d\phi\int^{\pi}_{0}\sin\theta~d\theta &=|\varphi|^2r^2dr\times(2\pi-0)\times(-\cos\pi+\cos 0)\\&=|\varphi|^2r^2dr\times(2\pi)\times 2\\&=4\pi r^2|\varphi|^2dr\end{array}$$

若將上式的 $$4\pi r^2|\varphi|^2$$（一般稱為徑向分佈函數）對 $$r$$ 作圖，則可觀察氫原子 $$1s$$、$$2s$$ 和 $$3s$$ 軌域電子在不同的徑向殼層中，其出現機率的總和分佈情形。

此時若利用 gnuplot 繪圖，只要將圖一的程式集稍加修加即可。首先在圖一中，原先已經輸入過 $$1s$$、$$2s$$ 和 $$3s$$ 的徑向函數分別為 P1s(x)、P2s(x) 和 P3s(x)，因此只要將最後二行的 Phi1s(x)、Phi1s(x) 和 Phi1s(x) 換成前述的徑向分佈函數，並重新設定 y 軸的範圍，將 set ylabel [0.,0.6] 加入圖一中的文字檔指令集中，再執行指令集便能得到圖七。由圖中可看出，紅色 $$1s$$ 線形出現 $$1$$ 個極大值，其位置最接近原子核，藍色 $$2s$$ 線形有 $$2$$ 個極值，後者的極值較大，即機率值較大。棕色的 $$3s$$ 則有 $$3$$ 個極值，也是最後一個極值較大。

圖七中會出現極值的原因，若以 $$1s$$ 為例，$$r$$ 值愈大時，其球殼層的體積也愈大，雖然在原子核的位置電子密度最大，但其體積為 $$0$$，因此在此處找到電子的機率仍為 $$0$$。隨著 $$r$$ 的增大，機率密度儘管下降，球殼層的體積卻增大，因此找到電子的機率會逐漸上升，但前者下降的速率比後者大，兩者相權的結果，極大值於焉產生。

圖七、利用 gnuplot 繪製氫原子 $$s$$ 軌域電子的徑向分佈函數圖（圖片來源：作者製作）

三、利用 gnuplot 求出 $$2s$$ 軌域區間的積分值

Gnuplot 軟體除了能畫出各類函數的專業圖形以外，還有許多其他功能，在此特別介紹其具有類似 C 語言的程式功能，因此使得該軟體如虎添翼，可以處理更多技術性的問題。例如圖七中，如果我們想要了解 $$2s$$ 軌域，第一個小山丘電子出現的總機率究竟多少？即 $$r$$ 從 $$0$$ 到 $$2a_0$$ 的範圍間，圖形下的總面積為何？相當於是算出下列式子的積分值：

 $$\displaystyle\int^{2a_0/Z}_{0}4\pi|\varphi_{2s}|^2r^2dr=\frac{Z^3}{8a^3_a}\int^{2a_0/Z}_{0}(2-\frac{Zr}{a_0})^2r^2e^{-\frac{Zr}{a_0}}dr$$

上式的確有精確解，可利用部分積分 (integral by parts) 的方法求解，有興趣的讀者可以參考高瞻平台的文章《似氫原子2s軌域的解析》，其值為 $$0.053$$，約佔整體電子出現機率的 $$5.3\%$$。在這裏擬用數值解的方式，利用 gnuplot 的程式指令求解。

積分的數值解方式很多，在這裏使用辛普森法 (Simpson method)，其求法如（式-4）所列，將上下區間 $$a=0$$ 和 $$b=2$$ 間分成數個等距的區域 $$(n)$$，每個區域的間距為 $$h=0.01$$，將 $$x_0=a$$、$$x_1=a+h$$、$$x_2=a+2h$$ 代入下式的第一個中括弧中，然後將 $$x_2$$、$$x_3=x_2+h$$、$$x_4=x_2+2h$$ 代入第二個中括弧中，依序將所有偶數點均帶入下式，直到 $$x_n=b$$ 為止，最後將所有中括弧的和相加後，乘於 $$h$$ 除於 $$3$$ 即能得到答案。

$$\begin{multline*}\int^b_{a}f(x)\approx\frac{h}{3}\{[f(x_0)+4f(x_1)+f(x_2)]+[f(x_2)+4f(x_3)+f(x_4)]+\\ \cdots+[f(x_{n-2})+4f(x_{n-1})+f(x_n)]\}\end{multline*}$$   （式子-4）

在 gnuplot 軟體中，應如何將上述的程序寫成程式集？首先將圖一的最後二行改成如下：

plot P2s(x) lw 2 lc “blue” title “2s”

因為我們只要畫 $$2s$$ 的徑向分佈函數，當然標題、y 軸的範圍也要一併修改，接下來將下列文字檔一併加入圖一的檔案中即可。

圖八、利用 gnuplot 指令集，以辛普森數值法算出氫原子 $$2s$$ 軌域在第一個節點前，電子出現的機率總和（圖片來源：作者製作）

圖八的第 2 行連續設定 $$4$$ 個變數，即積分的上限、下限、區間值 (h=0.01) 及初始面積 s=0，第 3 行 n=floor((b-a)/h)，設定 n 個小區間，floor() 函數乃取整數值的意思，譬如 floor(7/2) 所得的值為 3。接下來的這個區塊

do for [i=0:n-2] {
xi=a+h*i
if (i%2==0){ s=s+P2s(xi)+4*P2s(xi+h)+P2s(xi+2*h)}
}

do for 代表程式會依照 [i=0:n-2] 的條件來判斷是否繼續執行大括弧 {…} 內的指令，i 變數的起始值為 0，每次執行 {…} 內的指令一次，則 i 增加 1，直到 i=n-2 為止。

i 每改變一次，xi 也改變一次，讓 xi=a+h*i，此時進入 if 指令，若 i%2==0 成真時才會進行紅色 {…} 內的指令。i%2 是將 i 除以 2 後，所得的餘數，當 i 為偶數時，i%2 的餘數為 0，i 為奇數時餘數不為 0，此時作邏輯判斷，== 表示判斷前後兩數是否相等。i 為偶數時經 i%2 處理後餘數自然為 0，因此邏輯判斷會真，接下來便會執 {…} 內的指令，若為奇數則不處理。

處理時會依辛普森法的（式-4），把一小區域、一小區域的面積陸續累加起來，即 s=s+P2s(xi)+4*P2s(xi+h)+P2s(xi+2*h)，直至完成 do for 指令。此時再把總面積除於 3，再乘以 h，即 s=s*h/3.。此時區間的積分值便告完成，並存在 s 變數中。圖八接下來的指令如下：

set label sprintf(“Integral [0:2] = %10.8f”,s) at 3,0.04 font “Times,12” textcolor “red”

只是利用 sprintf 印出一個標題 (label)，內容為 Integral [0:2] = s，但 s 是一個變數，其值已由前述各行算出，%10.8f 代表印出的標題需要一個變數，在 s，此設定為 s 並以具有 8 位小數的浮點數方式顯示，另外顯示在座標軸為 3,0.04 的位置，字體為 Times 的 12 號字形，字的顏色為紅色。

最後一組 do for 指令是畫出 11 條垂直没有箭頭 (nohead arrow) 的藍線，x 軸每間隔 0.2 畫一條。由於前面已經使用過 plot 指令，此時由於要另外要附加一些繪圖，因此要使用 replot指令方能完成任務，繪出的圖形詳如圖九。

圖九、利用 gnuplot 的指令，計算電子在 $$2s$$ 軌域第一個小山丘出現的機率（圖片來源：作者製作）

圖九中可看出透過 gnuplot 使用辛普森數值法所求出，$$r$$ 從 $$0$$ 到 $$2a_0$$ 區間的積分值為 $$0.0527$$，和精確解的答案幾乎一樣。當然我們也可以改變不同的區間，改變不同的函數來做運算，但利用 gnuplot 同樣的可以輕易的求出準確的解答。

四、結論

上一篇相關的文章如前言所述，曾利用 gnuplot 軟體，繪製似氫原子各類軌域的立體圖形，本文則介紹其在平面繪圖上的應用，一面複習使用過的指令，一面增添常用及新的指令，尤其是類似程式語言的部分例如 do for、if 等，期能對此軟體的應用能力，有大幅提升的功效。對於熟悉程式語言的讀者，文中的示範說明，雖然是多此一舉，但說不定能激發其創意，並應用在其他不同的領域上也未可知；對於不熟悉程式語言的讀者，也可以按圖索驥，邊讀、邊做、邊學，一定會有意想不到的功效。

本文除了學習 gnuplot 軟體的應用以外，並利用它來繪製似氫原子的波函數、機率密度函數及徑向分佈函數對原子半徑的分佈圖，讓讀者可以親自將教科書中抽象的函數，轉化為具體的圖形，更可以藉由座標軸範圍的改變，局部放大或全面觀察各種圖形中關鍵性的變化。另外，也藉機釐清：為何電子密度最高的原子核，卻是最找不到電子的地方？更能利用該軟體以數值解的方式將任何軌域在不同區間中，電子出現的機率總和給核算出來。

參考文獻

1. gnuplot homepage. http://gnuplot.info/
2. 地圖/統計圖/3d 函數圖/實驗報告圖 — Gnuplot 純畫圖｜”資訊人權貴” 之家。 http://user.frdm.info/ckhung/b/ma/gnuplot.php
3. gnuplotスクリプトの解説｜米澤進吾 ホームページ http://www.ss.scphys.kyoto-u.ac.jp/person/yonezawa/contents/program/gnuplot/index.html
4. 马欢 (2012)。使用 gnuplot 科学作图｜Gnuplot 中文教程。 http://www.phy.fju.edu.tw/files/archive/876_ab57aed9.pdf
5. N. Levine (2008), Physical Chemistry (6th ed.). McGRAW-HILL Book Company. p637~647.
6. Spherical Polar Coordinates – HyperPhysics. http://hyperphysics.phy-astr.gsu.edu/hbase/sphc.html
7. 葉名倉、劉如熹、邱智宏、周芳妃、陳建華、陳偉民（2013 年）高級中學化學選修上冊。南一書局。第14～26頁。
8. Saputra, A., Canaval, L. R., Fadiawati, N., Diawati, C., Setyorini, M., Kadaritna, N., & Kadaryanto, B. (2015). Visualizing Three-Dimensional Hybrid Atomic Orbitals Using Winplot: An Application for Student Self Instruction. Journal of Chemical Education, 92(9), 1557-1558.
9. Chung, W. C. (2013). Three-dimensional atomic orbital plots in the classroom using Winplot. Journal of Chemical Education, 90(8), 1090-1092.
10. Moore, B. G. (2000). Orbital Plots Using Gnuplot. Journal of Chemical Education, 77, 785-789.