r"""
Low-level module containing miscalleneous mathematical functions.
Functions
---------
* :func:`normalize_orbs`: normalize KS orbitals within defined sphere
* :func:`int_sphere`: integral :math:`4\pi \int \mathrm{d}r r^2 f(r)`
* :func:`laplace`: compute the second-order derivative :math:`d^2 y(x) / dx^2`
* :func:`fermi_dirac`: compute the Fermi-Dirac occupation function for given order `n`
* :func:`ideal_entropy`: Define the integrand to be used in :func:`ideal_entropy_int`.
* :func:`ideal_entropy_int`: Compute the entropy for the ideal electron gas (no prefac).
* :func:`fd_int_complete`: compute the complete Fermi-Dirac integral for given order `n`
* :func:`chem_pot`: compute the chemical potential by enforcing charge neutrality
* :func:`f_root_id`: make root input fn for chem_pot with ideal apprx for free electrons
* :func: `lorentzian`: Lorentzian function
"""
# standard libraries
from math import sqrt, pi
import warnings
# external libraries
import numpy as np
from scipy import optimize, integrate
# internal libraries
from . import config
[docs]
def normalize_orbs(eigfuncs_x, xgrid, grid_type):
r"""
Normalize the KS orbitals within the chosen sphere.
Parameters
----------
eigfuncs_x : ndarray
The radial KS eigenfunctions :math:`X_{nl}^{\sigma}(x)`
xgrid : ndarray
The log / sqrt grid over which normalization is performed
Returns
-------
eigfuncs_x_norm : ndarray
The radial KS eigenfunctions normalized over the chosen sphere
"""
# initialize the normalized eigenfunctions
eigfuncs_x_norm = eigfuncs_x
# loop over the eigenfunctions
for n in range(np.shape(eigfuncs_x)[0]):
# compute the mod squared eigenvalues
eigfuncs_sq = eigfuncs_x[n].real ** 2 + eigfuncs_x[n].imag ** 2
# compute the intergal ampsq=4*pi*\int_dr r^2 |R(r)|^2
if grid_type == "log":
ampsq = int_sphere(np.exp(-xgrid) * eigfuncs_sq, xgrid, "log")
else:
ampsq = int_sphere(eigfuncs_sq, xgrid, "sqrt")
# normalize eigenfunctions
eigfuncs_x_norm[n] = eigfuncs_x[n] / sqrt(ampsq)
return eigfuncs_x_norm
[docs]
def int_sphere(fx, xgrid, grid_type):
r"""
Compute integral over sphere defined by input grid.
Parameters
----------
fx : array_like
The function (array) to be integrated
xgrid : ndarray
The log / sqrt radial grid
Returns
-------
I_sph : float
The value of the integrand
Notes
-----
The integral formula on the log grid is given by
.. math:: I = 4 \pi \int \mathrm{d}x\ e^{3x} f(x)
and on the sqrt grid is given by
.. math:: I = 4 \pi \int \mathrm{d}s\ s^5 f(s)
"""
if grid_type == "log":
func_int = 4.0 * pi * np.exp(3.0 * xgrid) * fx
else:
func_int = 8.0 * pi * xgrid**5 * fx
I_sph = np.trapz(func_int, xgrid)
return I_sph
[docs]
def laplace(y, x, axis=-1):
r"""
Compute the second-order derivative :math:`d^2 y(x) / dx^2`.
Derivative can be computed over any given axis.
Parameters
----------
y : ndarray
array y(x) on which laplacian is computed
x : ndarray
x array
axis : int, optional
axis over which derivatives are taken
default : -1
Returns
-------
grad2_y : ndarray
the laplacian of y
"""
# first compute the first-order gradient
grad1_y = np.gradient(y, x, edge_order=2, axis=axis)
# now compute the second-order gradient
grad2_y = np.gradient(grad1_y, x, edge_order=2, axis=axis)
return grad2_y
[docs]
def fermi_dirac(eps, mu, beta, n=0):
r"""
Compute the Fermi-Dirac function, see notes for functional form.
Parameters
----------
mu : array_like
the chemical potential
beta : float
the inverse potential
eps : array_like
the energies
n : int
energy is raised to power n/2 in the numerator (see notes)
Returns
-------
f_fd : array_like
the fermi dirac occupation(s)
Notes
-----
The FD function is defined as:
.. math:: f^{(n)}_{fd}(\epsilon, \mu, \beta) = \frac{\epsilon^{n/2}}{1+\exp(1+
\beta(\epsilon - \mu))}
"""
# dfn the exponential function
# ignore warnings here
with np.errstate(over="ignore"):
fn_exp = np.minimum(np.exp(beta * (eps - mu)), 1e12)
# fermi_dirac dist
f_fd = (eps) ** (n / 2.0) / (1 + fn_exp)
return f_fd
[docs]
def ideal_entropy(eps, mu, beta, n=0):
r"""
Define the integrand to be used in :func:`ideal_entropy_int` (see notes).
Parameters
----------
eps : array_like
the energies
mu : array_like
the chemical potential
beta : float
the inverse potential
n : int
energy is raised to power n/2 in the numerator (see notes)
Returns
-------
f_ent : array_like
the entropy integrand function
Notes
-----
The ideal entropy integrand is defined as
.. math::
f_n(\epsilon,\mu,\beta) = \epsilon^{n/2} (f_\mathrm{fd}\log{f_\mathrm{fd}}
+ (1-f_\mathrm{fd}) \log(1-f_\mathrm{fd}) ),
where :math:`f_\mathrm{fd}=f_\mathrm{fd}(\epsilon,\mu,\beta)` is the Fermi-Dirac
distribution.
"""
# dfn the exponential function
# ignore warnings here
with np.errstate(over="ignore"):
fn_exp = np.minimum(np.exp(beta * (eps - mu)), 1e12)
# the 'raw' Fermi-Dirac distribution
f_fd_raw = 1 / (1 + fn_exp)
# define high and low tolerances for the log function (to avoid nans)
tol_l = 1e-8
tol_h = 1.0 - 1e-8
# first replace the zero values
f_fd_mod = np.where(f_fd_raw > tol_l, f_fd_raw, tol_l)
# now replace the one values
f_fd = np.where(f_fd_mod < tol_h, f_fd_mod, tol_h)
# fermi_dirac dist
f_ent = (eps) ** (n / 2.0) * (f_fd * np.log(f_fd) + (1 - f_fd) * np.log(1 - f_fd))
return f_ent
[docs]
def fd_int_complete(mu, beta, n):
r"""
Compute complete Fermi-Dirac integral for given order (see notes for function form).
Parameters
----------
mu : float
chemical potential
beta : float
inverse temperature
n : int
order of Fermi-Dirac integral (see notes)
Returns
-------
I_n : float
the complete fermi-dirac integral
Notes
-----
Complete Fermi-Dirac integrals are of the form
.. math::
I_{n}(\mu,\beta)=\int_0^\infty\mathrm{d}\epsilon\ \epsilon^{n/2}f_{fd}
(\mu,\epsilon,\beta)
where n is the order of the integral
"""
# use scipy quad integration routine
limup = np.inf
# ignore integration warnings (omnipresent because of inf upper limit)
with warnings.catch_warnings():
warnings.filterwarnings("ignore")
I_n, err = integrate.quad(fermi_dirac, 0, limup, args=(mu, beta, n))
return I_n
[docs]
def ideal_entropy_int(mu, beta, n):
r"""
Compute the entropy for the ideal electron gas (without prefactor) - see notes.
Parameters
----------
mu : float
chemical potential
beta : float
inverse temperature
n : int
order of Fermi-Dirac integral (see notes)
Returns
-------
I_n : float
the complete fermi-dirac integral
Notes
-----
The entropy of an ideal electron gas is defined as
.. math::
I_n(\mu,\beta) = \int_0^\infty \mathrm{d}\epsilon\ \epsilon^{n/2}
(f_\mathrm{fd}\log{f_\mathrm{fd}} + (1-f_\mathrm{fd}) \log(1-f_\mathrm{fd}) ),
where :math:`f_\mathrm{fd}=f_\mathrm{fd}(\epsilon,\mu,\beta)` is the Fermi-Dirac
distribution.
"""
# use scipy quad integration routine
limup = np.inf
# ignore integration warnings (omnipresent because of inf upper limit)
with warnings.catch_warnings():
warnings.filterwarnings("ignore")
I_n, err = integrate.quad(ideal_entropy, 0, limup, args=(mu, beta, n))
return I_n
[docs]
def chem_pot(orbs):
r"""
Determine the chemical potential by enforcing charge neutrality (see notes).
Uses scipy.optimize.root_scalar with brentq implementation.
Parameters
----------
orbs : staticKS.Orbitals
the orbitals object
Returns
-------
mu : array_like
chemical potential (spin-dependent)
Notes
-----
Finds the roots of equation
.. math:: \sum_{nl} (2l+1) f_{fd}(\epsilon_{nl},\beta,\mu) +
N_{ub}(\beta,\mu) - N_e = 0.
The number of unbound electrons :math:`N_{ub}` depends on the implementation choice.
"""
mu = config.mu
mu0 = mu # set initial guess to existing value of chem pot
# so far only the ideal treatment for unbound electrons is implemented
for i in range(config.spindims):
# increase the bracket size until a solution is found
bracket = np.array([-10, 10])
while True:
try:
x0 = mu0[i]
args = (orbs.eigvals[:, i], orbs.occ_weight[:, i], config.nele[i])
maxiter = 1000
if config.unbound == "ideal":
f_root = f_root_id
elif config.unbound == "quantum":
f_root = f_root_qu
if config.nele[i] != 0:
soln = optimize.root_scalar(
f_root,
x0=x0,
args=args,
method="brentq",
bracket=bracket,
options={"maxiter": maxiter},
)
mu[i] = soln.root
# in case there are no electrons in one spin channel
else:
mu[i] = -np.inf
break
except ValueError:
bracket *= 10
return mu
[docs]
def f_root_id(mu, eigvals, lbound, nele):
r"""
Functional input for the chemical potential root finding function (ideal approx).
See notes for function returned, the ideal approximation is used for free electrons.
Parameters
----------
mu : array_like
chemical potential
eigvals : ndarray
the energy eigenvalues
lbound : ndarray
the lbound matrix :math:`(2l+1)\Theta(\epsilon_{nl}^\sigma)`
nele : union(int, float)
the number of electrons for given spin
Returns
-------
f_root : float
the difference of the predicted electron number with given mu
and the actual electron number
Notes
-----
The returned function is
.. math:: f = \sum_{nl} (2l+1) f_{fd}(\epsilon_{nl},\beta,\mu) +
N_{ub}(\beta,\mu) - N_e
"""
# caluclate the contribution from the bound electrons
occnums = lbound * fermi_dirac(eigvals, mu, config.beta)
contrib_bound = occnums.sum()
# now compute the contribution from the unbound electrons
# this function uses the ideal approximation
prefac = (2.0 / config.spindims) * config.sph_vol / (sqrt(2) * pi**2)
contrib_unbound = prefac * fd_int_complete(mu, config.beta, 1.0)
# return the function whose roots are to be found
f_root = contrib_bound + contrib_unbound - nele
return f_root
[docs]
def f_root_qu(mu, eigvals, occ_weight, nele):
r"""
Functional input for the chemical potential root finding function (ideal approx).
See notes for function returned, the ideal approximation is used for free electrons.
Parameters
----------
mu : array_like
chemical potential
eigvals : ndarray
the energy eigenvalues
lbound : ndarray
the lbound matrix :math:`(2l+1)\Theta(\epsilon_{nl}^\sigma)`
nele : union(int, float)
the number of electrons for given spin
Returns
-------
f_root : float
the difference of the predicted electron number with given mu
and the actual electron number
Notes
-----
The returned function is
.. math:: f = \sum_{nl} (2l+1) f_{fd}(\epsilon_{nl},\beta,\mu) +
N_{ub}(\beta,\mu) - N_e
"""
# caluclate the contribution from the bound electrons
occnums = occ_weight * fermi_dirac(eigvals, mu, config.beta)
# return the function whose roots are to be found
f_root = occnums.sum() - nele
return f_root
[docs]
def lorentzian(x, x0, gamma):
r"""
Compute the Lorentzian function.
Parameters
----------
x : array_like
the "x" variable
x0 : array_like
the position of the peak
gamma : float
half the peak width / smoothing factor
Returns
-------
lorentzian_ : array_like
the lorentzian function
Notes
-----
The Lorentzian function is defined as
.. math::
\mathcal{L}(x, x_0, \gamma)=\frac{\gamma}{\pi}\frac{1}{\gamma^2+(x-x_0)^2}
"""
lorentzian_ = (gamma / np.pi) * (1.0 / (gamma**2 + (x - x0) ** 2))
return lorentzian_
[docs]
def grad_func(den, xgrid, grid_type):
"""
Compute the gradient of a function on the log / sqrt grid.
Parameters
----------
den : ndarray
density array or any other function that is integrated
xgrid : ndarray
log / sqrt grid
grid_type : str
grid type, log or sqrt
Returns
-------
grad : ndarray
The gradient of the density w.r.t. the radial grid.
"""
if grid_type == "log":
grad = (np.exp(-xgrid)) * np.gradient(den, xgrid)
else:
grad = np.gradient(den, xgrid) / (2 * xgrid)
return grad