from collections import namedtuple
from functools import partial
import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
from scipy.special import eval_legendre
[docs]
def cubic_degeneracy(N, ny, nz):
"""Return the degeneracy of the point (N, ny, nz) on the cube.
Assumes N >= ny >= nz.
"""
assert N >= ny
assert ny >= nz
# Permutation degeneracy
if N > ny > nz:
# None equal
perm = 6
elif N == ny and ny == nz:
# All equal
perm = 1
else:
# Two equal
perm = 3
# Include sign degeneracy
degen = perm * 2 ** ((np.array([N, ny, nz]) != 0).sum())
return degen
[docs]
def get_n2s(N):
"""Return dictionary of `sum(n**2)` keys with degeneracies."""
n2s = {}
nx = N
for nz in range(N + 1):
for ny in range(nz, N + 1):
n2 = nx**2 + ny**2 + nz**2
n2s.setdefault(n2, 0)
n2s[n2] += int(cubic_degeneracy(nx, ny, nz))
return n2s
[docs]
class HomogeneousSLDA:
"""Class for computing parameters for the Homogeneous SLDA.
Attributes
----------
hbar2_2m : float
`hbar**2 / 2 / m`.
xi : float
Bertsch parameter.
alpha : float
Inverse effective mass.
eta : float
Gap in units of the Fermi energy: `eta = Delta / E_F`
aurel : bool
If `True`, then use Aurel's regularization scheme which converges faster.
"""
def __init__(self, **kw):
for key in kw:
if not hasattr(self, key):
raise ValueError(f"Unknown {key=}")
setattr(self, key, kw[key])
[docs]
def _quad(self, f, mu_alpha, delta_alpha, kc_k0):
"""Return the integral of `f(k)` over a sphere in momentum up to `kc`."""
k0 = np.sqrt(mu_alpha / self.hbar2_2m)
kc = kc_k0 * k0
integrand = lambda k: f(k, mu_alpha=mu_alpha, delta_alpha=delta_alpha) * (
k**2 / 2 / np.pi**2
)
return quad(integrand, 0, kc, points=[k0])[0]
[docs]
def _get_nk(self, k, mu_alpha, delta_alpha):
"""Carefully compute the density integrand to avoid round-off errors.
Uses a series expansion for large k.
"""
ep_alpha = self.hbar2_2m * k**2 - mu_alpha
e_alpha = np.sqrt(ep_alpha**2 + delta_alpha**2)
if np.allclose(1, ep_alpha / e_alpha):
x = (delta_alpha / ep_alpha) ** 2
# nk = 1-1/np.sqrt(1 + x)
nk = x * (1 / 2 - 3 / 8 * x + 5 / 16 * x**2 - 35 / 128 * x**3)
else:
nk = 1 - ep_alpha / e_alpha
return nk
[docs]
def _get_kF_EF_EFG(self, n):
"""Return `(k_F, E_F, E_FG)`. (DRY)
Arguments
---------
n : float or array-like
Total density n = n_a + n_b.
"""
k_F = (3 * np.pi**2 * n) ** (1 / 3)
E_F = self.hbar2_2m * k_F**2
E_FG = 3 / 5 * n * E_F
return k_F, E_F, E_FG
[docs]
def n_integrand(self, k, **kw):
return self._get_nk(k, **kw)
[docs]
def tau_integrand(self, k, **kw):
return k**2 * self.n_integrand(k, **kw)
[docs]
def nu_integrand(self, k, mu_alpha, delta_alpha):
ep_alpha = self.hbar2_2m * k**2 - mu_alpha
e_alpha = np.sqrt(ep_alpha**2 + delta_alpha**2)
return -delta_alpha / 2 / e_alpha
[docs]
def gap_eq(self, delta_alpha, **kw):
"""Return the gap equation objective (should be zero)."""
n = self._quad(self.n_integrand, delta_alpha=delta_alpha, **kw)
k_F, E_F, E_FG = self._get_kF_EF_EFG(n)
return self.alpha * delta_alpha / E_F - self.eta
[docs]
def compute(self, kc_k0=50, mu_alpha=1.0, delta_alpha=1):
"""Solve the gap equation, then compute the parameters.
Arguments
---------
kc_k0 : float
Momentum cutoff in units of k0.
mu_alpha : float
mu / alpha. Sets the scale for densities and k0.
delta_alpha : float
Initial guess for `delta / alpha`. If this is zero, then properties of the
normal state are computed.
Returns
-------
Res : tuple
Object with `alpha`, `beta` and `gamma` parameters.
Note, for debugging purposes, all variables are stored in `self._locals` which can be
used with e.g. `locals().update(s._locals)`.
"""
if delta_alpha > 0:
f = partial(self.gap_eq, mu_alpha=mu_alpha, kc_k0=kc_k0)
f0 = f1 = f(0)
while f0 * f1 > 0:
delta_alpha *= 2.0
f1 = f(delta_alpha)
delta_alpha = brentq(f, 0, delta_alpha)
n, tau, nu = [
self._quad(f, mu_alpha=mu_alpha, delta_alpha=delta_alpha, kc_k0=kc_k0)
for f in [self.n_integrand, self.tau_integrand, self.nu_integrand]
]
k_F, E_F, E_FG = self._get_kF_EF_EFG(n)
k0 = np.sqrt(mu_alpha / self.hbar2_2m)
kc = kc_k0 * k0
Lambda_c = kc / 2 / np.pi**2
if self.aurel:
Lambda_c *= 1 - k0 / 2 / kc * np.log((kc + k0) / (kc - k0))
if delta_alpha == 0:
g_alpha = nu = 0
else:
g_alpha = delta_alpha / nu
alpha = self.alpha
eta = self.eta
xi = self.xi
chi = alpha * (self.hbar2_2m * tau + delta_alpha * nu) / E_FG
beta = xi - chi
gamma = g_alpha * n ** (1 / 3) / (2 * self.hbar2_2m + g_alpha * Lambda_c)
# For comparison, we compute the homogenous infinite matter results here
y0 = mu_alpha / np.sqrt(mu_alpha**2 + delta_alpha**2)
f1_2 = self.f(l=1 / 2, y0=y0)
f3_2 = self.f(l=3 / 2, y0=y0)
f5_2 = self.f(l=5 / 2, y0=y0)
n_h = k0**3 / 4 / np.pi**2 * (y0 * f1_2 - f3_2) / y0 ** (3 / 2)
k_F_h, E_F_h, E_FG_h = self._get_kF_EF_EFG(n_h)
tau_h = k0**5 / 4 / np.pi**2 * (y0 * f3_2 - f5_2) / y0 ** (5 / 2)
nu_h = -delta_alpha / 2 * k0 / 4 / np.pi**2 / self.hbar2_2m * f1_2 / y0 ** (1 / 2)
chi_h = self.alpha * (self.hbar2_2m * tau_h + delta_alpha * nu_h) / E_FG_h
# gamma_c = -8*np.pi**2*eta / (3*np.pi**2
Res = namedtuple("Res", ["alpha", "beta", "gamma"])
res = Res(alpha=self.alpha, beta=beta, gamma=gamma)
self._locals = locals()
return res
[docs]
def f(self, y0, l):
"""Return the `-pi P_l(-y_0)/sin(pi*l)`."""
return -np.pi * eval_legendre(l, -y0) / np.sin(np.pi * l)
[docs]
def compute_lattice(self, y0, L):
pass
[docs]
def _quad_box(self, f, mu_alpha, delta_alpha):
"""Return the integral of `f(k)` over a box."""
def summand(N):
"""Return the sum over the cubic shell of "radius" `N`."""
ks = []
degen = []
for nx in range(N):
for ny in range(N):
for nz in range(N):
k = np.pi / L * np.sqrt(nx**2 + ny**2 + nz**2)
integrand = lambda k: f(k, mu_alpha=mu_alpha, delta_alpha=delta_alpha)
return quad(integrand, 0, kc, points=[k0])[0]