diff --git a/src/hwave/solver/uhfr.py b/src/hwave/solver/uhfr.py index 000d0e8..6c65032 100644 --- a/src/hwave/solver/uhfr.py +++ b/src/hwave/solver/uhfr.py @@ -881,12 +881,6 @@ def _calc_energy(self): _green_list = self.green_list Ene = {} Ene["band"] = 0 - # Zero temperature case - sum up energies of occupied states - if self.T == 0: - Ene["band"] = _calc_zero_temp_energy(_green_list) - else: - Ene["band"] = _calc_finite_temp_energy(_green_list) - def _calc_zero_temp_energy(self, green_list): """Calculate band energy for zero temperature case. @@ -915,6 +909,39 @@ def _calc_zero_temp_energy(self, green_list): energy += np.sum(eigenvalue[:occupied_number]) return energy + def _calc_log_terms(self, eigenvalue, mu): + """Calculate logarithmic terms in grand potential. + + Parameters + ---------- + eigenvalue : ndarray + Array of eigenvalues + mu : float + Chemical potential + + Returns + ------- + ndarray + Array of logarithmic terms + + Notes + ----- + Calculates ln(1 + exp(-(e-mu)/T)) for each eigenvalue. + Uses cutoff to avoid numerical overflow: + - If -(e-mu)/T < cutoff: use log1p for numerical stability + - If -(e-mu)/T >= cutoff: approximate as -(e-mu)/T + """ + ln_Ene = np.zeros(eigenvalue.shape) + for idx, value in enumerate(eigenvalue): + # Check if exponential will overflow + if -(value - mu) / self.T < self.ene_cutoff: + # Use log1p for numerical stability + ln_Ene[idx] = np.log1p(np.exp(-(value - mu) / self.T)) + else: + # For large negative arguments, approximate as -(e-mu)/T + ln_Ene[idx] = -(value - mu) / self.T + return ln_Ene + def _calc_finite_temp_energy(self, green_list): """Calculate band energy for finite temperature case. @@ -942,6 +969,7 @@ def _calc_finite_temp_energy(self, green_list): eigenvec = self.green_list[k]["eigenvector"] # Get eigenvectors mu = self.green_list[k]["mu"] # Chemical potential + # Calculate Fermi-Dirac occupations fermi = self._fermi(mu, eigenvalue) # Calculate logarithmic terms for entropy @@ -951,40 +979,16 @@ def _calc_finite_temp_energy(self, green_list): # Add mu*N term and entropy term energy += mu*np.sum(tmp_n) - self.T * np.sum(ln_Ene) - return energy - - def _calc_log_terms(self, eigenvalue, mu): - """Calculate logarithmic terms in grand potential. - - Parameters - ---------- - eigenvalue : ndarray - Array of eigenvalues - mu : float - Chemical potential - - Returns - ------- - ndarray - Array of logarithmic terms - - Notes - ----- - Calculates ln(1 + exp(-(e-mu)/T)) for each eigenvalue. - Uses cutoff to avoid numerical overflow: - - If -(e-mu)/T < cutoff: use log1p for numerical stability - - If -(e-mu)/T >= cutoff: approximate as -(e-mu)/T - """ - ln_Ene = np.zeros(eigenvalue.shape) - for idx, value in enumerate(eigenvalue): - # Check if exponential will overflow - if -(value - mu) / self.T < self.ene_cutoff: - # Use log1p for numerical stability - ln_Ene[idx] = np.log1p(np.exp(-(value - mu) / self.T)) - else: - # For large negative arguments, approximate as -(e-mu)/T - ln_Ene[idx] = -(value - mu) / self.T - return ln_Ene + return energy + + + # Zero temperature case - sum up energies of occupied states + if self.T == 0: + Ene["band"] = _calc_zero_temp_energy(_green_list) + else: + Ene["band"] = _calc_finite_temp_energy(_green_list) + + Ene["InterAll"] = 0 green_local = self.Green.reshape((2 * self.Nsize) ** 2)