用Python画氢原子的电子云

最新推荐文章于 2023-11-09 20:39:26 发布

微小冷

最新推荐文章于 2023-11-09 20:39:26 发布

阅读量2.6k

点赞数 3

分类专栏： # Python可视化文章标签： python 电子云薛定谔方程量子力学原子物理

本文链接：https://blog.csdn.net/m0_37816922/article/details/128788923

版权

Python可视化专栏收录该内容

72 篇文章

订阅专栏

本文介绍了如何使用Python来实现薛定谔方程，特别是在球坐标系下解决电子云的分布问题。通过结合球谐函数和拉盖尔多项式，计算不同量子数状态下的电子云，并利用matplotlib进行三维散点图的绘制，展示电子在原子中的概率分布。文章还展示了从简单的1s状态到更复杂的2p1和3d2电子云的生成过程。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

文章目录

薛定谔方程

请不要心脏骤停，这里只是复制粘贴一下，本文主要是画图，和解方程没有任何关系

$-\frac{\hbar^2}{2\mu}\nabla^2\psi+V(r)\psi=E\psi,\quad V(r)=-\frac{e^2}{r\pi\epsilon_0r}$

使用球坐标，将Laplace算子展开

$-\frac{\hbar^2}{2\mu r^2}\{\frac{\partial}{\partial r}(r^2\frac{\partial}{\partial r})+\frac{1}{\sin^2\theta}[\sin\theta\frac{\partial}{\partial\theta}(\sin\theta\frac{\partial}{\partial\theta})+\frac{\partial^2}{\partial\phi^2}]\}\psi-\frac{e^2}{r\pi\epsilon_0r}=E\psi$

其解为径向函数和球谐函数的乘积

$\psi(r,\theta,\phi)=R_{nl}(r)Y_{lm}(\theta,\phi)$

其中 $Y_{lm}$ 是球谐函数，在scipy中有定义，可直接调用sph_harm， $R_{nl}$ 可以继续展开为

$R_{nl}=\sqrt{(\frac{2}{na})^3\frac{(n-l-1)!}{2n[(n+l)!]^3}}e^{-\frac{r}{na}}(\frac{2r}{na})^lL^{2l+1}_{n-l-1}(\frac{2r}{na})$

其中 $a$ 是玻尔半径0.53 $\text{\r{A}}$ ，在写代码是带入0.53，表示整个模型以埃米为单位。

Python实现

$L^{2l+1}_{n-l-1}$ 是拉盖尔多项式，在scipy中同样有定义，可直接调用assoc_laguerre

import numpy as np
import scipy.special as ss
def hydrogen(n,l,m,r,th,phi):
    na = 0.53*n
    tmp = (2/na)**3
    tmp *= ss.factorial(n-l-1)
    tmp /= 2*n*ss.factorial(n+l)**3
    tmp2 = np.exp(-r/na) * (2*r/na)**l
    laguerre = ss.assoc_laguerre(2*r/na, n-l-1, 2*l+1)
    R_nl = np.sqrt(tmp)*tmp2*laguerre
    return R_nl*ss.sph_harm(m,l,th,phi)

接下来生成一个主量子数为1，角量子数、磁量子数皆为0的电子云，由于在随机生成坐标点时，越是靠近原点，坐标点越密集，所以通过距离调整概率

from itertools import product
rs = np.arange(50)/50
ths = np.arange(0, np.pi*2, 0.1)
phis = np.arange(0, np.pi, 0.1)

pts = []
for axis in product(rs, ths,phis):
    p = hydrogen(1,0,0,*axis)**2
    if p*axis[0]**2 > np.random.rand()/2:
        pts.append(axis)

绘图

import matplotlib.pyplot as plt
def pts2xyz(pts):
    r, th, phi = zip(*pts)
    x = r*np.sin(phi)*np.cos(th)
    y = r*np.sin(phi)*np.sin(th)
    z = r*np.cos(phi)
    return x,y,z

x,y,z = pts2xyz(pts)
ax = plt.subplot(projection='3d')
ax.scatter(x, y, z, marker='.')
plt.show()

效果如下，至少可以看出，随着半径的增加，电子出现的概率越来越小

在这里插入图片描述

更复杂的电子云

接下来看一下2p1的情况，考虑到能级越高越不稳定，即出现的概率越低，所以对随机数的量级进行调整。又考虑到等间隔生成坐标点，过于刻意，所以随机生成坐标点

from numpy.random import rand

def getPts(n, l, m, N, rMax, randMax):
    r = rand(N)*rMax
    th = rand(N)*np.pi*2
    phi = rand(N)*np.pi
    pts = []
    for axis in zip(r,th,phi):
        p = hydrogen(n, l, m,*axis)**2
        if p*axis[0]**2 > np.random.rand()*randMax:
            pts.append(axis)
    print(len(pts))
    return pts

pts = getPts(2, 1, 1, 200000, 5, 1/200)
x,y,z = pts2xyz(pts)
ax = plt.subplot(projection='3d')
ax.scatter(x, y, z, marker='.')
plt.show()

效果为

在这里插入图片描述

继续看更高的能级

pts = getPts(3, 2, 2, 200000, 15, 1/1e5)
x,y,z = pts2xyz(pts)
ax = plt.subplot(projection='3d')
ax.scatter(x, y, z, marker='.')
plt.show()

在这里插入图片描述

如果把主量子数为3，角量子数为2的所有电子云悉数画出，则效果如下

fig = plt.figure()
ax0 = fig.add_subplot(231,projection='3d')

ms = [0,-1,1,-2,2]
cs = ['r', 'g', 'b', 'c', 'y']
for i in range(5):
    ax = fig.add_subplot(2,3,i+2,projection='3d')
    pts = getPts(4, 2, ms[i], 50000, 25, 1/1e7)
    x,y,z = pts2xyz(pts)
    ax0.scatter(x, y, z, marker='.', c=cs[i])
    ax.scatter(x, y, z, marker='.', c=cs[i])
    plt.xticks([])
    plt.yticks([])

plt.show()