I'm trying to fit my thermal conductivity into the Debye-Callaway equation. However, one of my parameters is coming back negative. I've tried different initial guesses. So I'm attaching a code with thermal conductivity data from literature and the values that they got to understand where my model is going wrong. What concerns me is it's not very sensitive to initial parameters. I'd really appreciate some help regarding this. Thank you so much! The parameters that the paper reported were A=6.73E-43 B=5.38E-18 The parameters I got: Fitted parameters: 1.417070339751493e-44, B = 1.075485244260637e-15. The functions are defined exactly the same as in paper:   EDIT: I've added a second piece of code that still returns the second factor negatively. I'm not sure if it's because the factor is getting stuck at a local minimum somewhere? import numpy as np from scipy.integrate import quad from scipy.optimize import curve_fit import matplotlib.pyplot as plt # Constants kb = 1.380649e-23 # Boltzmann constant (J/K) hbar = 1.0545718e-34 # Reduced Planck's constant (J*s) theta_D = 280 # Debye temperature (K) v = 2300 # Average phonon group velocity (m/s) - # Define the inverse relaxation times def tau_PD_inv(omega, A): return A * omega**4 def tau_U_inv(omega, B, T): return B * omega**2 * np.exp(-theta_D / (3 * T)) def tau_GB_inv(omega): D = 0.7e-6 # Grain boundary characteristic length (m) return v / D def tau_total_inv(omega, A, B, T): return tau_PD_inv(omega, A) + tau_U_inv(omega, B, T) + tau_GB_inv(omega) # Debye Callaway model for thermal conductivity def kappa_lattice(T, A, B): def integrand(x): omega = x * kb * T / hbar tau_total = 1 / tau_total_inv(omega, A, B, T) return x**4 * np.exp(x) / (np.exp(x) - 1)**2 * tau_total integral, _ = quad(integrand, 0, theta_D / T) prefactor = kb**4 * T**3 / (2 * np.pi**2 * v * hbar**3) return prefactor * integral # Vectorize the model for curve fitting def kappa_lattice_vectorized(T, A, B ): return np.array([kappa_lattice(Ti, A, B ) for Ti in T]) # Example thermal conductivity and temperature data for Ge WQ T_data = np.array([3.0748730372823156, 3.1962737114707203, 3.4441443065638015, 3.785528026732028, 4.043999986888003, 4.524033859711473, 4.902329259251676, 5.270E+00, 5.902534836759174, 6.721645553370004, 7.690683028568884, 9.031E+00, 1.103E+01, 12.811754075061538, 16.07373842258417, 19.67675447459597, 23.951202854256202, 29.710517518825778, 3.684E+01, 45.030707655028536, 56.92469139605285, 6.929E+01, 85.62796640159988, 105.41803279520694, 129.5096325800293, 159.77670578646953, 200.9314386671018, 240.34130109997898, 291.63127697234717]) # Temperature in K kappa_data = np.array([ 0.5021456742832973, 0.5907837911587943, 0.7331309897165834, 0.9570228747931655, 1.1281866396926783, 1.4904675008427906, 1.8682223847710366, 2.1856825810901013, 2.8515011626042894, 4.057791636291224, 5.328111104734368, 7.645627510566391, 11.045409412114584, 14.277055520040301, 19.54065663430087, 24.738891215975602, 28.776422188514683, 32.562485211039274, 33.37572711709385, 32.455084248250884, 28.56802589740088, 24.937473590597776, 20.887192910896978, 17.37073097716312, 14.783541083880793, 12.26010669954394, 9.991973525044475, 8.531955488387288, 7.347906481962733]) # Thermal conductivity in W/mK # Initial guess for parameters A, B, and C initial_guess = [4.75E-42, 1.054E-17] # Curve fitting popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, maxfev=8_000) #popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, bounds=(0, np.inf)) # Extract fitted parameters A_fit, B_fit = popt # Print the fitted parameters print(f"Fitted parameters: A = {A_fit}, B = {B_fit}") # Generate curve fit data T_fit = np.linspace(min(T_data), max(T_data), 100) kappa_fit = kappa_lattice_vectorized(T_fit, *popt) # Plot the original data, fitted data, and extrapolated data plt.figure(figsize=(8, 6)) plt.plot(T_data, kappa_data, 'o', label='Experimental data') plt.plot(T_fit, kappa_fit, '-', label='Fitted model') #plt.plot(T_extrap, kappa_extrap, '--', label='Extrapolated (400-600K)') plt.xlabel('Temperature (K)') plt.ylabel('Thermal Conductivity (W/mK)') plt.legend() plt.title('Thermal Conductivity Fitting ') plt.show()  import numpy as np from scipy.integrate import quad from scipy.optimize import curve_fit import matplotlib.pyplot as plt # Constants kb = 1.380649e-23 # Boltzmann constant (J/K) hbar = 1.0545718e-34 # Reduced Planck's constant (J*s) theta_D = 483.0 # Debye temperature (K) v = 4618.0 # Average phonon group velocity (m/s) - # Define the inverse relaxation times def tau_PD_inv(omega, A): return A * omega**4 def tau_U_inv(omega, B, T): return B * omega**2 * T * np.exp(-theta_D / (3 * T)) def tau_GB_inv(omega): D = 1e-6 # Grain boundary characteristic length (m) return v / D def tau_new_inv(omega, C): return C * omega**2 def tau_total_inv(omega, A, B, C, T): return tau_PD_inv(omega, A) + tau_U_inv(omega, B, T) + tau_GB_inv(omega) + tau_new_inv(omega, C) # Debye Callaway model for thermal conductivity def kappa_lattice(T, A, B, C): def integrand(x): omega = x * kb * T / hbar tau_total = 1 / tau_total_inv(omega, A, B, C, T) return x**4 * np.exp(x) / (np.exp(x) - 1)**2 * tau_total integral, _ = quad(integrand, 0, theta_D / T) prefactor = kb**4 * T**3 / (2 * np.pi**2 * v * hbar**3) return prefactor * integral # Vectorize the model for curve fitting def kappa_lattice_vectorized(T, A, B, C): return np.array([kappa_lattice(Ti, A, B, C) for Ti in T]) # Example thermal conductivity and temperature data for Ge WQ T_data = np.array([31.21198032, 31.63863233, 33.88415583, 36.02698918, 38.39872854, 39.56437842, 43.09065509, 44.50006634, 46.62080745, 43.20705304, 47.77194375, 51.1374173, 53.33894737, 55.68342411, 57.7839973, 60.44084481, 63.31137245, 66.34735373, 69.16983511, 74.09615366, 77.63998236, 82.66616204, 87.95740211, 93.19748533, 99.59297955, 105.6903923, 113.1334233, 120.6776065, 127.0279836, 134.1693093, 141.3911575, 148.1783285, 155.5512459, 162.631427, 169.7085111, 177.378102, 183.5898785, 190.9367188, 197.6917958, 205.7338436, 212.8065032, 220.0698325, 226.7609477, 234.0571384, 241.1287763, 247.8786187, 255.2722521, 262.3425521, 269.4163916, 276.4872816, 283.5510925, 297.1613034, 304.7919793, 311.8668132, 318.9352276, 326.0927208, 333.0893085, 340.164332, 347.2375209, 354.9542919, 362.0246034, 369.0955138, 375.8247965, 383.5799572]) # Temperature in K kappa_data = np.array([ 1.085444844, 1.13592882, 1.176429149, 1.217086899, 1.258678341, 1.305290878, 1.357081371, 1.398857871, 1.26084173, 1.436622858, 1.483042548, 1.528031021, 1.561551609, 1.598710119, 1.632418606, 1.665656346, 1.703762245, 1.742520105, 1.774186181, 1.820109169, 1.855169951, 1.893680154, 1.932568555, 1.962172074, 1.999232142, 2.029595098, 2.0490363, 2.067859052, 2.092450212, 2.119053495, 2.131413546, 2.145488003, 2.158842932, 2.169084281, 2.169487956, 2.184424479, 2.196747698, 2.20620033, 2.217744549, 2.223537925, 2.236483003, 2.258522251, 2.273748326, 2.278514432, 2.28270469, 2.288247862, 2.291668986, 2.298648488, 2.302662675, 2.299560106, 2.311410793, 2.339970141, 2.347949378, 2.349474819, 2.354437976, 2.369869851, 2.378039597, 2.384365082, 2.39183344, 2.395266095, 2.399300875, 2.40766296, 2.409888354, 2.442162261]) # Thermal conductivity in W/mK # Initial guess for parameters A, B, and C initial_guess = [4.75E-43, 3.6979E-19, 1e-16] # Curve fitting popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess,maxfev=8_000) #popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, bounds=(0, np.inf)) # Extract fitted parameters A_fit, B_fit, C_fit = popt # Print the fitted parameters print(f"Fitted parameters: A = {A_fit}, B = {B_fit}, C = {C_fit}") # Generate curve fit data T_fit = np.linspace(min(T_data), max(T_data), 100) kappa_fit = kappa_lattice_vectorized(T_fit, *popt) # Plot the original data, fitted data, and extrapolated data plt.figure(figsize=(8, 6)) plt.plot(T_data, kappa_data, 'o', label='Experimental data') plt.plot(T_fit, kappa_fit, '-', label='Fitted model') #plt.plot(T_extrap, kappa_extrap, '--', label='Extrapolated (400-600K)') plt.xlabel('Temperature (K)') plt.ylabel('Thermal Conductivity (W/mK)') plt.legend() plt.title('Thermal Conductivity Fitting ') plt.show() 
Your tau_U_inv function forgot to take into account T. The corrected version should be this: def tau_U_inv(omega, B, T): return B * omega**2 * T * np.exp(-theta_D / (3 * T))  When I make this change the parameters show Fitted parameters: A = 1.935907419814169e-43, B = 9.770738747258265e-18  which appears to be much closer to the original.