diff --git a/nonrad/scaling.py b/nonrad/scaling.py index 90bd6a8..bb34cc0 100644 --- a/nonrad/scaling.py +++ b/nonrad/scaling.py @@ -10,6 +10,7 @@ from typing import Any, Optional, Union import numpy as np +from mpmath import whitw from numpy.polynomial.laguerre import laggauss from pymatgen.io.vasp.outputs import Wavecar from scipy import constants as const @@ -25,14 +26,81 @@ def _njit(func): return _njit +def _s_k( + k: Union[float, np.ndarray], + Z: int, + m_eff: float, + eps0: float, + dim: int = 3, + x0: float = 1e-3, +) -> Union[float, np.ndarray]: + """Compute the Sommerfeld parameter as a function of momentum. + + Parameters + ---------- + k : float, np.array(dtype=float) + wavevector + Z : int + Z = Q / q where Q is the charge of the defect and q is the charge of + the carrier. Z < 0 corresponds to attractive centers and Z > 0 + corresponds to repulsive centers + m_eff : float + effective mass of the carrier in units of m_e (electron mass) + eps0 : float + static dielectric constant + dim : int + dimension of the system (3, 2, or 1) + x0 : float + cutoff length for 1D evaluation + method : str + specify method for evaluating sommerfeld parameter ('Integrate' or + 'Analytic'). The default is recommended as the analytic equation may + introduce significant errors for the repulsive case at high T. + + Returns + ------- + float, np.array(dtype=float) + sommerfeld factor evaluated at the given temperature + """ + if Z == 0: + return 0. + + a = (eps0 / m_eff) * const.value('Bohr radius') + nu = Z / a / k + if dim == 1: + def _phi(Z, k, a, x, x0): + nu = np.abs(Z) / a / k + z0 = 2*1j*k*x0 + z = 2*1j*k*np.abs(x) + z0 + + W1 = lambda nu, z: np.exp(z0/2) * whitw(-1j*nu, 1/2, z) # noqa: E731 + W2 = lambda nu, z: np.exp(-z0/2) * whitw(1j*nu, 1/2, -z) # noqa: E731 + + D10 = (W1(nu, 1.01*z0) - W1(nu, 0.99*z0))/(2*0.01*z0) + D20 = (W2(nu, 1.01*z0) - W2(nu, 0.99*z0))/(2*0.01*z0) + + N = np.sqrt(np.exp(-np.pi * nu) / 2 / + (np.abs(D10)**2 + np.abs(D20)**2)) + + return N*(D20*W1(nu, z) - D10*W2(nu, z)) + + return np.vectorize(lambda q: np.abs(_phi(Z, q, a, 0., x0*a))**2)(k) + elif dim == 2: + return 2 / (1 + np.exp(2 * np.pi * nu)) + else: + return 2 * np.pi * nu / (np.exp(2 * np.pi * nu) - 1) + + def sommerfeld_parameter( T: Union[float, np.ndarray], Z: int, m_eff: float, eps0: float, + dim: int = 3, + x0: float = 1e-3, method: str = 'Integrate' ) -> Union[float, np.ndarray]: - """Compute the sommerfeld parameter. + """Compute the T-dependent Sommerfeld parameter. Computes the sommerfeld parameter at a given temperature using the definitions in R. Pässler et al., phys. stat. sol. (b) 78, 625 (1976). We @@ -50,6 +118,10 @@ def sommerfeld_parameter( effective mass of the carrier in units of m_e (electron mass) eps0 : float static dielectric constant + dim : int + dimension of the system (3, 2, or 1) + x0 : float + cutoff length for 1D evaluation method : str specify method for evaluating sommerfeld parameter ('Integrate' or 'Analytic'). The default is recommended as the analytic equation may @@ -64,29 +136,47 @@ def sommerfeld_parameter( return 1. if method.lower()[0] == 'i': - kT = const.k * T - m = m_eff * const.m_e - eps = (4 * np.pi * const.epsilon_0) * eps0 - f = -2 * np.pi * Z * m * const.e**2 / const.hbar**2 / eps - - def s_k(k): - return f / k / (1 - np.exp(-f / k)) + mkT = m_eff * const.m_e * const.k * T + x_to_k = np.sqrt(2*mkT)/const.hbar + + def _f(k, Z, m_eff, eps0, dim, x0): + if dim == 1: + return _s_k(k, Z, m_eff, eps0, dim=dim, x0=x0) / k + elif dim == 2: + return _s_k(k, Z, m_eff, eps0, dim=dim) + else: + return k * _s_k(k, Z, m_eff, eps0, dim=dim) + + def _norm(mkT, dim): + if dim == 1: + return np.sqrt(np.pi*mkT/2) / const.hbar + elif dim == 2: + return mkT / const.hbar**2 + else: + return np.sqrt(np.pi/2) * (mkT)**(3/2) / const.hbar**3 t = 0. x, w = laggauss(64) for ix, iw in zip(x, w): - t += iw * np.sqrt(ix) * s_k(np.sqrt(2 * m * kT * ix) / const.hbar) - return t / np.sum(w * np.sqrt(x)) + t += iw * _f(x_to_k*np.sqrt(ix), Z, m_eff, eps0, dim, x0) + t *= (mkT/const.hbar**2) + return t / _norm(mkT, dim) - # that 4*pi from Gaussian units.... theta_b = np.pi**2 * (m_eff * const.m_e) * const.e**4 / \ (2 * const.k * const.hbar**2 * (eps0 * 4*np.pi*const.epsilon_0)**2) zthetaT = Z**2 * theta_b / T - - if Z < 0: - return 4 * np.sqrt(zthetaT / np.pi) - return (8 / np.sqrt(3)) * \ - zthetaT**(2/3) * np.exp(-3 * zthetaT**(1/3)) + if dim == 1: + raise ValueError('Cannot analytically evaluate Sommerfeld parameter in 1D') + elif dim == 2: + if Z < 0: + return 2. + return np.sqrt(8*np.pi/3) * \ + (8*zthetaT)**(1/6) * np.exp(-3 * zthetaT**(1/3)) + else: + if Z < 0: + return 4 * np.sqrt(zthetaT / np.pi) + return (8 / np.sqrt(3)) * \ + zthetaT**(2/3) * np.exp(-3 * zthetaT**(1/3)) @njit(cache=True) diff --git a/nonrad/tests/test_scaling.py b/nonrad/tests/test_scaling.py index 607c1e5..8e6bb8c 100644 --- a/nonrad/tests/test_scaling.py +++ b/nonrad/tests/test_scaling.py @@ -1,9 +1,11 @@ # pylint: disable=C0114,C0115,C0116 import unittest +from typing import Union import numpy as np from scipy import constants as const +from numpy.polynomial.laguerre import laggauss from nonrad.scaling import ( charged_supercell_scaling, @@ -17,6 +19,42 @@ from nonrad.tests import TEST_FILES, FakeFig +def _old_sommerfeld_parameter( + T: Union[float, np.ndarray], + Z: int, + m_eff: float, + eps0: float, + method: str = 'Integrate' +) -> Union[float, np.ndarray]: + if Z == 0: + return 1. + + if method.lower()[0] == 'i': + kT = const.k * T + m = m_eff * const.m_e + eps = (4 * np.pi * const.epsilon_0) * eps0 + f = -2 * np.pi * Z * m * const.e**2 / const.hbar**2 / eps + + def s_k(k): + return f / k / (1 - np.exp(-f / k)) + + t = 0. + x, w = laggauss(64) + for ix, iw in zip(x, w): + t += iw * np.sqrt(ix) * s_k(np.sqrt(2 * m * kT * ix) / const.hbar) + return t / np.sum(w * np.sqrt(x)) + + # that 4*pi from Gaussian units.... + theta_b = np.pi**2 * (m_eff * const.m_e) * const.e**4 / \ + (2 * const.k * const.hbar**2 * (eps0 * 4*np.pi*const.epsilon_0)**2) + zthetaT = Z**2 * theta_b / T + + if Z < 0: + return 4 * np.sqrt(zthetaT / np.pi) + return (8 / np.sqrt(3)) * \ + zthetaT**(2/3) * np.exp(-3 * zthetaT**(1/3)) + + class SommerfeldTest(unittest.TestCase): def setUp(self): self.args = { @@ -61,7 +99,7 @@ def test_list(self): def test_compare_methods(self): self.args = { - 'T': 150, + 'T': 100, 'Z': -1, 'm_eff': 0.2, 'eps0': 8.9, @@ -80,6 +118,48 @@ def test_compare_methods(self): f1 = sommerfeld_parameter(**self.args) self.assertGreater(np.abs(f0-f1)/f1, 0.1) + def test_old_sommerfeld(self): + self.args = {'m_eff': 0.2, 'eps0': 8.9} + for m in ['i', 'a']: + self.args['method'] = m + for t in [100, 300, 700, 900]: + self.args['T'] = t + for z in [0, 1, -1]: + self.args['Z'] = z + f0 = _old_sommerfeld_parameter(**self.args) + f1 = sommerfeld_parameter(**self.args) + self.assertAlmostEqual(f0, f1, places=2) + + def test_sommerfeld_dim(self): + self.args = { + 'T': 200, + 'Z': -1, + 'm_eff': 0.2, + 'eps0': 8.9, + 'dim': 2, + 'method': 'Integrate' + } + + self.assertAlmostEqual(sommerfeld_parameter(**self.args), 2., places=2) + self.args['method'] = 'a' + self.assertAlmostEqual(sommerfeld_parameter(**self.args), 2., places=5) + self.args['method'] = 'i' + self.args['Z'] = 1 + self.assertLess(sommerfeld_parameter(**self.args), 1.) + self.args['method'] = 'a' + self.assertLess(sommerfeld_parameter(**self.args), 1.) + + self.args['dim'] = 1 + with self.assertRaises(ValueError): + self.assertLess(sommerfeld_parameter(**self.args), 1.) + self.args['method'] = 'i' + self.assertLess(sommerfeld_parameter(**self.args), 1.) + self.args['Z'] = -1 + self.assertLess(sommerfeld_parameter(**self.args), 1.) + + self.args['Z'] = 0 + self.assertEqual(sommerfeld_parameter(**self.args), 1.) + class ChargedSupercellScalingTest(unittest.TestCase): def test_find_charge_center(self):