如何修正Python代码绘制晶体场d轨道以匹配标准示意图？

我尝试使用Python绘制晶体场d轨道，编写了如下代码，但生成的d轨道3D图像与标准示意图存在差异，希望能得到代码修正方案，使绘制结果与标准图一致。

原代码 #

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(231, projection='3d')
ax2 = fig.add_subplot(232, projection='3d')
ax3 = fig.add_subplot(233, projection='3d')
ax4 = fig.add_subplot(234, projection='3d')
ax5 = fig.add_subplot(235, projection='3d')


theta = np.linspace(0, 2*np.pi, 100)
phi = np.linspace(0, np.pi, 100)
theta, phi = np.meshgrid(theta, phi)

# Define the radial functions for each orbital
dx2_minus_y2 = np.sqrt(3/np.pi) * np.cos(phi)**2
dz2 = np.sqrt(5/16*np.pi) * (3*np.cos(phi)**2 - 1)
dyz = np.sqrt(15/np.pi) * np.sin(phi) * np.cos(phi) * np.exp(1j*theta)
dxz = np.sqrt(15/np.pi) * np.sin(phi) * np.cos(phi) * np.exp(1j*theta)
dxy = np.sqrt(15/np.pi) * np.sin(phi)**2 * np.exp(2j*theta)


# Convert spherical coordinates to Cartesian coordinates
x = np.sin(phi) * np.cos(theta)
y = np.sin(phi) * np.sin(theta)
z = np.cos(phi)

# colormap
colormap = plt.cm.viridis

# Plot the d-orbitals in separate subplots
ax1.plot_surface(x * (dx2_minus_y2 - y**2), y * (dx2_minus_y2 - y**2), z * (dx2_minus_y2 - y**2), cmap=colormap, alpha=0.8)
ax1.set_title('d_x^2 - y^2')
ax1.axis('off')
ax1.plot_wireframe(x * (dx2_minus_y2 - y**2), y * (dx2_minus_y2 - y**2) , z * (dx2_minus_y2 - y**2) , color='gray', linewidth=0.1)



ax2.plot_surface(x * dz2, y * dz2, z * dz2, cmap=colormap, alpha=0.8)
ax2.set_title('d_z^2')
ax2.axis('off')
ax2.plot_wireframe(x * dz2, y * dz2, z * dz2 , color='gray', linewidth=0.1)



ax3.plot_surface(x * dyz.real, y * dyz.real, z * dyz.real, cmap=colormap, alpha=0.8)
ax3.set_title('d_yz')
ax3.axis('off')
ax3.plot_wireframe(x * dyz.real, y * dyz.real, z * dyz.real , color='gray', linewidth=0.1)


ax4.plot_surface(x * dxz.real, y * dxz.real, z * dxz.real, cmap=colormap, alpha=0.8)
ax4.set_title('d_xz')
ax4.axis('off')
ax4.plot_wireframe(x * dxz.real, y * dxz.real, z * dxz.real, color='gray', linewidth=0.1)



ax5.plot_surface(x * dxy.real, y * dxy.real, z * dxy.real, cmap=colormap, alpha=0.8)
ax5.set_title('d_xy')
ax5.axis('off')
ax5.plot_wireframe(x * dxy.real, y * dxy.real, z * dxy.real, color='gray', linewidth=0.1)



fig.suptitle('Spatial Orientation of d-Orbitals with XYZ Axes', y=1.02)


plt.tight_layout()


plt.show()

当前生成图像 #

标准d轨道示意图 #

代码修正方案 #

原代码的核心问题是球谐函数定义错误、复轨道实部处理错误以及绘图时的径向缩放逻辑错误，以下是修正步骤和完整代码：

错误分析 #

1. 球谐函数的系数和表达式完全错误：比如dz2的系数写成了√(5/16*np.pi)，正确应为√(5/(16π))；dx2_minus_y2、dyz、dxz、dxy的表达式均不符合d轨道的球谐函数形式。
2. 复轨道dyz和dxz定义重复，没有区分开方位角的三角函数项。
3. 绘图时错误地对径向函数做了额外运算（如dx2_minus_y2 - y**2），正确逻辑应为直接用球谐函数的实部/模作为径向距离，乘以球坐标转笛卡尔的x、y、z值。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(10, 8))
axes = [
    fig.add_subplot(231, projection='3d'),
    fig.add_subplot(232, projection='3d'),
    fig.add_subplot(233, projection='3d'),
    fig.add_subplot(234, projection='3d'),
    fig.add_subplot(235, projection='3d')
]

# 定义球坐标网格：theta为方位角(0-2π)，phi为极角(0-π)
theta = np.linspace(0, 2*np.pi, 100)
phi = np.linspace(0, np.pi, 100)
theta, phi = np.meshgrid(theta, phi)

# 球坐标转笛卡尔坐标
x = np.sin(phi) * np.cos(theta)
y = np.sin(phi) * np.sin(theta)
z = np.cos(phi)

# 正确定义d轨道的球谐函数（实部形式，符合标准示意图）
# d_{x²-y²}
d_x2y2 = np.sqrt(15/(16*np.pi)) * np.sin(phi)**2 * np.cos(2*theta)
# d_{z²}
d_z2 = np.sqrt(5/(16*np.pi)) * (3*np.cos(phi)**2 - 1)
# d_{yz}
d_yz = np.sqrt(15/(4*np.pi)) * np.sin(phi) * np.cos(phi) * np.sin(theta)
# d_{xz}
d_xz = np.sqrt(15/(4*np.pi)) * np.sin(phi) * np.cos(phi) * np.cos(theta)
# d_{xy}
d_xy = np.sqrt(15/(16*np.pi)) * np.sin(phi)**2 * np.sin(2*theta)

orbitals = [
    (d_x2y2, 'd$_{x^2-y^2}$'),
    (d_z2, 'd$_{z^2}$'),
    (d_yz, 'd$_{yz}$'),
    (d_xz, 'd$_{xz}$'),
    (d_xy, 'd$_{xy}$')
]

colormap = plt.cm.viridis

for ax, (orbital, title) in zip(axes, orbitals):
    # 用球谐函数值作为径向缩放因子，生成曲面
    ax.plot_surface(x * orbital, y * orbital, z * orbital, cmap=colormap, alpha=0.8)
    ax.plot_wireframe(x * orbital, y * orbital, z * orbital, color='gray', linewidth=0.1)
    ax.set_title(title)
    ax.axis('off')
    # 保持坐标轴比例一致，避免变形
    ax.set_box_aspect([1,1,1])

fig.suptitle('Spatial Orientation of d-Orbitals', y=1.02)
plt.tight_layout()
plt.show()

关键修正点 #

1. 使用标准d轨道球谐函数的实部表达式，确保空间分布与标准图一致。
2. 为每个复轨道选择正确的实部形式（比如d_yz对应sinθ项，d_xz对应cosθ项），避免重复定义。
3. 绘图时直接用球谐函数值缩放笛卡尔坐标，去掉错误的额外运算。
4. 添加set_box_aspect([1,1,1])保证3D坐标轴比例一致，避免轨道形状变形。

内容的提问来源于stack exchange，提问作者Aleja Carmento