From 69995b02f114079f4133e35aa060dd51655d5a04 Mon Sep 17 00:00:00 2001 From: Marvin Poul Date: Mon, 7 Oct 2024 16:58:27 +0200 Subject: [PATCH] Add Fermi resmearing to GenericDFTJob --- pyiron_atomistics/dft/job/generic.py | 1 + pyiron_atomistics/dft/waves/electronic.py | 136 ++++++++++++++++++ .../tests_sphinx_sphinx_check_all.py | 2 +- 3 files changed, 138 insertions(+), 1 deletion(-) diff --git a/pyiron_atomistics/dft/job/generic.py b/pyiron_atomistics/dft/job/generic.py index 44fa00a64..adfd08e71 100644 --- a/pyiron_atomistics/dft/job/generic.py +++ b/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, diff --git a/pyiron_atomistics/dft/waves/electronic.py b/pyiron_atomistics/dft/waves/electronic.py index 1b13bd576..2b9254fbe 100644 --- a/pyiron_atomistics/dft/waves/electronic.py +++ b/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: diff --git a/test_integration/tests_sphinx_sphinx_check_all.py b/test_integration/tests_sphinx_sphinx_check_all.py index b6cfab558..15346a041 100644 --- a/test_integration/tests_sphinx_sphinx_check_all.py +++ b/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()