Skip to content

Commit

Permalink
change order _calc_**_energy functions
Browse files Browse the repository at this point in the history
  • Loading branch information
k-yoshimi committed Nov 11, 2024
1 parent faf5296 commit a8f6434
Showing 1 changed file with 44 additions and 40 deletions.
84 changes: 44 additions & 40 deletions src/hwave/solver/uhfr.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down

0 comments on commit a8f6434

Please sign in to comment.