From b3545e3f1b604418c7360eb15b9fc46cc616f54e Mon Sep 17 00:00:00 2001 From: Deepak Date: Mon, 13 Jul 2020 15:49:15 -0500 Subject: [PATCH] Refactor baskets to use parameters --- pynam/baskets.py | 43 ++++++++++++++ pynam/dielectric/__init__.py | 3 +- pynam/dielectric/lindhard_dielectric.py | 58 +++++++++++++++++++ ...nam_dielectric_coefficient_approximator.py | 14 +++++ tests/dielectric/test_lindhard_dielectric.py | 24 ++++++++ ...nam_dielectric_coefficient_approximator.py | 23 ++++---- 6 files changed, 152 insertions(+), 13 deletions(-) create mode 100644 pynam/baskets.py create mode 100644 pynam/dielectric/lindhard_dielectric.py create mode 100644 tests/dielectric/test_lindhard_dielectric.py diff --git a/pynam/baskets.py b/pynam/baskets.py new file mode 100644 index 0000000..0cf069a --- /dev/null +++ b/pynam/baskets.py @@ -0,0 +1,43 @@ +class CalculationConstants(object): + """Holds physical constants, in SI units + """ + + def __init__(self, epsilon_0=8.854e-12, h_bar=1.0546e-34, c_light=3e8, electron_mass=9.10938356e-31): + """Initialises constants in SI units, with sensible defaults + + :param epsilon_0: + :param h_bar: + :param c_light: + :param electron_mass: + """ + + self.epsilon_0 = epsilon_0 + self.h_bar = h_bar + self.c_light = c_light + self.electron_mass = electron_mass + + +class CalculationParams(object): + """Holds the parameters describing a calculation, in SI units. + """ + + def __init__(self, omega: float = None, omega_p: float = None, tau: float = None, v_f: float = None, + t_rel: float = None, t_c: float = None, dipole_moment: float = None): + """Creates parameter object, SI units + + :param omega: + :param omega_p: + :param tau: + :param v_f: + :param t_rel: + :param t_c: + :param dipole_moment: + """ + + self.omega = omega + self.omega_p = omega_p + self.tau = tau + self.v_f = v_f + self.t_rel = t_rel + self.t_c = t_c + self.dipole_moment = dipole_moment diff --git a/pynam/dielectric/__init__.py b/pynam/dielectric/__init__.py index 23b81ed..7aa8786 100644 --- a/pynam/dielectric/__init__.py +++ b/pynam/dielectric/__init__.py @@ -1 +1,2 @@ -from pynam.dielectric.nam_dielectric_coefficient_approximator import get_nam_dielectric_coefficients +from pynam.dielectric.nam_dielectric_coefficient_approximator import get_nam_dielectric +from pynam.dielectric.lindhard_dielectric import get_lindhard_dielectric diff --git a/pynam/dielectric/lindhard_dielectric.py b/pynam/dielectric/lindhard_dielectric.py new file mode 100644 index 0000000..6ace91d --- /dev/null +++ b/pynam/dielectric/lindhard_dielectric.py @@ -0,0 +1,58 @@ +import numpy as np +from pynam.baskets import CalculationConstants, CalculationParams + +LINDHARD_SERIES_THRESHOLD = 1e4 + + +class LindhardDielectric(object): + def __init__(self, + params: CalculationParams, + constants: CalculationConstants = CalculationConstants(), + thres=LINDHARD_SERIES_THRESHOLD): + self.series_threshold = thres + self.omega = params.omega + self.v_f = params.v_f + self.omega_p = params.omega_p + self.tau = params.tau + self.c_light = constants.c_light + + self.s = 1 / (self.tau * self.omega) + self.prefactor = 3 * (self.omega_p ** 2) / (self.omega ** 2) + + def get_eps(self): + + def eps_lindhard(u_inverse_wavelength: float) -> complex: + """ the lindhard dielectric function + + :param u_inverse_wavelength: u is in units of the reciprocal vacuum wavelength (omega / c_light) + :return: returns the value of epsilon, dimensionless + """ + + # converts u from inverse vacuum wavelength to inverse mean free path + u = u_inverse_wavelength * self.v_f / self.c_light + + if u < LINDHARD_SERIES_THRESHOLD * self.c_light / self.omega: + return eps_series(u) + else: + return eps_full_lindhard(u) + + def eps_series(u: float) -> complex: + rel_value = ( + (1j / (3 * (self.s - 1j))) + u ** 2 * ((-9j + 5 * self.s) / (45 * (-1j + self.s) ** 3)) + ) + return 1 + (self.prefactor * rel_value) + + def eps_full_lindhard(u: float) -> complex: + + log_value = np.log((1 - u + (self.s * 1j)) / (1 + u + (self.s * 1j))) + + rel_value = (1 + ((1 + (self.s * 1j)) / (2 * u)) * log_value) / (1 + ((self.s * 1j) / (2 * u)) * log_value) + + return 1 + (self.prefactor / (u ** 2)) * rel_value + + return eps_lindhard + + +def get_lindhard_dielectric(params: CalculationParams, + constants: CalculationConstants = CalculationConstants()): + return LindhardDielectric(params, constants).get_eps() diff --git a/pynam/dielectric/nam_dielectric_coefficient_approximator.py b/pynam/dielectric/nam_dielectric_coefficient_approximator.py index f3a26db..2dad7a9 100644 --- a/pynam/dielectric/nam_dielectric_coefficient_approximator.py +++ b/pynam/dielectric/nam_dielectric_coefficient_approximator.py @@ -2,6 +2,8 @@ import numpy as np import pynam.dielectric.sigma_nam import pynam.dielectric.low_k_nam +from pynam.baskets import CalculationParams, CalculationConstants + from typing import Tuple FIXED_LARGE_MOMENTUM = 1e8 @@ -103,3 +105,15 @@ def get_nam_dielectric_coefficients( c, d = get_big_momentum_coefficients(dedim) return NamDielectricCoefficients(a, b, c, d) + + +def get_nam_dielectric(u_c: float, params: CalculationParams, constants: CalculationConstants = CalculationConstants()): + sigma_n = params.omega_p**2 * params.tau / (4 * np.pi) + coeffs = get_nam_dielectric_coefficients(params.omega, + sigma_n, + params.tau, + params.v_f, + params.t_rel * params.t_c, + params.t_c, + constants.c_light) + return coeffs.eps(u_c) diff --git a/tests/dielectric/test_lindhard_dielectric.py b/tests/dielectric/test_lindhard_dielectric.py new file mode 100644 index 0000000..8b6bf1d --- /dev/null +++ b/tests/dielectric/test_lindhard_dielectric.py @@ -0,0 +1,24 @@ +import pynam.dielectric +import numpy as np +import pytest +from pynam.baskets import CalculationParams + + +def get_common_lindhard_dielectric(): + params = CalculationParams(omega=1e9, omega_p=3.5e15, tau=1e-14, v_f=2e6) + return pynam.dielectric.get_lindhard_dielectric(params) + + +@pytest.mark.parametrize("test_input,expected", [ + (10, -1222.185185062794 + 1.2249999998777178e8j), + (1000, 16924.14814718176 + 1.2250000020552777e8j), + (1e8, 83.687499999706 + 0.00022417398943752126j) +]) +def test_lindhard_dielectric(test_input, expected): + + eps_to_test = get_common_lindhard_dielectric() + + np.testing.assert_almost_equal( + eps_to_test(test_input), expected, + decimal=6, err_msg='b function is off' + ) diff --git a/tests/dielectric/test_nam_dielectric_coefficient_approximator.py b/tests/dielectric/test_nam_dielectric_coefficient_approximator.py index 2a604a7..72aba86 100644 --- a/tests/dielectric/test_nam_dielectric_coefficient_approximator.py +++ b/tests/dielectric/test_nam_dielectric_coefficient_approximator.py @@ -1,6 +1,7 @@ import pytest import numpy as np import pynam.dielectric.nam_dielectric_coefficient_approximator +from pynam.baskets import CalculationParams @pytest.mark.parametrize("test_input,expected", [ @@ -14,7 +15,6 @@ import pynam.dielectric.nam_dielectric_coefficient_approximator ) ]) def test_dedimensionalise_parameters(test_input, expected): - actual_parameters = pynam.dielectric.nam_dielectric_coefficient_approximator.get_dedimensionalised_parameters( *test_input) @@ -47,12 +47,11 @@ def test_dedimensionalise_parameters(test_input, expected): # (a, b, c, d, u_l) # ) ( - (1e9, 1e16, 1e-14, 2e6, 0.8e11, 1e11, 3e8), - (3.789672906817707e10, 3.257134605133221e8, 2.655709897616547e18, 2.15e16, 7.007759408279888e7) + (1e9, 1e16, 1e-14, 2e6, 0.8e11, 1e11, 3e8), + (3.789672906817707e10, 3.257134605133221e8, 2.655709897616547e18, 2.15e16, 7.007759408279888e7) ) ]) def test_nam_coefficients(test_input, expected): - actual_coefficients = pynam.dielectric.nam_dielectric_coefficient_approximator.get_nam_dielectric_coefficients( *test_input) @@ -81,14 +80,14 @@ def test_nam_coefficients(test_input, expected): def test_nam_eps(): u_c = 1e15 - eps_to_test = pynam.dielectric.nam_dielectric_coefficient_approximator.get_nam_dielectric_coefficients( - 1e9, - 1e16, - 1e-14, - 2e6, - 0.8e11, - 1e11, - 3e8).eps(u_c) + eps_to_test = pynam.dielectric.nam_dielectric_coefficient_approximator.get_nam_dielectric(u_c, CalculationParams( + omega=1e9, + omega_p=3.54491e15, + tau=1e-14, + v_f=2e6, + t_rel=0.8, + t_c=1e11 + )) np.testing.assert_allclose( eps_to_test(10), -3.789672906817707e10 + 3.257134605133221e8j,