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

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

學習原子軌域時，經常由於太過抽象，常使學習者敬而遠之，因此一般教科書均會各種圖形加以輔助說明，例如原子軌域的3D形狀、徑向分佈函數(radial distribution function)對原子半徑的變化圖及波函數的電子出現機率之等高線圖(contour map)等，但是學子仍然無法理解，這些圖形是如何繪製出來的？倘若有一個適宜的繪圖軟圖工具，讓學子親自繪製及感受各類軌域間的特徵及變化，相信對於學習效果會有事半功倍的效果。

在網路上頗為盛行的免費繪圖軟圖：gnuplot，由於具有可跨平台、輕薄短小、易於操作及功能強大等優點，已廣泛的被使用在各個學科領域中。筆者曾撰寫二篇其應用在量子化學上的文章：《利用Gnuplot軟體繪製似氫原子的軌域形狀》、《利用Gnuplot軟體繪製似氫原子的徑向分佈函數》，目前均放在高瞻平台上以供參考，前者主要是介紹如何利用 gnuplot 軟體，繪製各類原子軌域的立體形狀，後者則是在平面上繪製電子出現機率與原子半徑的關係圖，此二者均為探討原子軌域的一種觀點，前者是觀看軌域外部的形狀，內部的結構則無法判斷，後者則是由圖形中觀想軌域內部的分佈情形，外部的具體樣貎則需要想像力。

本文則應用 gnuplot 的各項指令，繪製結合前二者觀點的軌域等高線圖。另外，讀者在閱讀本文時，可一邊按照文中的指示操作電腦，相信從做中學，會有更深刻的領悟及更豐富的學習效果。

一、原子軌域等高線圖之說明

首先來了解何謂等高線圖？現以圖一 gnuplot 所繪出的 $$3d_{xy}$$ 軌域的圖形及其等高線圖為例加以說明。第一步必須將 $$(x,y)$$ 點上，$$3d_{xy}$$ 軌域電子出現的機密度計算出來，即等於其波函數的平方 $$(|\varphi_{3d_{xy}}|^2)$$，然後以其值為 $$z$$ 座標，在特定區間內將各個 $$(x,y,z)$$ 點依序描繪，可得圖一的上半部圖形，和平常 $$3d_{xy}$$ 軌域形狀的認知相同，有二個節面，即電子不會出現的區域。

聰慧的讀者可能會發現，波函數不是也需要 $$x$$、$$y$$ 和 $$z$$ 等 $$3$$ 個座標變數嗎？現在將 $$z$$ 用做電子出現的機率密度，那原先波函數的 $$z$$ 值應如何處理呢？此點容後說明，現在將焦點轉移到下半部的等高線圖。

圖一、利用 gnuplot 軟體繪製 $$3d_{xy}$$ 軌域的電子出現機率密度之分佈情形及其等高線圖（作者繪製）

試著想像一下，拿一塊平板與 $$xy$$ 平面平行，在 $$z$$ 軸等於 $$0.0001$$ 的高度橫切過去，與圖形相交的點將其投影在 $$xy$$ 平面上，應可得到如圖一下半部中，$$4$$ 個橢圓形的藍色圈圈，若將平板等距 $$(0.0001)$$ 上升，再重覆上述動作 $$4$$ 次，則可完成圖一下半部圖形的全貎，此即為等高線圖。

若經過練習，單從下半部的等高線圖來判讀，便能反推上半部的 3D 形狀，例如由等高線圖可知其為有 $$4$$ 個山峰的形狀，靠近原點附近的 $$4$$ 個點，電子出現的機密最大。圖形外部的等高線之間較疏，代表其山坡斜度較平緩，而內側等高線間較密，代表其山坡較陡，另外，在  $$x$$ 軸及 $$y$$ 軸是電子不會出現的地方，3D 圖形中即為節面。

圖一是以電子出現機率密度（波函數的平方）做圖，所以 $$z$$ 值没有負值，如果使用波函數做圖，則 $$z$$ 值有正亦有負，了解等高線圖以後，接著來試試如何在 gnuplot 軟體中，實際繪製出等高線圖。

二、利用 gnuplot 繪製 $$3p_x$$ 軌域及其等高線圖

由薛丁格方程式 (Schrödinger equation) 解出的似氫原子 (hydrohen like atom) $$2p_x$$、$$3p_x$$、$$3d_{xy}$$ 的波函數分別如下：

