code: def get_heat_capacity(self, temperature, structure, n, u, cutoff=1e2):
          
          k = 1.38065e-23
          kt = k*temperature
          hbar_w = 1.05457e-34*self.omega(structure, n, u)
          if hbar_w > kt * cutoff:
              return 0.0
          c = k * (hbar_w / kt) ** 2
          c *= np.exp(hbar_w / kt) / (np.exp(hbar_w / kt) - 1)**2
          return c * 6.022e23
docs: Gets the directional heat capacity for a higher order tensor
          expansion as a function of direction and polarization.
  
          Args:
              temperature (float): Temperature in kelvin
              structure (float): Structure to be used in directional heat
                  capacity determination
              n (3x1 array-like): direction for Cv determination
              u (3x1 array-like): polarization direction, note that
                  no attempt for verification of eigenvectors is made
              cutoff (float): cutoff for scale of kt / (hbar * omega)
                  if lower than this value, returns 0