R语言绘制氢原子原子轨道波函数

Posted on | In ,

原文地址 zhuanlan.zhihu.com

解氢原子薛定谔方程

薛定谔方程
$$
\left(-\frac{h^{2}}{8 \pi^{2} m_{e}} \nabla^{2}-\frac{Z e^{2}}{4 \pi \varepsilon_{0} r}\right) \psi=E \psi
$$
坐标轴转换
$$
\begin{cases}
x=r \sin \theta \cos \varphi \
y=r \sin \theta \sin \varphi \
z=r \cos \theta
\end{cases}
$$

$$
\begin{cases}
r^{2}=x^{2}+y^{2}+z^{2} \
\cos \theta=\frac{z}{\left(x^{2}+y^{2}+z^{2}\right)^{\frac{1}{2}}} \
\tan \varphi=\frac{y}{x}
\end{cases}
$$

拉普拉斯方程球坐标表示
$$
\nabla^{2}=\frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^{2} \frac{\partial}{\partial r}\right)+\frac{1}{r^{2} \sin \theta} \frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial}{\partial \theta}\right)+\frac{1}{r^{2} \sin ^{2} \theta} \frac{\partial^{2}}{\partial \varphi^{2}}
$$
球坐标下的薛定谔方程
$$
\frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^{2} \frac{\partial \psi}{\partial r}\right)+\frac{1}{r^{2} \sin \theta} \frac{\partial}{\partial \theta}\left(\sin \theta \frac{\partial \psi}{\partial \theta}\right)+\frac{1}{r^{2} \sin ^{2} \theta} \frac{\partial^{2} \psi}{\partial \varphi^{2}}+\frac{8 \pi^{2} m}{h^{2}}\left(E+\frac{Z e^{2}}{4 \pi \varepsilon_{0} r}\right) \psi=0
$$
通过解薛定谔方程（解法略过）可得氢原子的定态波函数：
$$
\Psi(r, \theta, \phi)=c\left[\sum_{I=1}^{n-1} c_{i}\left(\frac{Z r}{a_{0}}\right)^{l+i-1} e^{-\frac{Z_{r}}{n a_{0}}}\right] P_{l}^{|m|}(\cos \theta) e^{i m \phi}
$$

$n$为主量子数，取值为$1,2,3,…$

$l$为角量子数，取值为$0,1,2,…,n-1$

$m$为磁量子数，取值为$-l,-l+1,…,0,…,l-1,l$

R 语言绘制波函数三维散点图

R 语言中可以绘制三维图像的包很多：plot3D,plotly,rgl,threejs等

目前小编还没有找到 R 语言直接绘制波函数原子轨道轮廓图的方法。这里我只展示绘制散点图来展示原子轨道的大致轮廓。

我们来绘制一下氢原子的$3d_{z^2}$轨道，波函数省去常数为：
$$
\psi = r^2 e^{-r/3} (3 \cos^2 \theta - 1)
$$

我的绘制思路就是在一个三维正方体域内计算波函数值，再求解得到的波函数值平方的平均值。然后绘制平均值上下范围内的散点即可。这里我用rgl包的plot3D函数做演示。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
library(rgl)
​
# 定义波函数
wavefunction <- function(r,theta,phi) r^2*exp(-r/3)*(3*cos(theta)^2-1)
​
# 设定超参，m不应小于20否则绘制图像不完整
num <- 100  # 反映散点密度
m <- 20 # 范围
​
# 初始化一个正方体域
xt <- rep(seq(-m,m,length = num),each = num^2)
yt <- rep(rep(seq(-m,m,length = num),each = num),num)
zt <- rep(seq(-m,m,length = num),num^2)
​
# 转为球坐标参数
r <- sqrt(xt^2+yt^2+zt^2)
theta <- acos(zt/r)
phi <- atan(yt/xt)
​
# 计算波函数
wv <- wavefunction(r,theta,phi)
​
mn <- mean(wv^2,na.rm = T)  # 波函数平方的平均值
drw <- !is.na(wv)&wv^2<mn&wv^2>0.8*mn   # 找出波函数平方平均值0.8~1倍的点集
​
cl <- ifelse(wv>0,"red","blue") # 波函数为正的用红色标注，负的用蓝色
cl <- cl[drw]
​
# 找出满足波函数的x,y,z坐标
x <- xt[drw]
y <- yt[drw]
z <- zt[drw]
​
# 绘制三维散点图
plot3d(x,y,z,size = 1,type = "s",col = cl,alpha = 0.1)

此外我们可以试试其他包函数绘制

1
2
3
library(threejs)
scatterplot3js(x, y, z, color =  cl, size = 0.1)

1
2
3
4
# 还可以尝试使用其他相关的3D绘图包进行绘制
library(plotly)
library(plot3D)

示例展示