利用 Gnuplot 軟體繪製似氫原子的軌域形狀（下）

利用 Gnuplot 軟體繪製似氫原子的軌域形狀（下）(Hydrogen-like orbital plots using gnuplot (II))
國立臺灣師範大學化學系兼任教授 邱智宏

連結：利用 Gnuplot 軟體繪製似氫原子的軌域形狀（上）

圖三的文字檔以 fourier.txt 存檔後，可在 gnuplot 視窗的游標列上輸入：load “fourier.txt”，讓 gnuplot 讀取檔案後，一次執行一系列的指令以完成任務，當然讀取檔案的路徑必須正確，若不正確則無法順利執行。當內定路徑無法找到檔案時，可利用視窗中 ChDir 的下拉式選項改變讀取檔案的路徑。

圖三的指令集只是一個示範，初學者肯定看不懂，經過一段時間的練習，或許便能達到此出神入化的境界，其中比較困難的指令先略而不談，挑較簡單的部分說明如下：第 1 行的 reset 指令是將先前定義過的各項設定還原為預設值，這一點在編緝指令集時很重要，不要忘了在第一行加入這個指令，因為使用者經常會忘記先前下過那些指令，如果没有 reset，難保不會出現不可預期的結果。第 2 行的 multiplot，是在同一個視窗（即畫布）中畫出多個 2D 圖形，如果没有設定 multiplot，當下達 plot 指令時只能繪製單個圖形（畫 3D 則使用 splot）。set multiplot layout 2,2，即表示在同一個畫布中畫出二行二列共 $$4$$ 個圖形。

接下來的指令集屬於函數設定及程式控制，程式控制屬於較高階的技巧，我們暫時不會用到，請大家不用擔心，只要經過一段時間的練習，便能駕輕就熟。其大意是做 $$4$$ 次不同階層的傅立葉轉換 (fourier transform)，其結果如圖四。

圖四、使用 multiplot 指令畫出多個 2D 的傅立葉轉換圖形（作者製作）

二、利用 gnuplot 繪製似氫原子的 $$s$$ 軌域

有了上述基本概念以後，我們開始將其應用在似氫原子軌域的畫法。首先需先了解似氫原子的波函數 (wave function)：

$$\varphi(r,\theta,\phi)=N\times Y(\theta,\phi)\times R(r)$$         （式-1）

波函數 $$(\varphi)$$ 由三部分組成，$$N$$ 為正規化係數 (normalization factor)，$$Y(\theta,\phi)$$ 為角度部分 (angular part)，$$R(r)$$ 為徑向部分，一些常見的波函數如表一所列。例如 $$2s$$ 軌或的波函數可表示為

$$\displaystyle\varphi_{2s}(r,\theta,\phi)=\frac{1}{\sqrt{4\pi}}\times 1\times \frac{1}{\sqrt{2}}(\frac{Z}{a_0})^{\frac{3}{2}}\left(1-\frac{Zr}{2a_0}\right)e^{-Zr/2a_0}$$         （式-2）

上式中的 $$Z$$ 為原子序，$$a_0$$ 為波耳半徑，本文為簡化起見，均以 $$Z=1$$ 處理，另外作圖時以 $$r/a_0$$ 作為座標尺度，因此在撰寫指令時 $$Z$$ 和 $$a_0$$ 均可省略不寫。由表一可知波函數的外部形狀只和角度函數有關，徑向函數則和形狀大小及內部結構有關，若僅想了解各類波函數的形狀，可專注於（式-1）的前二項即可。除此之外，（式-1）是以極座標表示，必須將其轉換為大家熟悉的直角座標。

表一、一些常見似氫原子波函數的角度及徑向函數（作者整理）

有關直角座標和極座標的轉換方式，如圖五所示，$$(x,y,z)$$ 轉成 $$(r,\theta,\phi)$$ 可透過下列公式完成：

$$\begin{array}{ll}x&=r\sin\theta\cos\phi\\y&=r\sin\theta\sin\phi\\z&=r\cos\theta\end{array}$$               （式-3）

其中 $$\phi$$ 角為 $$r$$ 在 $$xy$$ 平面上的投影，其繞 $$z$$ 軸旋轉和 $$x$$ 軸間的夾角，$$\theta$$ 角為 $$r$$ 和 $$z$$ 軸間的夾角。

圖五、直角座標和極座標間的關係圖（來源：參考資料6）

有了上述基本資料以後，我們可以開始編寫一些指令集詳如圖六，並以 gnuplot 能讀取的文字檔 (.txt) 形式存檔，例如以檔名為 orbital_ns.txt 的指令集，來練習畫出似氫原子的 $$s$$ 軌域。

圖六各行指令的意義說明如下，首先，以「#」符號開頭的行，代表此行為註解，實際繪圖時並不會執行，「#」也不一定要寫在最前頭，若寫在中間即表示該行在此之後的部分為說明，僅執行前面的指令，當然其有效性不會影響到下一行。