$$\varphi_{2p_x}=\displaystyle\frac{1}{4\sqrt{2\pi}}(\frac{Z}{a_0})^{\frac{5}{2}}re^{-Zr/2a_0}\sin\theta\cos\phi$$          （式-1）

$$\varphi_{3p_x}=\displaystyle\frac{\sqrt{2}}{81\sqrt{\pi}}(\frac{Z}{a_0})^{\frac{5}{2}}(6-\frac{Zr}{a_0})re^{-Zr/3a_0}\sin\theta\cos\phi$$          （式-2）

$$\varphi_{3d_{xy}}=\displaystyle\frac{1}{81\sqrt{2\pi}}(\frac{Z}{a_0})^{\frac{7}{2}}r^2e^{-Zr/3a_0}\sin^2\theta\sin 2\phi$$          （式-3）

上式中 $$a_0$$ 為波耳半徑等於 $$0.53$$ Å，$$Z$$ 為原子核的正電數，對於氫原子而言 $$Z=1$$，若作圖時原子半徑再以 $$r/a_0$$ 做單位，則上列式子的 $$Z$$ 及 $$a_0$$ 均可省略。

為了要繪製上列波函數的圖形，首先必須將極座標轉換成我們熟悉的直角座標。有關直角座標和極座標的轉換方式，如圖二所示，$$(x,y,z)$$ 轉成 $$(r,\theta,\phi)$$ 可透過下列公式完成：

$$x=r\sin\theta\cos\phi$$

$$y=r\sin\theta\sin\phi$$

$$z=r\cos\theta$$

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

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

接下來，代表 $$3p_x$$ 的（式-2）在軟體中如何繪製波函數 $$\varphi_{3p_x}(r,\theta,\phi)$$ 隨 $$(r,\theta,\phi)$$ 的分佈圖？

這裏出現 $$4$$ 個變數，但一般 3D 繪圖只能有 $$3$$ 個變數，因此只能選擇一個特殊的視角，本文以 $$xy$$ 平面上來看波函數的變化，即 $$z=0$$，由圖二可知，此時 $$\theta=\pi/2$$ 為一固定值，此時改變 $$r$$ 值，從 $$0$$ 到 $$20$$，改變 $$\phi$$ 值，從 $$0$$ 到 $$2\pi$$，並將每一個 $$(r,\frac{\pi}{2},\phi)$$ 代入波函數出求出其值，然後以 $$x=r\cdot\sin\frac{\pi}{2}\cos\phi$$、$$y=r\cdot\sin\frac{\pi}{2}\sin\phi$$、$$z=\varphi_{3p_x}(r,\frac{\pi}{2},\phi)$$，將此 $$3$$ 個變數當成 3D 空間的一點，隨著 $$r$$、$$\phi$$ 值的變化再一點一點描繪上去即可完成任務。

有了上述繪圖邏輯以後，即可將其想法移到 gnuplot 軟體中，開始撰寫指令集的文字檔如圖三：大部分的指令已由前言中所提的二篇文章中說明過，讀者可自行前往高瞻平台參考，但為了方便了解及查詢，仍在圖三中加註說明，其中只要有「#」符號出現，該行從這個位置以後的文字均為說明或註解，並不會真正執行。

圖三、gnuplot 指令集，藉以繪出氫原子 $$2p_x$$、$$3p_x$$、$$3d_{xy}$$ 軌域之波函數在 $$xy$$ 平面上的分佈圖（作者繪製）

在圖三中的第 $$6$$ 行：set urange[0:20.] ; set vrange[0:2.0*pi]，即為設定 gnuplot 的內定變數 $$u$$ 及 $$v$$ 的範圍，相當於文中所提極座標的 $$r$$ 及 $$\phi$$，分別從 $$0$$ 到 $$20$$ 及 $$0$$ 到 $$2\pi$$，程式會依所設定的解析度，一點一點代入圖三的最後一行：splot u*cos(v) ,u*sin(v) ,(p3xorbital(u,v)) title “”，將極座標轉換成直角座標，並在 $$xy$$ 平面上描繪 $$3p_x$$ 波函數的分佈圖，詳如下篇圖四。另外圖三倒數第 3-5 行波函數的設定，即依（式-1）、（式-2）及（式-3）撰寫成 p2xorbital(x,y)、p3xorbital(x,y) 及 dxyorbital(x,y)。

連結：利用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.