绘制距离/角度函数形式的谐振势时绘图异常，求代码排查

我帮你梳理下代码里的几个关键疏漏，这些应该就是导致绘图异常的核心原因：

1. 错误将原子索引当作坐标计算距离 #

从bonds.dat读取的b_i和b_j应该是原子的编号（索引），而不是原子的空间坐标。你直接用r = b_j - b_i计算得到的不是原子间的实际距离，只是索引的差值，这完全不符合谐振势公式里r（原子间距离）的物理意义，自然会得到异常的绘图结果。

2. 数组维度不匹配，未针对每个键单独计算势能 #

r_0和k_b是从文件中读取的数组（每个键对应一组平衡距离和力常数），你直接将整个数组传入U_bonds函数时，会和r（哪怕是正确的距离数组）产生维度不匹配的问题，NumPy的广播机制可能会生成错误的结果，甚至直接抛出错误。每个键的势能需要对应自己的r_0和k_b来计算。

3. 缺少原子坐标数据的读取与处理 #

要计算原子间的实际距离，你需要一个包含原子空间坐标的数据文件（比如coords.dat，通常每行包含原子索引、x、y、z坐标），然后根据b_i和b_j的索引去提取对应原子的坐标，再计算欧几里得距离。

修正后的代码示例 #

假设你有一个coords.dat文件，格式如下（每行是原子索引、x坐标、y坐标、z坐标）：

1 0.0 0.0 0.0
2 1.0 0.0 0.0
3 0.5 1.0 0.0
...

修正后的代码如下：

import numpy as np
import matplotlib.pyplot as plt

# 读取键数据：原子对索引、平衡距离、力常数
b_i, b_j, r_0, k_b = np.loadtxt('bonds.dat', unpack=True)
# 读取原子坐标：原子索引、x、y、z
atom_indices, x, y, z = np.loadtxt('coords.dat', unpack=True)

# 定义键势能函数
def U_bonds(r, r0, kb):
    return 0.5 * kb * (r - r0)**2

# 计算每个键的实际距离
bond_distances = []
for i, j in zip(b_i, b_j):
    # 找到对应原子的坐标（注意索引可能从1开始，需要转成0索引）
    idx_i = np.where(atom_indices == i)[0][0]
    idx_j = np.where(atom_indices == j)[0][0]
    # 计算欧几里得距离
    dx = x[idx_j] - x[idx_i]
    dy = y[idx_j] - y[idx_i]
    dz = z[idx_j] - z[idx_i]
    dist = np.sqrt(dx**2 + dy**2 + dz**2)
    bond_distances.append(dist)

bond_distances = np.array(bond_distances)
# 计算每个键的势能
bond_potentials = U_bonds(bond_distances, r_0, k_b)

# 绘制键距离与势能的关系
plt.scatter(bond_distances, bond_potentials, label='Bond Potentials')
plt.xlabel('Interatomic Distance (r)')
plt.ylabel('Resonant Potential (U)')
plt.title('Bond Resonant Potential vs Interatomic Distance')
plt.legend()
plt.show()

关于角度势能的补充 #

如果要绘制角度的谐振势，你需要用类似的方法：读取原子坐标，根据a_i、a_j、a_k三个原子的坐标计算键角（theta），再代入U_angles函数计算势能，最后绘图。计算键角可以用向量点积公式：

def calculate_angle(coord_i, coord_j, coord_k):
    # 计算向量 ja 和 jk
    vec_ja = coord_i - coord_j
    vec_jk = coord_k - coord_j
    # 点积计算角度（弧度转角度可选）
    cos_theta = np.dot(vec_ja, vec_jk) / (np.linalg.norm(vec_ja) * np.linalg.norm(vec_jk))
    theta = np.arccos(np.clip(cos_theta, -1.0, 1.0))
    return theta

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