利用Gnuplot軟體繪製似氫原子軌域的等高線圖（下）

利用Gnuplot軟體繪製似氫原子軌域的等高線圖（下）(Using Gnuplot plot the contour map of hydrogen-like orbital (II))
國立臺灣師範大學化學系兼任教授 邱智宏

連結：利用Gnuplot軟體繪製似氫原子軌域的等高線圖（上）

圖四、利用 gnuplot 軟體繪製 $$3p_x$$ 軌域波函數的形狀及其等高線圖（作者繪製）

圖四中 $$z$$ 軸即為 $$3p_x$$ 波函數的值，由圖中可看出其值有正亦有負，正值代表在 $$xy$$ 平面上方有凸起的部分，有一個山峰及一個較平緩的小山丘，凹下去的部分為負值，有一個深的山谷及一平緩的小山溝。一般教科書常用等高線圖的方式，用平面圖形表達三度空間的訊息，在圖三中的第 15 到 17 行即為在 gnuplot 中繪製等高線圖的指令：

set contour
set cntrparam levels 9
set cntrparam levels incremental -2.e-02,5.0e-03

第一行設定要畫等高線圖，第二行要畫 $$9$$ 條等高線圖，第 $$3$$ 行代表 $$z$$ 軸的大小從 $$-0.02$$ 開始畫第 $$1$$ 條，每增加 $$0.005$$ 增加一條，可以想見的在第 $$5$$ 條的等高線是畫在 $$z=0$$ 即波函數為 $$0$$ 的地方，從圖四中的下半部可看出圍住中間二團橢圓圏圏的淡藍色圓圏和 $$y$$ 軸，即為節面的所在。圖四中的座標軸格式及最右邊的圖例，在 gnuplot 中是用下列指令來控制：

set tics font “Times,11” ; set key at screen 0.95,0.84

tics 為座標軸上數字的設定，字體設為 Times，大小為 11 號字。Key 則為圖例的設定，程式有內定的值包括大小、顏色等，我們可以自行變換，在這裏只設定其擺放在螢幕的位置 (at screen 0.95,0.84)。如果我們想繪製 $$2p_x$$ 軌域，其電子出現機率密度在 $$xy$$ 平面上的分佈圖及等高線圖，僅需將圖三的後一行修改為下式即可：

splot u*cos(v),u*sin(v),(p2xorbital(u,v)**2) title “”

由於 $$2p_x$$ 及 $$3d_{xy}$$ 的波函數在圖三中均已設定過，所以直接代入上式即可，其中 p2xorbital(u,v)**2 表示為 $$2p_x$$ 波函數的平方，即為電子出現的機率密度。另外圖三的第 17 行，要修改如下：

set cntrparam levels incremental 1.0e-4,1.0e-03

重新執行程式後所得的結果如圖五，由於電子出現機率密度没有負值，所以圖五只有山峰没有山谷，而機率密度為 $$0$$ 的地方，在圖五中並没有顯現出來。

圖五、利用 gnuplot 軟體繪製 $$2p_x$$ 軌域電子出現機率密度分佈及其等高線圖（作者繪製）

如果想要繪製圖一 $$3d_{xy}$$ 軌域的電子出現機率密度之分佈情形及其等高線圖，只要將圖三的最後一行改成：splot u*cos(v),u*sin(v),(dxyorbital(u,v)**2) title “”，當然還要修改一些細節，例如標題、等高線的範圍等，讀者可自行練習。

三、等高線圖及極大值

目前為止 gnuplot 軟體似乎很神奇，只要給它一個函數，及適當的參數及座標軸範圍，便能迅速的畫出精緻實用的圖形，然而使用者仍有霧裏看花的感覺，這個軟體究竟是如何運作的？

事實上它是一種專注於繪圖的軟體工具，適用於各學科，其核心機制就是輸入一堆數據，畫出最佳圖形，因此只要依照它的邏輯，給它適當的 $$u$$、$$v$$ 參數範圍，程式就依解析度的設定（以3D繪圖為例 isosample 50,50），自動將 $$u$$ 及 $$v$$ 的區間各切割為 $$50$$ 點，然後將一點一點的參數值代入函數中，例如 $$f(u,v)$$，算出其值，然後按照各點的座標 $$(u,v,f(u,v))$$ 將該點描繪出來，最後輸出圖形。

依照這個邏輯，圖形中的任一點的函數值，程式已經運算處理完畢，但是我們不知道程式是如何儲存？存在何處？還好我們依然可以依照其邏輯，在 gnuplot 中寫一段指令集，取出我們想要特殊的位置座標，例如想要截取圖四 $$3p_x$$ 函數的極大值及極小值？可以在圖三最後一行畫出圖形的指令之前，加入下列指令集：

wmax=-1.0 ; wmin=1.0

do for [i=0:50] {

do for [j=0:50] {

wx=(20./50.)*i

wy=(2.*pi/50.)*j

if (p3xorbital(wx,wy)>=wmax) {wmax= p3xorbital(wx,wy);xmax=wx;ymax=wy}

if (p3xorbital(wx,wy)<=wmin) {wmin=p3xorbital(wx,wy);xmin=wx;ymin=wy}

}

}

上列指令集的大意為：先設定 $$2$$ 個變數 wmax、wmin 分別存放波函數的極大值及極小值，先預設各為 $$-1.0$$ 及 $$1.0$$，顯然此值並不合理，經過後面程式的檢查，一定會被取代掉。

