diff --git a/pynam/dielectric/__init__.py b/pynam/dielectric/__init__.py new file mode 100644 index 0000000..23b81ed --- /dev/null +++ b/pynam/dielectric/__init__.py @@ -0,0 +1 @@ +from pynam.dielectric.nam_dielectric_coefficient_approximator import get_nam_dielectric_coefficients diff --git a/pynam/low_k_nam.py b/pynam/dielectric/low_k_nam.py similarity index 89% rename from pynam/low_k_nam.py rename to pynam/dielectric/low_k_nam.py index bed90d0..2363303 100644 --- a/pynam/low_k_nam.py +++ b/pynam/dielectric/low_k_nam.py @@ -1,7 +1,7 @@ import numpy as np from numpy.lib.scimath import sqrt as csqrt -import pynam.complex_quad +import pynam.util.complex_quad def g(w, wp): @@ -35,7 +35,7 @@ def i2(w, wp, k, v): def a(w, k, v, t): - return pynam.complex_quad.complex_quadrature( + return pynam.util.complex_quad.complex_quadrature( lambda wp: np.tanh((w + wp) / (2 * t)) * (i1(w, wp, k, v)), 1 - w, 1 )[0] @@ -46,7 +46,7 @@ def b_int(wp, w, k, v, t): def b(w, k, v, t, b_max=np.inf): - return pynam.complex_quad.complex_quadrature( + return pynam.util.complex_quad.complex_quadrature( lambda wp: b_int(wp, w, k, v, t), 1, b_max )[0] diff --git a/pynam/dielectric/nam_dielectric_coefficient_approximator.py b/pynam/dielectric/nam_dielectric_coefficient_approximator.py new file mode 100644 index 0000000..deedb80 --- /dev/null +++ b/pynam/dielectric/nam_dielectric_coefficient_approximator.py @@ -0,0 +1,105 @@ +import numpy as np +import pynam.dielectric.sigma_nam +import pynam.dielectric.low_k_nam + +from typing import Tuple + +FIXED_LARGE_MOMENTUM = 1e8 + + +class DedimensionalisedParameters(object): + def __init__( + self, + omega: float, + sigma_n: float, + tau: float, + v_f: float, + temp: float, + critical_temp: float, + c_light: float): + gap = 0 + if temp < critical_temp: + # else, problems will happen + gap = 3.06 * np.sqrt(critical_temp * (critical_temp - temp)) + self.xi = omega / gap + self.nu = 1 / (tau * gap) + self.t = temp / gap + self.a = omega * v_f / (c_light * gap) + self.b = sigma_n / omega + + +class NamDielectricCoefficients(object): + def __init__(self, a: float, b: float, c: float, d: float): + self.a = a + self.b = b + self.c = c + self.d = d + self.u_l = np.real((-self.c + 1j * self.d) / (-self.a + 1j * self.b)) + + def eps(self): + + def piecewise_eps(u: float, u_c: float): + # todo add check for u_c vs u_l + if u < self.u_l: + return -self.a + 1j * self.b + elif self.u_l < u < u_c: + return 1 + (-self.c + 1j * self.d) / u + else: + return 1 + + return piecewise_eps + + +def get_dedimensionalised_parameters( + omega: float, + sigma_n: float, + tau: float, + v_f: float, + temp: float, + critical_temp: float, + c_light: float) -> DedimensionalisedParameters: + return DedimensionalisedParameters(omega, sigma_n, tau, v_f, temp, critical_temp, c_light) + + +def get_small_momentum_coefficients(dedim_params: DedimensionalisedParameters) -> Tuple[float, float]: + prefactor = 4j * np.pi * dedim_params.b + s = pynam.dielectric.low_k_nam.sigma_nam_alk(dedim_params.xi, 0, dedim_params.nu, dedim_params.t) + conductivity = prefactor * s + return -np.real(conductivity), np.imag(conductivity) + + +def get_big_momentum_coefficients(dedim_params: DedimensionalisedParameters) -> Tuple[float, float]: + prefactor = 4j * np.pi * dedim_params.b * FIXED_LARGE_MOMENTUM / dedim_params.a + s = pynam.dielectric.sigma_nam.sigma_nam(dedim_params.xi, + FIXED_LARGE_MOMENTUM, + dedim_params.nu, + dedim_params.t) + conductivity = prefactor * s + return -np.real(conductivity), np.imag(conductivity) + + +def get_nam_dielectric_coefficients( + omega: float, + sigma_n: float, + tau: float, + v_f: float, + temp: float, + crit_temp: float, + c_light: float) -> NamDielectricCoefficients: + """Gets a NamDielectricCoefficients object, using SI unit parameters + + :param omega: frequency + :param sigma_n: normal state conductivity + :param tau: tau in Hz + :param v_f: Fermi velocity, in m/s + :param temp: temperature in Hz + :param crit_temp: critical temperature, in Hz + :param c_light: speed of light, meters per second + :return: + """ + + dedim = get_dedimensionalised_parameters(omega, sigma_n, tau, v_f, temp, crit_temp, c_light) + a, b = get_small_momentum_coefficients(dedim) + c, d = get_big_momentum_coefficients(dedim) + + return NamDielectricCoefficients(a, b, c, d) diff --git a/pynam/sigma_nam.py b/pynam/dielectric/sigma_nam.py similarity index 88% rename from pynam/sigma_nam.py rename to pynam/dielectric/sigma_nam.py index 740f94c..ec83275 100644 --- a/pynam/sigma_nam.py +++ b/pynam/dielectric/sigma_nam.py @@ -1,7 +1,7 @@ import numpy as np from numpy.lib.scimath import sqrt as csqrt -import pynam.complex_quad +import pynam.util def g(w, wp): @@ -41,10 +41,13 @@ def i2(w, wp, k, v): def a(w, k, v, t): - return pynam.complex_quad.complex_quadrature( + result = pynam.util.complex_quad.complex_quadrature( lambda wp: np.tanh((w + wp) / (2 * t)) * (i1(w, wp, k, v)), - 1 - w, 1 - )[0] + 1 - w, 1, + epsabs=1e-10 + ) + + return result[0] def b_int(wp, w, k, v, t): @@ -52,7 +55,7 @@ def b_int(wp, w, k, v, t): def b(w, k, v, t, b_max=np.inf): - return pynam.complex_quad.complex_quadrature( + return pynam.util.complex_quadrature( lambda wp: b_int(wp, w, k, v, t), 1, b_max )[0] diff --git a/pynam/util/__init__.py b/pynam/util/__init__.py new file mode 100644 index 0000000..ef592eb --- /dev/null +++ b/pynam/util/__init__.py @@ -0,0 +1 @@ +from pynam.util.complex_quad import complex_quadrature diff --git a/pynam/complex_quad.py b/pynam/util/complex_quad.py similarity index 100% rename from pynam/complex_quad.py rename to pynam/util/complex_quad.py diff --git a/tests/dielectric/__init__.py b/tests/dielectric/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_low_k_nam.py b/tests/dielectric/test_low_k_nam.py similarity index 77% rename from tests/test_low_k_nam.py rename to tests/dielectric/test_low_k_nam.py index 38bb966..73b64db 100644 --- a/tests/test_low_k_nam.py +++ b/tests/dielectric/test_low_k_nam.py @@ -1,4 +1,4 @@ -import pynam.low_k_nam +import pynam.dielectric.low_k_nam import numpy as np import pytest @@ -9,7 +9,7 @@ import pytest ]) def test_g(test_input, expected): np.testing.assert_almost_equal( - pynam.low_k_nam.g(*test_input), expected, + pynam.dielectric.low_k_nam.g(*test_input), expected, decimal=7, err_msg='g function is off' ) @@ -21,7 +21,7 @@ def test_g(test_input, expected): ]) def test_f(test_input, expected): np.testing.assert_almost_equal( - pynam.low_k_nam.f(*test_input), expected, + pynam.dielectric.low_k_nam.f(*test_input), expected, decimal=7, err_msg='f function is off' ) @@ -32,7 +32,7 @@ def test_f(test_input, expected): ]) def test_i1(test_input, expected): np.testing.assert_almost_equal( - pynam.low_k_nam.i1(*test_input), expected, + pynam.dielectric.low_k_nam.i1(*test_input), expected, decimal=7, err_msg='i1 function is off' ) @@ -43,7 +43,7 @@ def test_i1(test_input, expected): ]) def test_i2(test_input, expected): np.testing.assert_almost_equal( - pynam.low_k_nam.i2(*test_input), expected, + pynam.dielectric.low_k_nam.i2(*test_input), expected, decimal=7, err_msg='i1 function is off' ) @@ -52,7 +52,7 @@ def test_i2(test_input, expected): ((1, 2, 3, 4), 0.228292), ]) def test_a(test_input, expected): - actual = np.real_if_close(pynam.low_k_nam.a(*test_input)) + actual = np.real_if_close(pynam.dielectric.low_k_nam.a(*test_input)) np.testing.assert_almost_equal( actual, expected, decimal=6, err_msg='a function is off' @@ -64,7 +64,7 @@ def test_a(test_input, expected): ((100, 1, 2, 3, 4), -2.62529549942e-6 + 2.60588e-12j), ]) def test_b_int(test_input, expected): - actual = np.real_if_close(pynam.low_k_nam.b_int(*test_input)) + actual = np.real_if_close(pynam.dielectric.low_k_nam.b_int(*test_input)) np.testing.assert_almost_equal( actual, expected, decimal=6, err_msg='b int function is off' @@ -76,7 +76,7 @@ def test_b_int(test_input, expected): ((1, 2, 3, 4), -0.0598057 + 0.437146j), ]) def test_b(test_input, expected): - actual = np.real_if_close(pynam.low_k_nam.b(*test_input)) + actual = np.real_if_close(pynam.dielectric.low_k_nam.b(*test_input)) np.testing.assert_almost_equal( actual, expected, decimal=6, err_msg='b function is off' @@ -88,7 +88,7 @@ def test_b(test_input, expected): ((1, 2, 3, 4), 0.98358 + 0.648221j), ]) def test_sigma_alk(test_input, expected): - actual = np.real_if_close(pynam.low_k_nam.sigma_nam_alk(*test_input)) + actual = np.real_if_close(pynam.dielectric.low_k_nam.sigma_nam_alk(*test_input)) np.testing.assert_almost_equal( actual, expected, decimal=6, err_msg='b function is off' @@ -96,7 +96,7 @@ def test_sigma_alk(test_input, expected): def test_sigma_alk_benchmark(benchmark): - result = benchmark(pynam.low_k_nam.sigma_nam_alk, 1, 2, 3, 4) + result = benchmark(pynam.dielectric.low_k_nam.sigma_nam_alk, 1, 2, 3, 4) np.testing.assert_almost_equal( result, 0.98358 + 0.648221j, 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 new file mode 100644 index 0000000..48e6f74 --- /dev/null +++ b/tests/dielectric/test_nam_dielectric_coefficient_approximator.py @@ -0,0 +1,104 @@ +import pytest +import numpy as np +import pynam.dielectric.nam_dielectric_coefficient_approximator + + +@pytest.mark.parametrize("test_input,expected", [ + # ( + # (omega, sigma_n, tau, v_f, T, T_c, c_light), + # (xi, nu, t, A, B) + # ) + ( + (1e9, 1e16, 1e-14, 2e6, 0.8e11, 1e11, 3e8), + (0.007307411691175783, 730.7411691175784, 0.5845929352940626, 0.00004871607794117188, 10000000) + ) +]) +def test_dedimensionalise_parameters(test_input, expected): + + actual_parameters = pynam.dielectric.nam_dielectric_coefficient_approximator.get_dedimensionalised_parameters(*test_input) + + np.testing.assert_almost_equal( + actual_parameters.xi, expected[0], + decimal=6, err_msg='xi incorrectly calculated' + ) + + np.testing.assert_almost_equal( + actual_parameters.nu, expected[1], + decimal=6, err_msg='nu incorrectly calculated' + ) + np.testing.assert_almost_equal( + actual_parameters.t, expected[2], + decimal=6, err_msg='t incorrectly calculated' + ) + np.testing.assert_almost_equal( + actual_parameters.a, expected[3], + decimal=6, err_msg='A incorrectly calculated' + ) + np.testing.assert_almost_equal( + actual_parameters.b, expected[4], + decimal=6, err_msg='B incorrectly calculated' + ) + + +@pytest.mark.parametrize("test_input,expected", [ + # ( + # (omega, sigma_n, tau, v_f, T, T_c, c_light), + # (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) + ) +]) +def test_nam_coefficients(test_input, expected): + + actual_coefficients = pynam.dielectric.nam_dielectric_coefficient_approximator.get_nam_dielectric_coefficients(*test_input) + + np.testing.assert_allclose( + actual_coefficients.a, expected[0], + rtol=1e-6, err_msg='a incorrectly calculated' + ) + + np.testing.assert_allclose( + actual_coefficients.b, expected[1], + rtol=1e-6, err_msg='b incorrectly calculated' + ) + np.testing.assert_allclose( + actual_coefficients.c, expected[2], + rtol=1e-2, err_msg='c incorrectly calculated' + ) + np.testing.assert_allclose( + actual_coefficients.d, expected[3], + rtol=1e-2, err_msg='d incorrectly calculated' + ) + np.testing.assert_allclose( + actual_coefficients.u_l, expected[4], + rtol=1e-5, err_msg='u_l incorrectly calculated' + ) + + +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() + + np.testing.assert_allclose( + eps_to_test(10, u_c), -3.789672906817707e10 + 3.257134605133221e8j, + rtol=1e-3, err_msg='below u_l bad' + ) + + np.testing.assert_allclose( + eps_to_test(1e10, u_c), -2.655709887616547e8 + 2.302290450767144e6j, + rtol=1e-3, err_msg='linear region bad' + ) + + np.testing.assert_allclose( + eps_to_test(1e17, u_c), 1, + rtol=1e-6, err_msg='above cutoff bad' + ) diff --git a/tests/test_sigma_nam.py b/tests/dielectric/test_sigma_nam.py similarity index 72% rename from tests/test_sigma_nam.py rename to tests/dielectric/test_sigma_nam.py index 70bb7fc..b54f743 100644 --- a/tests/test_sigma_nam.py +++ b/tests/dielectric/test_sigma_nam.py @@ -1,4 +1,4 @@ -import pynam.sigma_nam +import pynam.dielectric.sigma_nam import numpy as np import pytest @@ -9,7 +9,7 @@ import pytest ]) def test_g(test_input, expected): np.testing.assert_almost_equal( - pynam.sigma_nam.g(*test_input), expected, + pynam.dielectric.sigma_nam.g(*test_input), expected, decimal=7, err_msg='g function is off' ) @@ -21,7 +21,7 @@ def test_g(test_input, expected): ]) def test_s(test_input, expected): np.testing.assert_almost_equal( - pynam.sigma_nam.s(*test_input), expected, + pynam.dielectric.sigma_nam.s(*test_input), expected, decimal=7, err_msg='s function is off' ) @@ -36,7 +36,7 @@ def test_s(test_input, expected): ]) def test_f(test_input, expected): np.testing.assert_almost_equal( - pynam.sigma_nam.f(*test_input), expected, + pynam.dielectric.sigma_nam.f(*test_input), expected, decimal=7, err_msg='f function is off' ) @@ -48,7 +48,7 @@ def test_f(test_input, expected): ]) def test_i1(test_input, expected): np.testing.assert_almost_equal( - pynam.sigma_nam.i1(*test_input), expected, + pynam.dielectric.sigma_nam.i1(*test_input), expected, decimal=7, err_msg='i1 function is off' ) @@ -60,19 +60,20 @@ def test_i1(test_input, expected): ]) def test_i2(test_input, expected): np.testing.assert_almost_equal( - pynam.sigma_nam.i2(*test_input), expected, + pynam.dielectric.sigma_nam.i2(*test_input), expected, decimal=7, err_msg='i1 function is off' ) @pytest.mark.parametrize("test_input,expected", [ ((1, 2, 3, 4), 0.228396), + ((0.007307411691175783, 1e8, 730.7411691175784, 0.5845929352940626), 1.37272e-7) ]) def test_a(test_input, expected): - actual = np.real_if_close(pynam.sigma_nam.a(*test_input)) - np.testing.assert_almost_equal( + actual = np.real_if_close(pynam.dielectric.sigma_nam.a(*test_input)) + np.testing.assert_allclose( actual, expected, - decimal=6, err_msg='a function is off' + rtol=1e-5, err_msg='a function is off' ) @@ -81,7 +82,7 @@ def test_a(test_input, expected): ((100, 1, 2, 3, 4), -2.62529549976e-6 + 2.60765e-12j), ]) def test_b_int(test_input, expected): - actual = np.real_if_close(pynam.sigma_nam.b_int(*test_input)) + actual = np.real_if_close(pynam.dielectric.sigma_nam.b_int(*test_input)) np.testing.assert_almost_equal( actual, expected, decimal=6, err_msg='b int function is off' @@ -93,9 +94,9 @@ def test_b_int(test_input, expected): ((1, 2, 3, 4), -0.0595819 + 0.437385j), ]) def test_b(test_input, expected): - actual = np.real_if_close(pynam.sigma_nam.b(*test_input)) + actual = np.real_if_close(pynam.dielectric.sigma_nam.b(*test_input)) np.testing.assert_almost_equal( - actual, actual, + actual, expected, decimal=6, err_msg='b function is off' ) @@ -103,17 +104,18 @@ def test_b(test_input, expected): @pytest.mark.parametrize("test_input,expected", [ ((1, 2, 0, 4), 0), ((1, 2, 3, 4), 0.984117 + 0.647951j), + ((0.007307411691175783, 1e8, 730.7411691175784, 0.5845929352940626), 0.00008925294700016892 + 0.0102953966846717j) ]) def test_sigma_nam(test_input, expected): - actual = np.real_if_close(pynam.sigma_nam.sigma_nam(*test_input)) - np.testing.assert_almost_equal( + actual = np.real_if_close(pynam.dielectric.sigma_nam.sigma_nam(*test_input)) + np.testing.assert_allclose( actual, expected, - decimal=6, err_msg='sigma_nam function is off' + rtol=1e-3, err_msg='sigma_nam function is off' ) def test_sigma_nam_benchmark(benchmark): - result = benchmark(pynam.sigma_nam.sigma_nam, 1, 2, 3, 4) + result = benchmark(pynam.dielectric.sigma_nam.sigma_nam, 1, 2, 3, 4) np.testing.assert_almost_equal( result, 0.984117 + 0.647951j, decimal=6, err_msg='sigma nam benchmrak function is off' diff --git a/tests/util/__init__.py b/tests/util/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_complex_quad.py b/tests/util/test_complex_quad.py similarity index 65% rename from tests/test_complex_quad.py rename to tests/util/test_complex_quad.py index e29b479..fb69d04 100644 --- a/tests/test_complex_quad.py +++ b/tests/util/test_complex_quad.py @@ -1,9 +1,9 @@ import numpy as np -import pynam.complex_quad +import pynam.util.complex_quad def test_complex_quad(): - actual = pynam.complex_quad.complex_quadrature(lambda x: x**2 + 1j * x**3, 0, 6)[0] + actual = pynam.util.complex_quad.complex_quadrature(lambda x: x ** 2 + 1j * x ** 3, 0, 6)[0] # int_1^6 dx x^2 + i x^3 should equal (1/3)6^3 + (i/4)6^4 np.testing.assert_almost_equal( actual, (6**3)/3 + 1j*(6**4)/4,