为何使用scipy curve_fit得到的势能面拟合结果偏差极大？

用1991年《J. Chem. Phys.》的《Potential energy surface of H⋅⋅⋅H2O》文章数据测试拟合代码时出现问题：采用和文章完全相同的函数，通过scipy curve_fit拟合后，结果和文章数据、文章拟合曲线偏差极大，甚至不具备这类势能面拟合关键的极小值。

绘图时，多数几何结构的首个数据点因能量值过大导致图表失效被舍去，仅最后一种结构保留全点（该结构下文章和我的拟合效果都很差）。

已经检查过输入数据、重构了mapping函数，甚至设置了严格的参数范围（未知预期结果时没法用这种设置），但都没有效果。现在问题锁定在代码层面，需要让拟合出的能量曲线和文章数据或文章拟合曲线吻合。

from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import math
import scipy
from numpy import array

data = pd.read_excel('data.xlsx', header=None)
values_E = array(data[3])
values_R = array(data[0])
values_theta = array([math.radians(data[1][i]) for i in range(len(data[1]))])
values_phi = array([math.radians(data[2][i]) for i in range(len(data[2]))])
values = array([values_R, values_theta, values_phi])

def mapping(values, eps00, eps10, eps20, eps2_2, sigma00, sigma10, sigma20, sigma2_2):
    return [math.fsum([
    4*eps00*(sigma00/values[0][i])**6*((sigma00/values[0][i])**6 - 1)* 1/2 * np.sqrt(1/np.pi),
    4*eps10*(sigma10/values[0][i])**6*((sigma10/values[0][i])**6 - 1)* 1/2 * np.sqrt(3/np.pi) * np.cos(values[1][i]),
    4*eps20*(sigma20/values[0][i])**6*((sigma20/values[0][i])**6 - 1)* 1/4 * np.sqrt(5/np.pi) * (3 * (np.cos(values[1][i]))**2 - 1),
    4*eps2_2*(sigma2_2/values[0][i])**6*((sigma2_2/values[0][i])**6 - 1)* 1/4 * np.sqrt(15/np.pi) * (np.sin(values[1][i]))**2 * (np.sin(values[2][i]))**2
    ]) for i in range(len(values_R))]
param_bounds=([0,-20, 0, 0, 0, 0, 0, 0], [150,10,10,100,10,10,10,10])
args, _ = curve_fit(mapping, values, values_E, bounds=param_bounds, maxfev=10000000)
eps00, eps10, eps20, eps2_2, sigma00, sigma10, sigma20, sigma2_2 = args

#energy by approximation from the article and by my
article = mapping(values, 107.9, -9.5, 8.6, 35.6, 3.0, 3.3, 2.98, 2.92) 
new_E = mapping(values, eps00, eps10, eps20, eps2_2, sigma00, sigma10, sigma20, sigma2_2)

sp = plt.subplot(231) #90_00
plt.plot(values_R[0:14], values_E[0:14], label='E_90_00')
plt.plot(values_R[1:14], article[1:14], label='article')
plt.plot(values_R[0:14], new_E[0:14], label='new_E')
plt.xlabel('R_90_00')
plt.ylabel('E_90_00')
plt.legend(loc='best')

sp = plt.subplot(232) #90_90
plt.plot(values_R[14:28], values_E[14:28], label='E_90_90')
plt.plot(values_R[15:28], article[15:28], label='article')
plt.plot(values_R[14:28], new_E[14:28], label='new_E')
plt.xlabel('R_90_90')
plt.ylabel('E_90_90')
plt.legend(loc='best')

sp = plt.subplot(233) #00_00
plt.plot(values_R[28:42], values_E[28:42], label='E_00_00')
plt.plot(values_R[29:42], article[29:42], label='article')
plt.plot(values_R[28:42], new_E[28:42], label='new_E')
plt.xlabel('R_00_00')
plt.ylabel('E_00_00')
plt.legend(loc='best')

sp = plt.subplot(234) #180_00
plt.plot(values_R[42:56], values_E[42:56], label='E_180_00')
plt.plot(values_R[43:56], article[43:56], label='article')
plt.plot(values_R[42:56], new_E[42:56], label='new_E')
plt.xlabel('R_180_00')
plt.ylabel('E_180_00')
plt.legend(loc='best')

sp = plt.subplot(235) #127_90
plt.plot(values_R[56:62], values_E[56:62], label='E_127_90')
plt.plot(values_R[56:62], article[56:62], label='article')
plt.plot(values_R[56:62], new_E[56:62], label='new_E')
plt.xlabel('R_127_90')
plt.ylabel('E_127_90')
plt.legend(loc='best')

plt.show()

内容的提问来源于stack exchange，提问作者One-eight greek