然後將 $$x$$、$$y$$ 軸的座標指定到 wx、wy，將 $$x$$ 軸的區間 $$20$$，分成 $$51$$ 點，每點的間隔 $$20/50$$，將 $$y$$ 軸的區間 $$2\pi/50$$，分成 $$50$$ 點，每點的間隔 $$2\pi/50$$，接下來由 i、j 二個迴圏，來檢查波函數的極大值及極小值，i 值從 $$0$$ 開始到 $$50$$ 結束，每次 wx 改變為 wx=(20/50)*i，wy 的改變由 j 控制：wy=(2π/50)*j。

當固定一個 i 值時，分別代入不同的 $$51$$ 個 j 值，每次均以以 if 指令檢查 p3xorbital(wx,wy)>=wmax 是否為真？如果為真，則 wmax 就以 p3xorbital(wx,wy) 取代，順便將此時的座標 wx、wy 以 xmax、ymax 記錄下來，如果 if 指令邏輯檢查後若為偽，則波函數的極大值不做取代，直至 i、j 二個迴圏均完成。

另外，極小值也是以相同的方法檢查，執行完畢時，就會二組數據，分別儲存波函數的極大值及極小值的座標。我們也可以下列指令，將其座標值以藍色字體，在螢幕的適當位置列印出來：

set label sprintf(“wavefunction max = %7.4f, %7.4f, %7.4f”,xmax*cos(ymax), xmax*sin(ymax), \ wmax) at screen 0.5,0.4 textcolor “blue”

set label sprintf(“wavefunction min = %6.4f, %6.4f, %6.4f”,xmin*cos(ymin), xmin*sin(ymin), \ wmin) at screen 0.5,0.35 textcolor “blue”

有了上列觀念以後，便可以進一步了解，等高線圖是以何種邏輯繪製出來的。

想要在二度空間表達三度空間的訊息，首先必須了解三度空間的所有座標值，接下來選擇 $$xy$$ 平面，此平面的每個 $$(x,y)$$ 座標均有一個 $$z$$ 值，此時我們可以選擇 $$z$$ 值為特定值的所有點，將其連結成一條線即為等高線，此說法和先前的敘述相同，此等高線相當於 $$z$$ 等於某特定值時和 $$xy$$ 平行的平面與波函數的相交點所連結的線，再將其投影到 $$xy$$ 平面上，取幾組 $$z$$ 值，就有幾條等高線。

依照這個邏輯，我們也可以在 gnuplot 內寫一個小程式，單獨將等高線圖畫出來，幸好軟體內有現成指令可以輕易完作任務，以圖四、圖五所列，如欲單獨呈現等高線圖可在圖三的指令集中，在 set contour 之後加上一行指令 set view map 及 unset surface 即可得到如圖六的等高線圖。

圖六可以更清楚的看出等高線圖的面貎，線間愈密集的地方山坡愈徒，而且可以看到二個節面，一個為 $$r=6$$ 的棕色圏圏，一個為 $$y$$ 軸。顯然完成圖六的圖形尚經過一些修飾的指令，例如相同的等高線用一樣的顏色，波函數的值為正時，用實線繪製，為負時用虛線等。

另外，$$x$$ 軸、$$y$$ 軸的比例，圖形的標題亦須經過調整，由於無關宏旨，在這裏就不加說明，讀者可以自行嘗試。

圖六、利用 gnuplot 單獨繪製 $$3p_x$$ 的等高線圖（作者繪製）

四、結論

上二篇相關的文章如前言所述，曾介紹 gnuplot 軟體的基本指令操作、及類似程式語言指令集的撰寫，並應用在繪製似氫原子各類軌域的立體圖形、機率密度函數及徑向分佈函數對原子半徑的分佈圖等。本文則偏向介紹軟體運作的邏輯，如何處理數據及繪圖，使讀者知其然以外，更能知其所以然，最重要的還是希望能將文中的所有圖形，均能親自繪製出來，經過實做的磨練，會有如人飲水，了然於胸的感覺，甚至激發出不同的想法和應用。

本文除了學習 gnuplot 軟體的應用以外，並利用它來繪製似氫原子的波函數在固定平面上，其波函數或波函數平方的分佈圖，使讀者除了明白軌域的形狀以外，更能解析其內部的結構，也能了解繪製等高線圖的原理，並實際在軟體中繪製出來，甚至可以撰寫小程式，將波函數的極大值或極小值求出來。文中僅以 $$2p_x$$、$$3p_x$$、$$3d_{xy}$$ 的波函數為例，事實上只要知道波函數，不管主量子數、角動量量子數如何改變，其繪製的方法還是相同，利用 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. Spherical to Cartesian Coordinates Calculator – Learning about Electronics. http://www.learningaboutelectronics.com/Articles/Spherical-to-cartesian-rectangular-coordinate-converter-calculator.php
6. N. Levine (2008), Physical Chemistry (6th ed.). McGRAW-HILL Book Company. p637~647.
7. 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.
8. Chung, W. C. (2013). Three-dimensional atomic orbital plots in the classroom using Winplot. Journal of Chemical Education, 90(8), 1090-1092.
9. Moore, B. G. (2000). Orbital Plots Using Gnuplot. Journal of Chemical Education, 77, 785-789.