氢原子电子3D散点分布图制作介绍，通过概率密度生成随机点代码实现介绍-CSDN博客

本文链接：https://blog.csdn.net/qq_37668266/article/details/136700858

1.通过概率密度生成随机点

由量子力学可知，波函数 $\psi$ 的模长平方 $|\psi|^2$ 就是概率密度，下面介绍如何用概率密度生成随机点，以Python为例。假设现在有概率密度（1维情况）：

$\rho=e^{-x^2}$

import numpy as np
import matplotlib.pyplot as plt

#定义函数
def f(x):
    y = np.exp(-x**2)
    return y

#随机点个数
n = 300

#生成均匀分布的随机点
rx = np.random.random(n)#随机点x坐标
rx = (rx-0.5)*10
ry = np.random.random(n)#随机点y坐标
ry = ry*1.5

#存储正确的随机点
correct_random_x = np.array([])
correct_random_y = np.array([])

#存储错误的随机点
error_random_x = np.array([])
error_random_y = np.array([])

#判断所有随即点位置
for i in range(300):
    #！！！！！！！关键代码！！！！！！！
    if f(rx[i]) >= ry[i]:
        correct_random_x = np.append(correct_random_x,rx[i])
        correct_random_y = np.append(correct_random_y,ry[i])
    else:
        error_random_x = np.append(error_random_x,rx[i])
        error_random_y = np.append(error_random_y,ry[i])

x = np.linspace(-5,5,200)
y = f(x)

#画图
plt.scatter(correct_random_x,correct_random_y,c = 'b')
plt.scatter(error_random_x,error_random_y,c = 'r')
plt.plot(x,y)
plt.show()

由上述代码可以生成：蓝色点就是所要的随机点，即变量名为correct_random_x和correct_random_y所储存的值就是所要的随机点的x和y坐标。

2.绘制氢原子图

氢原子电子分布概率密度是3维，下面就以氢原子4p为例绘制:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import scipy.special as sp

#生成随机点个数
n = 10000

rx = np.random.random(n)#随机点x坐标组
rx = 80*(rx-0.5)
ry = np.random.random(n)#随机点y坐标组
ry = 80*(ry-0.5)
rz = np.random.random(n)#随机点z坐标组
rz = 80*(rz-0.5)
t = np.random.random(n)#波函数比较点

#径向波函数
def R(r):
    R = (192-144*r+24*r**2-r**3)*np.exp(-r/4)
    return R

#球谐函数
def Y(l,m,phi,theta):
    return sp.sph_harm(m,l,phi,theta)

#概率密度函数
def f(l,m,r,phi,theta):
    psi = R(r)*Y(l,m,phi,theta)#波函数
    x = np.conj(psi)*(psi)#概率密度函数
    return x.real

#将随机点的x,y,z坐标带入
ruslt = np.array([])
for i in range(n):
    #下面三行将直角坐标系转化为球坐标系
    r = (rx[i]**2+ry[i]**2+rz[i]**2)**0.5
    theta = np.arccos(rz[i]/((rx[i]**2+ry[i]**2+rz[i]**2)**0.5))
    phi = np.arctan(ry[i]**2/rx[i]**2)
    
    ruslt = np.append(ruslt,f(2,0,r,phi,theta))
#归一化
ruslt = ruslt-ruslt.min()/(ruslt.max() - ruslt.min())
length_ruslt = np.size(ruslt)

position_x = np.array([])
position_y = np.array([])
position_z = np.array([])
position_t = np.array([])

#判断哪些随机点可以保存！！！！！！！！！！！！！！
for i in range(length_ruslt):
    if ruslt[i] > t[i]:
         position_x = np.append(position_x,rx[i])
         position_y = np.append(position_y,ry[i])
         position_z = np.append(position_z,rz[i])
         position_t = np.append(position_t,t[i])
print(np.size(position_x))#查看有几个随机点被保存

#设置画布基础属性
fig = plt.figure(dpi = 150)
ax = fig.add_subplot(111, projection='3d')

#绘制散点图
picture = ax.scatter(position_x,position_y,position_z,s = 1.8, c=position_t,cmap='bwr')

ax.set(facecolor='black',alpha = 1)#设置为黑色
fig.set(facecolor='black')#画布颜色为黑色
bar = plt.colorbar(picture, ax=ax)#颜色条表示
ax.grid(False)#关闭网格
ax.set_axis_off()#关闭轴面

def update(angle):
    #设置轴向旋转:elev沿y周旋转，azim沿z周旋转
    ax.view_init(elev = 6,azim=angle)
    
ani = FuncAnimation(fig, update, frames=np.linspace(0,360,100), interval=50)

#打开下面一行用来保存动态图
#ani.save("3d氢原子散点图.gif",writer='pillow')
plt.show()

可以得到如下动态图：越红则电子出现的概率越大，越蓝则电子出现的概率越小。