您的位置: 主页 > 改变头脑 >利用Gnuplot软体绘製似氢原子的径向分布函数(下) >

利用Gnuplot软体绘製似氢原子的径向分布函数(下)


2020-06-19


连结:利用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$$

利用Gnuplot软体绘製似氢原子的径向分布函数(下)

图五、极座标中单位体积 $$(dV)$$ 的表示法示意图(图片来源:参考资料6)

由于 $$1s$$、$$2s$$ 和 $$3s$$ 轨域和 $$\theta$$、$$\phi$$ 无关,因此可单看径度变化的影响如图六所示,探索在从 $$r$$ 到 $$r+dr$$ 间圆壳 (spherical shell) 内找到电子的机率有多少?

利用Gnuplot软体绘製似氢原子的径向分布函数(下)

图六、$$1s$$ 轨域径度从 $$r$$ 到 $$r+dr$$ 间圆壳体积的示意图(图片来源

若将 $$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软体绘製似氢原子的径向分布函数(下)

图七、利用 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软体绘製似氢原子的径向分布函数(下)

图八、利用 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软体绘製似氢原子的径向分布函数(下)

图九、利用 gnuplot 的指令,计算电子在 $$2s$$ 轨域第一个小山丘出现的机率(图片来源:作者製作)

图九中可看出透过 gnuplot 使用辛普森数值法所求出,$$r$$ 从 $$0$$ 到 $$2a_0$$ 区间的积分值为 $$0.0527$$,和精确解的答案几乎一样。当然我们也可以改变不同的区间,改变不同的函数来做运算,但利用 gnuplot 同样的可以轻易的求出準确的解答。

四、结论

上一篇相关的文章如前言所述,曾利用 gnuplot 软体,绘製似氢原子各类轨域的立体图形,本文则介绍其在平面绘图上的应用,一面複习使用过的指令,一面增添常用及新的指令,尤其是类似程式语言的部分例如 do for、if 等,期能对此软体的应用能力,有大幅提升的功效。对于熟悉程式语言的读者,文中的示範说明,虽然是多此一举,但说不定能激发其创意,并应用在其他不同的领域上也未可知;对于不熟悉程式语言的读者,也可以按图索骥,边读、边做、边学,一定会有意想不到的功效。

本文除了学习 gnuplot 软体的应用以外,并利用它来绘製似氢原子的波函数、机率密度函数及径向分布函数对原子半径的分布图,让读者可以亲自将教科书中抽象的函数,转化为具体的图形,更可以藉由座标轴範围的改变,局部放大或全面观察各种图形中关键性的变化。另外,也藉机釐清:为何电子密度最高的原子核,却是最找不到电子的地方?更能利用该软体以数值解的方式将任何轨域在不同区间中,电子出现的机率总和给核算出来。


参考文献



上一篇:
下一篇: