code: def integrate_gamma(v, v0, gamma0, q0, q1, theta0):
      
      def f_integrand(v):
          gamma = gamma0 * np.exp(q0 / q1 * ((v / v0) ** q1 - 1.))
          return gamma / v
  
      theta_term = quad(f_integrand, v0, v)[0]
      return theta_term
docs: internal function to calculate Debye temperature
  
      :param v: unit-cell volume in A^3
      :param v0: unit-cell volume in A^3 at 1 bar
      :param gamma0: Gruneisen parameter at 1 bar
      :param q0: logarithmic derivative of Gruneisen parameter
      :param q1: logarithmic derivative of Gruneisen parameter
      :param theta0: Debye temperature at 1 bar in K
      :return: Debye temperature in K