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)