code: def vibrational_internal_energy(self, temperature, volume):
          
          y = self.debye_temperature(volume) / temperature
          return self.kb * self.natoms * temperature * (9./8. * y +
                                                        3*self.debye_integral(y))
docs: Vibrational internal energy, U_vib(V, T).
          Eq(4) in doi.org/10.1016/j.comphy.2003.12.001
  
          Args:
              temperature (float): temperature in K
              volume (float): in Ang^3
  
          Returns:
              float: vibrational internal energy in eV