code: def energy_density(self, strain, convert_GPa_to_eV=True):
          
          e_density = np.sum(self.calculate_stress(strain)*strain) / self.order
          if convert_GPa_to_eV:
              e_density *= self.GPa_to_eV_A3  # Conversion factor for GPa to eV/A^3
          return e_density
docs: Calculates the elastic energy density due to a strain