我有一段Python代码,它从文本文件中读取两列数据,并执行曲线拟合过程。代码基于numpy amd数学库,第一次运行代码时,我得到了一组结果。在此运行期间,机器还运行另一个Python代码来进行不同的分析。当我重新运行代码时,我发现得到了截然不同的结果。
为了确认这些结果的正确性,我在MATLAB中重新运行了计算,并且能够重现第二次运行Python代码获得的结果。任何关于为什么会发生这种情况的指示将不胜感激。如果有人需要测试,我可以私下分享数据和代码。
import pandas as pd
import numpy as np
import typing as tp
import math as mt
import matplotlib.pyplot as plt
def curveFitting(energyData: tp.Dict[float, float], order: int) -> None:
pesDataFrame = pd.DataFrame(columns = ["Distance", "Energy"])
for key in energyData:
pesDataFrame.loc[len(pesDataFrame)] = [key, energyData[key]] # Create a DataFrame to drop all energy values above 0 cm-1
pesDataFrame = pesDataFrame[pesDataFrame["Energy"] <= 0] # Drop all energies greater than 0 cm-1
# pesDataFrame = pesDataFrame[pesDataFrame["Distance"] > 1.4]
distances = pesDataFrame["Distance"].to_numpy()
energies = pesDataFrame["Energy"].to_numpy()
rEquilibrium = distances[energies.argmin()]
distances = np.delete(distances, energies.argmin())
energies = np.delete(energies, energies.argmin())
xRvalues = np.zeros(distances.size)
print(rEquilibrium)
for index, value in enumerate(distances):
xRvalues[index] = (pow(value, 3) - pow(rEquilibrium, 3)) / (pow(value, 3) + pow(rEquilibrium, 3)) # Calculate x(r) values according to Eq 14
xMatrix = np.zeros((distances.size, order+1))
for row, value in enumerate(xRvalues):
for ordercount in range(0, order+1):
xMatrix[row][ordercount] = pow(value, ordercount) # Calculate xMatrix according to Eq 19
print(xMatrix.shape)
# term1 = np.zeros(distances.size)
term2 = np.zeros(distances.size)
rescaledEnergy = np.zeros(distances.size)
minEnergy = energies[energies.argmin()]
dissocEnergy = 0 - minEnergy # Calculate dissociation energy as difference between minimum and zero energy.
for index, value in enumerate(energies):
rescaledEnergy[index] = (value - minEnergy) / dissocEnergy
for index, value in enumerate(distances):
# term1[index] = np.emath.log(1 - mt.sqrt(rescaledEnergy[index])) # Eq 11
term2[index] = np.emath.log(1 - mt.sqrt(rescaledEnergy[index])) # Eq 12
# term1[index] = term1[index] / (rEquilibrium - value)
term2[index] = term2[index] / (rEquilibrium - value)
# First calculate beta from Eq 11
beta, _, _, _ = np.linalg.lstsq(xMatrix, term2, rcond = None)
print(beta.shape)
calculateBetas = np.dot(xMatrix, beta)
errorVector = term2 - calculateBetas
print(np.linalg.norm(errorVector))
# calculateBetas = np.zeros(distances.shape)
# for row in range(xMatrix.shape[0]):
# calculateBetas[row] = np.dot(xMatrix[row], beta)
calculateEnergies1 = np.zeros(distances.shape)
# errorVector = np.zeros(distances.shape)
for index, value in enumerate(distances):
distance = value - rEquilibrium
betaterm = mt.exp(-1*calculateBetas[index]*distance)
nextterm = pow(1-betaterm, 2) * dissocEnergy
calculateEnergies1[index] = nextterm + minEnergy
# errorVector[index] = energies[index] - calculateEnergies1[index]
error = np.linalg.norm(errorVector)
plt.plot(distances, energies, "-r", label = "UCCSD(T) Energy")
plt.plot(distances, calculateEnergies1, "-b", label = "Fitted Energies")
plt.title("Fitting Order : {} and Error : {:10.5f}".format(order, error))
plt.xlabel("Interatomic Distance")
plt.ylabel("Energy (cm-1)")
plt.legend()
plt.savefig("uccsdt_fitting_order_{}.png".format(order), dpi = 600)
with open("energies_neutral_631g.dat") as energyFile:
uhfEnergy = {}
uccsdEnergy = {}
uccsdtEnergy = {}
for line in energyFile:
line = line.split()
dist = float(line[0]) # Distance
uhf = float(line[1]) # UHF Energy
uccsd = float(line[2]) # UCCSD Energy
uccsdt = float(line[3]) # UCCSD Energy with UCCSD(T) corrections
uhfEnergy[dist] = uhf
uccsdEnergy[dist] = uccsd
uccsdtEnergy[dist] = uccsdt
for order in range(0,20):
curveFitting(energyData = uccsdtEnergy, order = order)
plt.clf()
here提供了用于拟合的原始数据
1条答案
按热度按时间kb5ga3dv1#
要解决此问题,您可以尝试以下步骤:
验证库的版本和基础环境是否相同。
独立运行代码,而不同时运行其他计算密集型进程。
考虑实施检查或记录,以跟踪曲线拟合过程中的中间结果或任何潜在问题。