-
Notifications
You must be signed in to change notification settings - Fork 125
Description
I am trying to use the GAP model to predict the energy and force of the structure, but unfortunately my python script gives me the following error. I would like to know how to deal with it?This is my python script:
import ase.io
import numpy as np
import matplotlib.pyplot as plt
from quippy.potential import Potential
from quippy.convert import ase_to_quip
from sklearn.metrics import mean_squared_error
import re
读取所有构型
frames = ase.io.read('test.xyz', index=':')
true_energies = []
pred_energies = []
true_forces = []
pred_forces = []
初始化GAP势
try:
gap_pot = Potential(param_str='IP GAP', param_filename='gap_GaO_20220921.xml')
except Exception as e:
raise RuntimeError(f"Failed to load GAP potential: {str(e)}") from e
for atoms in frames:
# 从comment中提取DFT能量
comment = atoms.info.get('comment', '')
energy_match = re.search(r'energy=([-\d.]+)', comment)
if not energy_match:
raise ValueError(f"Energy not found in comment: {comment}")
true_energy = float(energy_match.group(1))
# 获取DFT计算的力
dft_forces = atoms.arrays['forces']
# 转换ASE Atoms为QUIP Atoms对象
quip_atoms = ase_to_quip(atoms)
# 使用GAP势计算
try:
gap_pot.calc(quip_atoms)
pred_energy = quip_atoms.energy
gap_forces = quip_atoms.forces
except Exception as e:
raise RuntimeError(f"GAP calculation failed: {str(e)}") from e
# 存储结果
true_energies.append(true_energy)
pred_energies.append(pred_energy)
true_forces.append(dft_forces.ravel())
pred_forces.append(gap_forces.ravel())
转换为numpy数组
true_forces = np.concatenate(true_forces)
pred_forces = np.concatenate(pred_forces)
计算RMSE
energy_rmse = np.sqrt(mean_squared_error(true_energies, pred_energies))
force_rmse = np.sqrt(mean_squared_error(true_forces, pred_forces))
print(f"Energy RMSE: {energy_rmse:.4f} eV")
print(f"Force RMSE: {force_rmse:.4f} eV/Å")
绘制对角线图
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(true_energies, pred_energies, alpha=0.6)
min_e = min(min(true_energies), min(pred_energies))
max_e = max(max(true_energies), max(pred_energies))
plt.plot([min_e, max_e], [min_e, max_e], 'k--', lw=1)
plt.xlabel('DFT Energy (eV)')
plt.ylabel('GAP Predicted Energy (eV)')
plt.title(f'Energy Comparison (RMSE={energy_rmse:.4f} eV)')
plt.subplot(1, 2, 2)
plt.scatter(true_forces, pred_forces, alpha=0.6, s=2)
min_f = min(min(true_forces), min(pred_forces))
max_f = max(max(true_forces), max(pred_forces))
plt.plot([min_f, max_f], [min_f, max_f], 'k--', lw=1)
plt.xlabel('DFT Force (eV/Å)')
plt.ylabel('GAP Predicted Force (eV/Å)')
plt.title(f'Force Comparison (RMSE={force_rmse:.4f} eV/Å)')
plt.tight_layout()
plt.savefig('gap_vs_dft_comparison.png')
plt.show()