reset 指令很重要已如前述，接下來 set view 60, 30,1.,1.5 指令，是設定圖形畫好以後觀看的視角，第一個數字 $$60$$，代表逆時針翻轉 $$x$$ 軸讓圖形旋轉 $$60$$ 度，第二個數字 $$30$$，則為逆時針旋轉 $$z$$ 軸 $$30$$ 度， 1. 代表圖形不放大，以原尺寸顯示（1. 為實數 (real)，1 則為整數 (integer) 的表示法），而 1.5 則為 $$z$$ 軸放大 $$1.5$$ 倍，當圖形視窗出現時，可以實際使用滑鼠測試、感受一下圖形旋轉不同角度的效果，可以仰視、俯視或任何角度觀看同一圖形。

圖六、繪製似氫原子 $$s$$ 軌域形狀的指令集（作者製作）

Gnuplot 也可以設定變數，例如第 6 行的 box=0.6，即將 box 當 0.6 使用。接下來將 x、y、z 的範圍加以設定，如果没有設定，gnuplot 會自動幫你設定最佳的區間，在此 x、y 的範圍均設定為 [-box:box]，即從 -0.6 到 0.6，至於 z 軸則為不同的設定方法，區間為 -0.4 到 0.4。另外 set ytics -0.6,0.3,0.6 是設定 y 軸的刻度及刻度間距，其設定為 y 軸從 -0.6 開始標刻度，每間隔 0.3 標一個，直到 0.6 為止，如果没有設定則系統將會自動設定，例如 x 軸及 z 軸。再下來 set nokey 代表圖形中不顯示畫出圖形方程式的名稱及圖例。

接下來的第 14 行是設定繪圖參數 (set parametric)，3D 圖形的內定參數為 u 和 v，2D 圖形則為 x，即 gnuplot 是根據 u 和 v 的範圍，將 u 和 v 的值依序代入函數中，一點一點繪出圖形。這裏要特別提醒，u、v 的範圍和座標軸範圍的意義不同，前者是真正繪圖時使用，後者是畫好圖以後，選擇要觀看的範圍，若選擇不當也許只能看到一部分的圖形。我們將 u 的範圍設定為 $$0$$ 到 $$\pi$$ (set urange [0:pi] )，即為 $$\theta$$ 角的區間，y 的範圍設定為 $$-\pi$$ 到 $$\pi$$ (set vrange [-pi:pi] )，即為 $$\phi$$ 角的區間。

接著下面三行函數的設定 (x(u,v)=sin(u)*cos(v)、y(u,v)=sin(u)*sin(v)、z(u,v)=cos(u) )，主要是利用（式-3）將極座標轉換為直角座標。再下來的幾行，有關解析度及圖形顯現方式前面已解釋過，接著直接跳到設定 $$s$$ 軌域的部分：s(u,v)=sqrt(1.0/4.0/pi)，相當於是 $$s$$ 軌域的正規化係數乘於角度函數：$$s(u,v)=\frac{1}{\sqrt{4\pi}}$$，從這個函數可看出，不管 u、v 代入何值 s(u,v) 始終為定值，可見其圖形必為圓球體。

最後開始繪圖，先設定圖形的標題為：s-orbital 即利用下列指令 set title ‘s orbital’，再設定另一個函數 orbp(u,v)，使其等於 s(u,v)，這樣作的目的是後續要畫其他軌域時較為方便並無必然性。

接下來的二行指令較難理解，請耐心領悟一下，posp(u,v)= orbp(u,v)<0 ? 0 : orbp(u,v)，代表先判斷一下 orbp(u,v) 是否小於0，若是則設定 posp(u,v)=0，若否則 posp(u,v)=orbp(u,v)，此舉即保證 posp(u,v) 函數所得的值均為正數或為 0。相同地，negp(u,v)= orbp(u,v)>0 ? 0 : -orbp(u,v) 這一行的意思也是一樣，即保証 negp(u,v) 的值均為正值或為 0。

讀者也許會覺得困惑，s(u,v) 不是為常數嗎？做此判斷是否多此一舉，但是我們若要繪 $$p$$、$$d$$、$$f$$ 軌域時，這二行指令就顯得格外重要。最終以 splot 開始繪圖，必須給系統 x、y、z 的三個座標值，因為 posp(u,v) 是極座標必須分別乘上 x(u,v)、y(u,v)、z(u,v) 以進行轉換如右：splot posp(u,v)*x(u,v), posp(u,v)*y(u,v), posp(u,v)*z(u,v)，輸入此行指令即能得到圖七的 $$s$$ 軌域圖形。如果指令一行寫不完時，可以倒斜線 $$(\backslash)$$ 相連，圖六的最後一行，代表除了描繪第一點，接續描繪第二點其座標為 (negp(u,v)*x(u,v), negp(u,v)*y(u,v), negp(u,v)*z(u,v))，而且使用不同的顏色，當然在 $$s$$ 軌域中看不出任何差異，但對於角動量量子數大於 $$1$$ 的軌域，就會有明顯的效果。

圖七、利用 gnuplot 畫出 $$s$$ 軌域的 3D 圖形（作者製作）

