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'

def thermal_expansion_coeff(self, structure, temperature, mode="debye"):
        """
        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'
        """
        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)