要计算 ASE 分子动力学轨迹文件中的径向分布函数（RDF）和势能均值函数（PMF），可以使用 ASE 本身结合 NumPy 和 Matplotlib 来完成。以下是一个完整的 Python 示例代码： ```python import numpy as np import matplotlib.pyplot as plt from ase.io import read from ase import Atoms from ase.geometry import get_distances def compute_rdf_and_pmf(traj_file, atom_pair, r_max=5.0, nbins=100, temperature=300): """ 计算 ASE 分子动力学轨迹文件中指定原子对的 RDF 和 PMF。 参数: traj_file (str): ASE 轨迹文件路径（如 'trajectory.xyz'）。 atom_pair (tuple): 原子对（如 ('O', 'H')）。 r_max (float): 最大距离（Å）。 nbins (int): 直方图的 bin 数量。 temperature (float): 温度（单位为 K）。 返回: r (np.array): 径向距离。 rdf (np.array): 径向分布函数。 pmf (np.array): 势能均值函数。 """ # 读取轨迹 traj = read(traj_file, index=':') # 初始化距离列表 distances = [] # 遍历每一帧 for atoms in traj: positions = atoms.positions symbols = np.array(atoms.symbols) # 获取指定原子类型的索引 idx1 = np.where(symbols == atom_pair[0])[0] idx2 = np.where(symbols == atom_pair[1])[0] # 计算不同原子之间的距离 for i in idx1: for j in idx2: if i != j: # 避免同一个原子之间的距离 d = np.linalg.norm(positions[i] - positions[j]) distances.append(d) # 计算 RDF bins = np.linspace(0, r_max, nbins) hist, bin_edges = np.histogram(distances, bins=bins) bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2 volume = atoms.cell.volume # 假设盒子体积恒定 density = len(np.where(symbols == atom_pair[1])[0]) / volume # 第二种原子的密度 rdf = hist / (4 * np.pi * bin_centers**2 * (r_max / nbins) * density * len(traj)) # 计算 PMF kb = 1.380649e-23 # 玻尔兹曼常数 (J/K) pmf = -kb * temperature * np.log(rdf + 1e-10) # 加一个小的常数防止除零错误 return bin_centers, rdf, pmf # 示例使用 traj_file = 'trajectory.xyz' # ASE 轨迹文件路径 atom_pair = ('O', 'H') # 原子对 r, rdf, pmf = compute_rdf_and_pmf(traj_file, atom_pair) # 绘制 RDF 和 PMF plt.figure(figsize=(12, 6)) plt.subplot(1, 2, 1) plt.plot(r, rdf, label='RDF') plt.title('Radial Distribution Function (RDF)') plt.xlabel('Distance (Å)') plt.ylabel('RDF') plt.legend() plt.subplot(1, 2, 2) plt.plot(r, pmf, label='PMF', color='r') plt.title('Potential of Mean Force (PMF)') plt.xlabel('Distance (Å)') plt.ylabel('PMF (J)') plt.legend() plt.tight_layout() plt.show() ``` ### 说明： - 使用 `ase.io.read` 读取轨迹文件（如 `xyz` 或 `traj` 格式）。 - 利用 `numpy.linalg.norm` 计算原子间距离。 - 计算 RDF 时考虑了体积和密度的影响。 - PMF 通过 RDF 计算得出，公式为： $$ \text{PMF}(r) = -k_B T \ln(\text{RDF}(r)) $$ - 使用 `matplotlib` 绘制 RDF 和 PMF 曲线。 ---