code: def debye_temperature(self, volume):
          
          term1 = (2./3. * (1. + self.poisson) / (1. - 2. * self.poisson))**1.5
          term2 = (1./3. * (1. + self.poisson) / (1. - self.poisson))**1.5
          f = (3. / (2. * term1 + term2))**(1. / 3.)
          debye = 2.9772e-11 * (volume / self.natoms) ** (-1. / 6.) * f * \
                 np.sqrt(self.bulk_modulus/self.avg_mass)
          if self.anharmonic_contribution:
              gamma = self.gruneisen_parameter(0, self.ev_eos_fit.v0)  # 0K equilibrium Gruneisen parameter
              return debye * (self.ev_eos_fit.v0 / volume) ** (gamma)
          else:
              return debye
docs: Calculates the debye temperature.
          Eq(6) in doi.org/10.1016/j.comphy.2003.12.001. Thanks to Joey.
  
          Eq(6) above is equivalent to Eq(3) in doi.org/10.1103/PhysRevB.37.790
          which does not consider anharmonic effects. Eq(20) in the same paper
          and Eq(18) in doi.org/10.1016/j.commatsci.2009.12.006 both consider
          anharmonic contributions to the Debye temperature through the Gruneisen
          parameter at 0K (Gruneisen constant).
  
          The anharmonic contribution is toggled by setting the anharmonic_contribution
          to True or False in the QuasiharmonicDebyeApprox constructor.
  
          Args:
              volume (float): in Ang^3
  
          Returns:
              float: debye temperature in K