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.

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()

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.

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()

• python
• scipy
• curve-fitting
• scipy-optimize
• quad

Share Improve this question Follow edited Feb 23 at 15:05 desertnaut 60.5k3232 gold badges155155 silver badges182182 bronze badges asked Jan 29 at 18:37 MaterialsgirlMaterialsgirl 4311 silver badge44 bronze badges 3

• 1 Your tau_U_inv is incorrect. You fot to account for T with the initial multiplication. – Chrispresso Commented Jan 29 at 18:42
• 1 You should also use np.expm1(x)**2 which is a lot more accurate than the explicit algebraic formula when x is small. Numerical analysis is all about the details and handling edge cases very carefully. More info on np.expm1 – Martin Brown Commented Jan 29 at 22:19
• It seems like you did take a stab at adding bounds, but then commented them out. Why? If it's a problem that your decision variable goes negative, then... bound it. – Reinderien Commented Jan 30 at 0:31

Add a comment  | 

2 Answers 2

Sorted by: Reset to default 7

Your tau_U_inv function fot 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.

An easy way to fix your issue (sensitive to initial guess, negative fit non-physical fits) is to fit the log of your parameters - i.e fit lA = log(A), lB = log(B), lC = log(C)

This solves:
- Large exponent variations, which are numerically unstable
- Guaranteed positive A, B, C values
- No initial guess necessary (i.e. scipy.optimize default works [1,1,1])

Here you go (your second code part)

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

# Already fixed 
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])

def kappa_lattice_vectorized_log(T, lA, lB, lC):
    return np.log(kappa_lattice_vectorized(T, np.exp(lA), np.exp(lB), np.exp(lC)))

# 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

# No need for initial guess
### Initial guess for parameters A, B, and C
### initial_guess = [np.log(4.75E-41), np.log(3.6979E-19), np.log(3e-15)]

# Curve fitting
popt, pcov = curve_fit(kappa_lattice_vectorized_log, T_data, kappa_data, maxfev=8_000)
# popt, pcov = curve_fit(kappa_lattice_vectorized_log, T_data, kappa_data, p0=initial_guess, maxfev=8_000)
# popt, pcov = curve_fit(kappa_lattice_vectorized_log, T_data, kappa_data, p0=initial_guess, bounds=(0, 1e-13))

# Extract fitted parameters
lA_fit, lB_fit, lC_fit = popt

# Print the fitted parameters
result_text = f"Fitted parameters: A = {np.exp(lA_fit):.2e}, B = {np.exp(lB_fit):.2e}, C = {np.exp(lC_fit):.2e}"
print(result_text)

# Generate curve fit data
T_fit = np.linspace(min(T_data), max(T_data), 100)
kappa_fit = kappa_lattice_vectorized_log(T_fit, *popt)
# kappa_initial = kappa_lattice_vectorized_log(T_fit, *initial_guess)

# 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_initial, '--g', label='Initial guess')
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\n' + result_text)
plt.show()

Comes out beautifully (Fit parameters in the title):