From b0895ac1fa76ea46d4cffe040826982178018d05 Mon Sep 17 00:00:00 2001 From: Marvin Poul Date: Tue, 8 Oct 2024 13:39:49 +0200 Subject: [PATCH 1/2] Add new class to calculate electronic F from DFT The new class FermiSmearing is accessible via job.get_electronic_structure().fermi for all GenericDFTJobs. it calculates the free energy of the electronic system and can output it as a pandas dataframe. It also calculates the density at the fermi level for comparison with the Sommerfeld method, though the method itself is not implemented. --- pyiron_atomistics/dft/waves/electronic.py | 165 ++++++++++++++++++++++ 1 file changed, 165 insertions(+) diff --git a/pyiron_atomistics/dft/waves/electronic.py b/pyiron_atomistics/dft/waves/electronic.py index 1b13bd576..fdc804327 100644 --- a/pyiron_atomistics/dft/waves/electronic.py +++ b/pyiron_atomistics/dft/waves/electronic.py @@ -4,7 +4,16 @@ from __future__ import print_function +from typing import Iterable +import dataclasses + import numpy as np +import numpy.typing as npt +import pandas as pd +import scipy.optimize as so +from scipy.constants import Boltzmann, eV +kB = Boltzmann / eV +import scipy.special from pyiron_atomistics.atomistics.structure.atoms import ( Atoms, @@ -748,6 +757,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 +868,156 @@ 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 energies to obtain electronic free + energies. + + See [1] for an introduction to the theory. + + `eigenvalues` and `weights` must be in the following shape (spin channels, kpoints, bands). + + [1]: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.95.165126 + + Attributes: + eigenvalues (numpy.array of floats): the energy eigenvalues from a DFT calculations + weights (numpy.array of floats): the k point weights from a DFT calculations + n_electrons (int): the number of electrons in the underlying structure + """ + eigenvalues: npt.NDArray[float] # shape of (#spin channels, #kpoints, #bands) + weights: npt.NDArray[float] # shape of (#spin channels, #kpoints, #bands) + n_electrons: int + + def __post_init__(self): + self.eigenvalues.setflags(write=False) + self.weights.setflags(write=False) + + @classmethod + def from_vasp(cls, job: "Vasp") -> "Self": + if not job.status.finished: + raise ValueError(f"Can only instantiate {cls.__name__} from finished jobs!") + 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()) + + @classmethod + def from_sphinx(cls, job: "Sphinx") -> "Self": + if not job.status.finished: + raise ValueError(f"Can only instantiate {cls.__name__} from finished jobs!") + eigenvals = job.content['output/electronic_structure/eig_matrix'] + weights = job.content['output/electronic_structure/k_weights'].reshape(1, -1, 1) + weights *= 3 - eigenvals.shape[0] + occ = job.content['output/electronic_structure/occ_matrix'] + n_elect = int((weights * occ).sum()) + return cls(eigenvals, weights, n_elect) + + @classmethod + def from_electronic_structure(cls, electronic_structure: ElectronicStructure) -> "Self": + eigenvals = electronic_structure.eigenvalue_matrix.copy() + weights = electronic_structure.kpoint_weights.copy().reshape(1, -1, 1) + weights *= 3 - eigenvals.shape[0] + n_elect = int((weights * electronic_structure.occupancy_matrix).sum()) + return cls(eigenvals, weights, n_elect) + + @staticmethod + def _efermi(e, T, mu): + if T != 0: + # avoid overflow warnings by restricting range of exponent, 1/1+e-/+100 is 0/1 anyway. + exponent = np.clip((e-mu)/kB/T, -100, 100) + return 1 / (1 + np.exp(exponent)) + 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.root_scalar(lambda mu: self.occupation(T, mu) - self.n_electrons, + bracket=(self.eigenvalues.min(), self.eigenvalues.max()), + xtol=1e-5) + ef = ret.root + occ = self.occupation(T, ef) + if not ret.converged or (ret.converged and abs(occ - self.n_electrons) > 0.01): + logger.warning(f"Failed to find fermi level at {T}K! Got {self.occupation(T, ef)} electrons " + f"instead of {self.n_electrons}!") + return ef + + def density_at_fermi_level( + self, T: float | Iterable[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_fermi(Ti, smear) for Ti in T]) + ef = self.find_fermi_level(T) + smear = kB*T + # avoid overflow warnings by restricting range of exponent, 1/1+cosh-/+100 is 0 anyway. + ediff = np.clip((self.eigenvalues - ef)/smear, -100, 100) + # 1/cosh+1 is the "delta" function corresponding to Fermi smearing, see: + # https://doi.org/10.1103/PhysRevB.107.195122 + return 1/smear * (self.weights * ( 1/2 * 1/(np.cosh(ediff) + 1) )).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*(scipy.special.entr(f) + scipy.special.entr(1-f)) ).sum() + F = U - T * S + return { + 'F': F, + 'U': U, + 'S': S, + 'T': T, + 'mu': ef, + } def electronic_structure_dict_to_hdf(data_dict, hdf, group_name): with hdf.open(group_name) as h_es: From cdd9391c71864cacc3b02bc1ffd5ebce3514de6c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:35:28 +0000 Subject: [PATCH 2/2] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- pyiron_atomistics/dft/waves/electronic.py | 72 ++++++++++++++--------- 1 file changed, 45 insertions(+), 27 deletions(-) diff --git a/pyiron_atomistics/dft/waves/electronic.py b/pyiron_atomistics/dft/waves/electronic.py index fdc804327..a551c362f 100644 --- a/pyiron_atomistics/dft/waves/electronic.py +++ b/pyiron_atomistics/dft/waves/electronic.py @@ -4,14 +4,15 @@ from __future__ import print_function -from typing import Iterable import dataclasses +from typing import Iterable import numpy as np import numpy.typing as npt import pandas as pd import scipy.optimize as so from scipy.constants import Boltzmann, eV + kB = Boltzmann / eV import scipy.special @@ -868,6 +869,7 @@ def resolved_dos_matrix(self): def resolved_dos_matrix(self, val): self._resolved_dos_matrix = val + @dataclasses.dataclass(frozen=True) class FermiSmearing: """ @@ -885,8 +887,9 @@ class FermiSmearing: weights (numpy.array of floats): the k point weights from a DFT calculations n_electrons (int): the number of electrons in the underlying structure """ - eigenvalues: npt.NDArray[float] # shape of (#spin channels, #kpoints, #bands) - weights: npt.NDArray[float] # shape of (#spin channels, #kpoints, #bands) + + eigenvalues: npt.NDArray[float] # shape of (#spin channels, #kpoints, #bands) + weights: npt.NDArray[float] # shape of (#spin channels, #kpoints, #bands) n_electrons: int def __post_init__(self): @@ -897,9 +900,9 @@ def __post_init__(self): def from_vasp(cls, job: "Vasp") -> "Self": if not job.status.finished: raise ValueError(f"Can only instantiate {cls.__name__} from finished jobs!") - eigenvals = job.content['output/generic/dft/bands/eig_matrix'] + 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) + 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()) @@ -908,15 +911,17 @@ def from_vasp(cls, job: "Vasp") -> "Self": def from_sphinx(cls, job: "Sphinx") -> "Self": if not job.status.finished: raise ValueError(f"Can only instantiate {cls.__name__} from finished jobs!") - eigenvals = job.content['output/electronic_structure/eig_matrix'] - weights = job.content['output/electronic_structure/k_weights'].reshape(1, -1, 1) + eigenvals = job.content["output/electronic_structure/eig_matrix"] + weights = job.content["output/electronic_structure/k_weights"].reshape(1, -1, 1) weights *= 3 - eigenvals.shape[0] - occ = job.content['output/electronic_structure/occ_matrix'] + occ = job.content["output/electronic_structure/occ_matrix"] n_elect = int((weights * occ).sum()) return cls(eigenvals, weights, n_elect) @classmethod - def from_electronic_structure(cls, electronic_structure: ElectronicStructure) -> "Self": + def from_electronic_structure( + cls, electronic_structure: ElectronicStructure + ) -> "Self": eigenvals = electronic_structure.eigenvalue_matrix.copy() weights = electronic_structure.kpoint_weights.copy().reshape(1, -1, 1) weights *= 3 - eigenvals.shape[0] @@ -927,7 +932,7 @@ def from_electronic_structure(cls, electronic_structure: ElectronicStructure) -> def _efermi(e, T, mu): if T != 0: # avoid overflow warnings by restricting range of exponent, 1/1+e-/+100 is 0/1 anyway. - exponent = np.clip((e-mu)/kB/T, -100, 100) + exponent = np.clip((e - mu) / kB / T, -100, 100) return 1 / (1 + np.exp(exponent)) else: return 1 * (e <= mu) @@ -955,18 +960,22 @@ def find_fermi_level(self, T: float) -> float: Returns: float: chemical potential or fermi level """ - ret = so.root_scalar(lambda mu: self.occupation(T, mu) - self.n_electrons, - bracket=(self.eigenvalues.min(), self.eigenvalues.max()), - xtol=1e-5) + ret = so.root_scalar( + lambda mu: self.occupation(T, mu) - self.n_electrons, + bracket=(self.eigenvalues.min(), self.eigenvalues.max()), + xtol=1e-5, + ) ef = ret.root occ = self.occupation(T, ef) if not ret.converged or (ret.converged and abs(occ - self.n_electrons) > 0.01): - logger.warning(f"Failed to find fermi level at {T}K! Got {self.occupation(T, ef)} electrons " - f"instead of {self.n_electrons}!") + logger.warning( + f"Failed to find fermi level at {T}K! Got {self.occupation(T, ef)} electrons " + f"instead of {self.n_electrons}!" + ) return ef def density_at_fermi_level( - self, T: float | Iterable[float] + self, T: float | Iterable[float] ) -> float | npt.NDArray[float]: """ Calculate the electronic density at the fermi level. @@ -974,14 +983,16 @@ def density_at_fermi_level( if isinstance(T, Iterable): return np.array([self.density_at_fermi_level_fermi(Ti, smear) for Ti in T]) ef = self.find_fermi_level(T) - smear = kB*T + smear = kB * T # avoid overflow warnings by restricting range of exponent, 1/1+cosh-/+100 is 0 anyway. - ediff = np.clip((self.eigenvalues - ef)/smear, -100, 100) + ediff = np.clip((self.eigenvalues - ef) / smear, -100, 100) # 1/cosh+1 is the "delta" function corresponding to Fermi smearing, see: # https://doi.org/10.1103/PhysRevB.107.195122 - return 1/smear * (self.weights * ( 1/2 * 1/(np.cosh(ediff) + 1) )).sum() + return 1 / smear * (self.weights * (1 / 2 * 1 / (np.cosh(ediff) + 1))).sum() - def energy_entropy(self, T: float | Iterable[float]) -> dict[str, float] | pd.DataFrame: + def energy_entropy( + self, T: float | Iterable[float] + ) -> dict[str, float] | pd.DataFrame: """ Calculate the electronic free energy properties. @@ -1008,17 +1019,24 @@ def energy_entropy(self, T: float | Iterable[float]) -> dict[str, float] | pd.Da 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*(scipy.special.entr(f) + scipy.special.entr(1-f)) ).sum() + U = ( + ((f - (self.eigenvalues <= ef) * 1) * (self.eigenvalues - ef)) + * self.weights + ).sum() + S = ( + kB + * (self.weights * (scipy.special.entr(f) + scipy.special.entr(1 - f))).sum() + ) F = U - T * S return { - 'F': F, - 'U': U, - 'S': S, - 'T': T, - 'mu': ef, + "F": F, + "U": U, + "S": S, + "T": T, + "mu": ef, } + def electronic_structure_dict_to_hdf(data_dict, hdf, group_name): with hdf.open(group_name) as h_es: for k, v in data_dict.items():