Add Fermi resmearing to GenericDFTJob

pyiron_atomistics/dft/job/generic.py

@@ -5,6 +5,7 @@
 import warnings
 
 import numpy as np
+import pandas as pd
 
 from pyiron_atomistics.atomistics.job.atomistic import (
     AtomisticGenericJob,

pyiron_atomistics/dft/waves/electronic.py

@@ -4,7 +4,11 @@
 
 from __future__ import print_function
 
+from typing import Iterable
+import dataclasses
+
 import numpy as np
+import numpy.typing as npt
 
 from pyiron_atomistics.atomistics.structure.atoms import (
     Atoms,
@@ -748,6 +752,12 @@ def __str__(self):
     def __repr__(self):
         return self.__str__()
 
+    @property
+    def fermi(self):
+        """Resmear the single particle energies with a Fermi distribution for
+        different temperatures."""
+        return FermiSmearing.from_electronic_structure(self)
+
 
 class Kpoint(object):
     """
@@ -853,6 +863,132 @@ def resolved_dos_matrix(self):
     def resolved_dos_matrix(self, val):
         self._resolved_dos_matrix = val
 
+@dataclasses.dataclass(frozen=True)
+class FermiSmearing:
+    """
+    Resmears a fixed set of single particle energy spectra to obtain free
+    energies.
+
+    See [1] for a 
+    """
+    eigenvalues: npt.NDArray[float] # shape of (#spin channels, #kpoints, #bands)
+    weights:     npt.NDArray[float] # shape of (#spin channels, #kpoints, #bands)
+    n_electrons: int
+    n_atoms:     int
+    initial_fermi_level: float
+
+    def __post_init__(self):
+        self.eigenvalues.setflags(write=False)
+        self.weights.setflags(write=False)
+
+    @classmethod
+    def from_vasp(cls, job: "Vasp") -> "Self":
+        eigenvals = job.content['output/generic/dft/bands/eig_matrix']
+        # bring the weights in the same shape as eig values
+        weights = job.content['output/electronic_structure/k_weights'].reshape(1, -1, 1)
+        # adjust for (non-) spin polarized
+        weights *= 3 - eigenvals.shape[0]
+        return cls(
+            eigenvals, weights, 
+            job.get_nelect(),
+            len(job.structure),
+            job.content['output/electronic_structure/efermi']
+        )
+
+    def _efermi(e, T, mu):
+        if T != 0:
+            return 1 / (1 + np.exp(( (e-mu)/kB/T)))
+        else:
+            return 1 * (e <= mu)
+
+    def occupation(self, T: float, mu: float) -> float:
+        """
+        Calculate the number of electrons in the system.
+
+        Args:
+            T (float): temperature
+            mu (float): chemical potential of electrons
+
+        Returns:
+            float: number of electrons, can be fractional!
+        """
+        return (self._efermi(self.eigenvalues, T, mu) * self.weights).sum()
+
+    def find_fermi_level(self, T: float) -> float:
+        """
+        Find the chemical potential that fills the system with the correct number of electrons.
+
+        Args:
+            T (float): temperature
+
+        Returns:
+            float: chemical potential or fermi level
+        """
+        ret = so.minimize_scalar(lambda mu: (self.occupation(T, mu) - self.n_electrons)**2, 
+                                 bracket=(self.eigenvalues.min(), self.eigenvalues.max()), 
+                                 tol=1e-7)
+        ef = ret.x
+        if ret.success:
+            if abs(self.occupation(T, ef) - self.n_electrons) > 0.001:
+                print(f'Got {self.occupation(T, ef)} instead of {self.n_electrons}!')
+        else:
+            ef = self.fermi
+        return ef
+
+    def density_at_fermi_level(
+            self, T: float | Iterable[float], smear: float
+    ) -> float | npt.NDArray[float]:
+        """
+        Calculate the electronic density at the fermi level.
+        """
+        if isinstance(T, Iterable):
+            return np.array([self.density_at_fermi_level(Ti, smear) for Ti in T])
+        ef = self.find_fermi_level(T)
+        prefactor = 1/np.sqrt(np.pi)/smear
+        ediff = self.eigenvalues - ef
+        return prefactor * (self.weights * np.exp(-(ediff/smear)**2)).sum()
+
+    def energy_entropy(self, T: float | Iterable[float]) -> dict[str, float] | pd.DataFrame:
+        """
+        Calculate the electronic free energy properties.
+
+        Returns a dictionary with the keys:
+            F: electronic free energy per structure
+            U: electronic internal energy per structure
+            S: electronic entropy per structure
+            f, u, s: corresponding per atomic quantities
+            N: number of atoms
+            T: temperature
+            mu: fermi level
+
+        If `T` is an iterable, call this method for every individual temperature
+        and wrap it in a data frame with the corresponding columns.
+
+        Args:
+            T (float, iterable of float): temperature or list of them
+
+        Returns:
+            dict: electronic properties at given temperature
+            DataFrame: electronic properties at given temperatures
+        """
+        if isinstance(T, Iterable):
+            return pd.DataFrame([self.energy_entropy(Ti) for Ti in T])
+        ef = self.find_fermi_level(T)
+        f = self._efermi(self.eigenvalues, T, ef)
+        U = ( ((f - (self.eigenvalues <= ef)*1) * (self.eigenvalues - ef)) * self.weights ).sum()
+        S = kB * ( self.weights*(entr(f) + entr(1-f)) ).sum()
+        F = U - T * S
+        return {
+            'F': F,
+            'U': U,
+            'S': S,
+            'f': F/self.n_atoms,
+            'u': U/self.n_atoms,
+            's': S/self.n_atoms,
+            'N': self.n_atoms,
+            'T': T,
+            'mu': ef,
+        }
 
 def electronic_structure_dict_to_hdf(data_dict, hdf, group_name):
     with hdf.open(group_name) as h_es:

test_integration/tests_sphinx_sphinx_check_all.py

@@ -68,7 +68,7 @@ def test_Fe_ferro_C(self):
         job = self.project.create_job(self.project.job_type.Sphinx, "spx_Fe_ferro_C")
         job.structure = self.project.create.structure.ase.bulk("Fe", a=self.a_Fe)
         job.structure.set_initial_magnetic_moments(len(job.structure) * [2])
-        job.structure += self.project.create_atoms(
+        job.structure += self.project.create.structure.atoms(
             elements=["C"], positions=[[0, 0, 0.5 * self.a_Fe]], magmoms=[0]
         )
         job.calc_static()