三、利用 gnuplot 繪製似氫原子的各種軌域

有了畫 $$s$$ 軌域的經驗，要畫其他軌域便不是難事，首先試試 $$2p$$ 的 $$3$$ 個軌域，我們在圖六中加入下列 $$4$$ 行 gnuplot 能認識的波函數如下（可直接由表一的正規化係數乘以角度函數得出）：

### 設定 p orbitals 有關角度部分的函數 ###

pz(u,v)=sqrt(3.0/4.0/pi)*cos(u)                         # m = 0
px(u,v)=sqrt(3.0/4.0/pi)*sin(u)*cos(v)          # m =-1
py(u,v)=sqrt(3.0/4.0/pi)*sin(u)*sin(v)              # m = 1

並將圖六中 orbp(u,v)= s(u,v) 這一行，換成 orbp(u,v)= pz(u,v)，即可繪出 $$p_z$$ 軌域的圖形如圖八。若要繪製 $$p_x$$ 或 $$p_y$$ 軌域，只要將圖六的 orbp(u,v) 換成等於 px(u,v) 或 px(u,v) 即可，讀者可自行嘗試。

眼尖的讀者會發現 x、y、z 軸的座標已重新設定過，圖形的標題也變成：$$p_z$$ orbital，至於如何更動，讀者可從圖六中稍作修改便能完成任務。

另外，圖形中多了 $$3$$ 個箭頭的座標指示，其實我們可以在圖形的任何地方繪製箭頭，只要給定箭頭的標號、開始及結束的座標即可，加上文字更簡單，只要給座標的位置及文字即可，以下即為再加入圖六中的片斷 (box=0.9)：

set arrow 1 from -box, 0, 0 to box, 0, 0
set arrow 2 from 0, -box, 0 to 0, box, 0
set arrow 3 from 0, 0, -box to 0, 0, box-0.2
set label ‘x’ at box,-.25*box,0
set label ‘y’ at 0,1.1*box,-.15*box
set label ‘z’ at 0.10*box,0,box-0.2

圖八、利用 gnuplot 畫出 $$p_z$$ 軌域的 3D 圖形（作者製作）

既然能繪製 $$p$$ 軌域，要繪出 $$d$$、$$f$$ 軌域自然不是難事，以下圖九是利用表一組合出 $$3d_{z^2}$$ 軌域的波函數所畫出的圖形，指令應怎麼寫？就留給讀者自行練習。

圖九、利用 gnuplot 畫出 $$d_{z^2}$$ 軌域的 3D 圖形（作者製作）

四、結論

Gnuplot 軟體具有強大的繪圖功能，除了免費、簡單易學及友善的溝通界面以外，還是一種跨平台的軟體，能在不同的電腦系統中使用，包括 window、osX及linux，且能產生各種通用的圖檔，例如 png、svg、ps、hpgl, … 等，， 提供文書處理、簡報、試算表… 等軟體匯入。因此深受各界的喜愛，廣泛的應用在各類學術領域上。另外，它適合各種不同程度的使用者，不懂程式語言的人，可用簡單互動式的輸入，將指令一步一步的執行。熟悉程式語言者則能在彈指間，快速精簡的完成各式各樣的繪圖。基於上述優點，使它具有高回報率的學習價值，因此強力推薦給所有理工科系的同學做為統計、繪圖之用。

本文即為 gnuplot 在化學科目上的一種應用，協助初學原子軌域的學子，能親自依據似氫原子的波函數輕易的繪製各種軌域的形狀，使軌域、波函數由抽象的數學式子，變成具體、可變大小、可由不同視角觀察的 3D 圖形。文中雖只畫出一部分的圖形，讀者若能從閱讀中邊做邊學，畫出其餘的軌域，甚至繪製 $$4f$$ 軌域，相信這個過程，會帶給你很高的成就感及投資報酬率。

參考文獻

1. Saputra, A.; Canaval, L. R.; Sunyono; Fadiawati, N.; Diawati, C.;Setyorini, M.; Kadaritna, N.; Budi Kadaryanto, K. Visualizing Three-Dimensional Hybrid Atomic Orbitals Using Winplot: An Applicationfor Student Self Instruction. Chem. Educ. 2015, 92, 1557−1558.
2. Chung, W. C. Three−Dimensional Atomic Orbital Plots in the Classroom Using Winplot. Chem. Educ. 2013, 90, 1090−1092.
3. Moore, B. G. Orbital Plots Using Gnuplot. Chem. Educ. 2000, 77, 785-789.
4. http://gnuplot.info/
5. http://user.frdm.info/ckhung/b/ma/gnuplot.php
6. http://www.ss.scphys.kyoto-u.ac.jp/person/yonezawa/contents/program/gnuplot/in dex.html
7. http://www.phy.fju.edu.tw/files/archive/876_ab57aed9.pdf
8. N. Levine(2008), Physical Chemistry (6th ed.). p637~647, McGRAW-HILL Book Company.
9. http://wps.prenhall.com/wps/media/objects/340/348272/Instructor_Resources/Chapter_02/Text_Images/FG02_14.JPG