前几天看b站up主solara570的视频，并想尝试其思路实现氢原子的电子云可视化，具体原理不是很复杂，也就作为一个尝试和学习思考的过程，记录下来。
下面附上代码，具体公式的分析，有时间再写吧。
简单写一下思路：
其实这个问题分为以下几个部分
1.氢原子的薛定谔方程是有解析解的，可以计算其s，p，d，f原子轨道
2.用py的numpy可以很容易的写出这些公式
3.可视化的方法，matplotlib是可以的，我就举了两种可行的方法，但是不是很优美的方法，其实最优美的方法应该是blender渲染，但我不会。

先添加所需要的库，如numpy，scipy中的factorial（阶乘），genlaguerre（连带拉盖尔函数）和sph_harm（球谐函数）当然还有一些常见的绘图库，
这个skimage是后续要计算等势面需要的库

1
2
3
4
5
import numpy as np
from scipy.special import  factorial, genlaguerre, sph_harm
from skimage.measure import marching_cubes
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

确定玻尔半径，将氢原子的公式用代码表示出来，但我并不很熟练latex，所以我就先不放这个公式了，等我有时间我再放这个公式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
a0 = 1
def hydrogen_wave_function(n,l,m):
    #注意连带拉盖尔函数的输入参数顺序

    def R(r):
        factor = np.sqrt((2. / (n * a0))**3 * \
                     factorial(n - l - 1) / (2 * n * factorial(n + l))) 
        rho = 2 * r / (n * a0)  # 中间量

        return factor * (rho**l) * np.exp(-rho / 2) * \
           genlaguerre(n - l - 1, 2 * l + 1)(rho)

    
    def Y(theta,phi):
        return  sph_harm(m, l, phi, theta) #注意这里theta和phi是相反的变量
    return lambda r, theta, phi: R(r) * Y(theta, phi)

坐标系的变化-转换到球坐标系中

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#定义原子轨道
n,l,m = 2,1,0 #这是量子数，我选择了2，1，0 量子数，也就是p轨道，
psi = hydrogen_wave_function(n,l,m)
print(f"{psi = }")#psi是一个复数是正确的， 
#空间网格化 
limit = 10
n_points=50
vec = np.linspace(-limit,limit,n_points)
print(f"{len(vec)= }")
#生成三维网格，记录每个点坐标
X, Y, Z = np.meshgrid(vec, vec, vec)
#坐标变化到柱坐标系
R = np.sqrt(X**2 + Y**2 + Z**2) 
THETA = np.arccos(Z / R) 
PHI = np.arctan2(Y, X)
psi_values = psi(R, THETA, PHI) 
print(f"{psi_values.shape = }")

后续就是寻找等势面，计算出有相同概率的空间位置，

1
2
3
4
5
6
7
8
9
10
11
12
#等概率面寻找
prob_dens = np.abs(psi_values)**2 #概率密度函数
iso_value = 4e-4 
#_这个的意思是这个数据我们丢弃，我们并不关系这个数值
step = 2 * limit / (n_points - 1) 
verts, faces, _, _ = marching_cubes( 
    prob_dens, 
    level=iso_value, 
    spacing=(step, step, step), 
)
verts -= limit 
verts[:, [0, 1]] = verts[:, [1, 0]]

后面是可视化的部分，视频给出了三种方法，但我研究了其中两种

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
36
37
38
39
40
41
#这是一个画布
def new_fig_and_ax(plot_range=10): 
    fig = plt.figure() 
    ax = fig.add_subplot(1, 1, 1, projection="3d")
    # 设置坐标比例关系，让图像变得方正
 
    ax.set_box_aspect([1, 1, 1]) 
     
    # 设置合适的视图俯仰角和方位角，让坐标轴方向符合直观
 
    # x轴对着屏幕外方向，y轴朝右，z轴朝上
 
    ax.view_init(elev=30, azim=20)
    ax.set_xlabel("$x$") 
    ax.set_ylabel("$y$") 
    ax.set_zlabel("$z$")
    ax.set_xlim(-plot_range, plot_range) 
    ax.set_ylim(-plot_range, plot_range) 
    ax.set_zlim(-plot_range, plot_range) 
    ticks = np.linspace(-plot_range, plot_range, 5) 
    ax.set_xticks(ticks) 
    ax.set_yticks(ticks) 
    ax.set_zticks(ticks)
    return fig, ax 
fig, ax = new_fig_and_ax() 
 
iso_surface = ax.plot_trisurf(
    verts[:, 0], verts[:, 1], faces, verts[:, 2],
    lw = 0,
    cmap = "coolwarm_r") 

plt.show
#第二种方法，当然也可以用blender做渲染
fig, ax = new_fig_and_ax() 
mesh = Poly3DCollection(verts[faces], lw=0.1) # 
mesh.set_facecolor("yellow") # 黄色的表面
 
mesh.set_edgecolor("grey") #
 
ax.add_collection3d(mesh) #
plt.show()

亲测有效，可以拿去跑一跑，试一试。改变量子数可以观察出不同轨道的可视化形状

经过计算可以看出来有一些问题，p轨道和高中学的不一样，d轨道也和高中学的不一样，这是正常的，因为这是物理学家和化学家采用了不同的轨道表达

化学家在研究化学键和分子轨道的时候会有以下两个问题：
1.函数是复数域内的函数，并不直观；
2。对称性考虑上有一些轨道对称性直观不同（就是长得不一样），但其实本质相同（能级间并，而且如果没有磁效应，磁量子数不起作用，没有zeeman效应）
基于这个观点考量，以为轨道之间正交且线性无关，因此正交基可以选择其他的，化学家通过线性变换找到三组x，y，z对称的p轨道，来方便解释化学键的交叠。
这一点在d轨道上尤其明显（这个表示一部分是考量了对称性，另一部分其实是实数化后e指数上的数值来命名的），这就涉及好多推导。以后可能有兴趣，我自己进行推导，来给自己加深印象。
关于这个网站使用问题，我未来可能会把这个当作我的垃圾站，做一些随想和一些机器学习深度学习代码实现，仿真实现，然后留着给自己当备份了。
下一个要更新的我基本想好了。一个是k近邻算法的实现，另一个是循环神经网络预测股票。