code: def thermal_expansion_coeff(self, structure, temperature, mode="debye"):
          
          soec = ElasticTensor(self[0])
          v0 = (structure.volume * 1e-30 / structure.num_sites)
          if mode == "debye":
              td = soec.debye_temperature(structure)
              t_ratio = temperature / td
              integrand = lambda x: (x**4 * np.exp(x)) / (np.exp(x) - 1)**2
              cv = 9 * 8.314 * t_ratio**3 * quad(integrand, 0, t_ratio**-1)[0]
          elif mode == "dulong-petit":
              cv = 3 * 8.314
          else:
              raise ValueError("Mode must be debye or dulong-petit")
          tgt = self.get_tgt(temperature, structure)
          alpha = np.einsum('ijkl,ij', soec.compliance_tensor, tgt)
          alpha *= cv / (1e9 * v0 * 6.022e23)
          return SquareTensor(alpha)
docs: Gets thermal expansion coefficient from third-order constants.
  
          Args:
              temperature (float): Temperature in kelvin, if not specified
                  will return non-cv-normalized value
              structure (Structure): Structure to be used in directional heat
                  capacity determination, only necessary if temperature
                  is specified
              mode (string): mode for finding average heat-capacity,
                  current supported modes are 'debye' and 'dulong-petit'