From 616ca5abfc19a614c086e5bf9360b348d144f5d4 Mon Sep 17 00:00:00 2001 From: Prefetch Date: Sun, 23 Jul 2023 14:48:58 +0200 Subject: Initial commit for publication --- xc.py | 196 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 196 insertions(+) create mode 100644 xc.py (limited to 'xc.py') diff --git a/xc.py b/xc.py new file mode 100644 index 0000000..51fe381 --- /dev/null +++ b/xc.py @@ -0,0 +1,196 @@ +from math import * + +import numpy as np +import numpy.fft +import numpy.linalg + + +############################# +##### LDA XC potentials ##### +############################# + +# Standard X potential of homogeneous electron gas (C is neglected) +def get_x_potential_lda_slater(s, rho): + return -np.cbrt((3 / pi) * rho) + + +# XC potential used for the GTH and HGH pseudopotentials. +# See https://arxiv.org/abs/mtrl-th/9512004 +def get_xc_potential_lda_teter93(s, rho): + a0 = 0.4581652932831429 + a1 = 2.217058676663745 + a2 = 0.7405551735357053 + a3 = 0.01968227878617998 + b1 = 1.0 + b2 = 4.504130959426697 + b3 = 1.110667363742916 + b4 = 0.02359291751427506 + + rs = np.cbrt(3 / (4 * pi * rho)) + + # Energy density (epsilon) as a function of rs + eps_num = a0 + rs * (a1 + rs * (a2 + a3 * rs)) + eps_den = rs * (b1 + rs * (b2 + rs * (b3 + rs * b4))) + eps = -eps_num / eps_den + + # Derivative of energy density with respect to rs + deps_num = a0 * (b1 + rs * (2*b2 + 3*b3*rs + 4*b4*rs**2)) + deps_num += rs**3 * (2*a1*b3 + 3*a1*b4*rs - 2*a3*b1 - a3*b2*rs + a3*b4*rs**3) + deps_num += rs**2 * (a1*b2 - a2*b1 + a2*rs**2 * (b3 + 2*b4*rs)) + deps = deps_num / eps_den**2 + + # Chain rule: derivative of rs with respect to rho + chain = - 1 / np.cbrt(36 * sqrt(pi) * rho**4) + + return eps + rho * deps * chain + + +############################# +##### GGA XC potentials ##### +############################# + +# Because our grid usually has a non-orthogonal rectilinear basis, +# we first need some helper functions to calculate derivatives. + + +# Concept from tensor calculus, we'll need this later +def inverse_metric_tensor(s): + basis = s.crystal.a / s.Nres + + g = np.zeros((3, 3)) + for i in range(3): + for j in range(3): + g[i, j] = np.dot(basis[i], basis[j]) + + # Never fails if basis is linearly independent? + ginv = np.linalg.inv(g) + + return ginv + + +# Compute Laplacian of f on a grid with non-orthogonal basis vectors. +# Shamelessly taken from https://math.stackexchange.com/questions/2632598/ +def laplacian(s, f): + g = inverse_metric_tensor(s) + + lf = np.zeros_like(f) + for i in range(s.Nres): + for j in range(s.Nres): + for k in range(s.Nres): + # Wrap indices at the boundaries + inext = (i + 1) % s.Nres + jnext = (j + 1) % s.Nres + knext = (k + 1) % s.Nres + # Calculate a lot of finite differences + lf[i, j, k] += -2 * (g[0, 0] + g[1, 1] + g[2, 2]) * f[i, j, k] + lf[i, j, k] += g[0, 0] * (f[i - 1, j, k] + f[inext, j, k]) + lf[i, j, k] += g[1, 1] * (f[i, j - 1, k] + f[i, jnext, k]) + lf[i, j, k] += g[2, 2] * (f[i, j, k - 1] + f[i, j, knext]) + lf[i, j, k] += 0.5 * g[0, 1] * (f[i - 1, j - 1, k] - f[i - 1, jnext, k]) + lf[i, j, k] += 0.5 * g[0, 1] * (f[inext, jnext, k] - f[inext, j - 1, k]) + lf[i, j, k] += 0.5 * g[0, 2] * (f[i - 1, j, k - 1] - f[i - 1, j, knext]) + lf[i, j, k] += 0.5 * g[0, 2] * (f[inext, j, knext] - f[inext, j, k - 1]) + lf[i, j, k] += 0.5 * g[1, 2] * (f[i, j - 1, k - 1] - f[i, j - 1, knext]) + lf[i, j, k] += 0.5 * g[1, 2] * (f[i, jnext, knext] - f[i, jnext, k - 1]) + + return lf + + +# Compute gradient of f on a grid with non-orthogonal basis vectors +def gradient(s, f): + basis = s.crystal.a / s.Nres + g = inverse_metric_tensor(s) + + dfdx = np.zeros_like(f) + dfdy = np.zeros_like(f) + dfdz = np.zeros_like(f) + for i in range(s.Nres): + for j in range(s.Nres): + for k in range(s.Nres): + # Technically, we have periodic boundary conditions, but + # in practice we should avoid crossing the boundaries, + # hence the following verbose derivative calculations. + + # Derivative along direction indexed by i + if i == 0: + df0 = 1.0 * (f[i + 1, j, k] - f[i, j, k]) + elif i == s.Nres - 1: + df0 = 1.0 * (f[i, j, k] - f[i - 1, j, k]) + else: + df0 = 0.5 * (f[i + 1, j, k] - f[i - 1, j, k]) + # Derivative along direction indexed by j + if j == 0: + df1 = 1.0 * (f[i, j + 1, k] - f[i, j, k]) + elif j == s.Nres - 1: + df1 = 1.0 * (f[i, j, k] - f[i, j - 1, k]) + else: + df1 = 0.5 * (f[i, j + 1, k] - f[i, j - 1, k]) + # Derivative along direction indexed by k + if k == 0: + df2 = 1.0 * (f[i, j, k + 1] - f[i, j, k]) + elif k == s.Nres - 1: + df2 = 1.0 * (f[i, j, k] - f[i, j, k - 1]) + else: + df2 = 0.5 * (f[i, j, k + 1] - f[i, j, k - 1]) + + # Get the local gradient at [i, j, k] + df = [df0, df1, df2] + gradf = np.zeros(3, dtype=f.dtype) + for l in range(3): + for m in range(3): + gradf += g[l, m] * df[l] * basis[m] + + dfdx[i, j, k] = gradf[0] + dfdy[i, j, k] = gradf[1] + dfdz[i, j, k] = gradf[2] + + return dfdx, dfdy, dfdz + + +# C potential by Lee, Yang and Parr, very popular for DFT. +# See Eqs.(23)&(24) in https://doi.org/10.1103/PhysRevB.37.785 +def get_c_potential_gga_lyp(s, rho): + lr = laplacian(s, rho) + drdx, drdy, drdz = gradient(s, rho) + gr2 = drdx**2 + drdy**2 + drdz**2 + + a = 0.04918 + b = 0.132 + c = 0.2533 + d = 0.349 + CF = 0.3 * (3 * pi**2)**(2/3) + + # Abbreviations + r = rho + rr = r**(1/3) + er = np.exp(-c / rr) + + F1 = 1 / (d/rr + 1) + dF1 = d / (3*rr**2 * (d + rr)**2) + ddF1 = 2*d**2 / (9*rr**8 * (d/rr + 1)**3) - 4*d / (9*rr**7 * (d/rr + 1)**2) + + G1 = F1 * er / rr**5 + dG1 = (F1 * (c - 5*rr) + 3*dF1*rr**4) * er / (3*r**3) + ddG1 = (F1 * (c - 10*rr) * (c - 4*rr) + 6*dF1*rr**4 * (c - 5*rr) + 9*ddF1*rr**8) * er / (9*rr**13) + + V = -a * (dF1*r + F1) + V += -a*b*CF * rr**5 * (dG1*r + (8/3)*G1) + V += -(1/ 4)*a*b * ( ddG1*r*gr2 + dG1 * (3*gr2 + 2*r*lr) + 4*G1*lr) + V += -(1/72)*a*b * (3*ddG1*r*gr2 + dG1 * (5*gr2 + 6*r*lr) + 4*G1*lr) + + return V + + +# X potential by Van Leeuwen and Baerends, inspired by Becke88. +# See Eqs.(50)&(54) in https://doi.org/10.1103/PhysRevA.49.2421 +def get_x_potential_gga_lb(s, rho): + drdx, drdy, drdz = gradient(s, rho) + x = np.sqrt(drdx**2 + drdy**2 + drdz**2) / rho**(4/3) + + beta = 0.05 + fx = - beta * x**2 / (1 + 3 * beta * x * np.arcsinh(x)) + + V = get_x_potential_lda_slater(s, rho) + return V + fx * rho**(1/3) + + -- cgit v1.2.3