93 lines
2.5 KiB
Python
93 lines
2.5 KiB
Python
"""
|
||
LDA exchange-correlation functional.
|
||
|
||
Exchange: Dirac/Slater
|
||
ε_x(ρ) = −(3/4)(3/π)^(1/3) ρ^(1/3)
|
||
v_x(ρ) = (4/3) ε_x(ρ)
|
||
|
||
Correlation: Perdew-Zunger (1981) parametrization of Ceperley-Alder QMC data.
|
||
ε_c and v_c as a function of r_s = (3/(4π ρ))^(1/3).
|
||
|
||
Reference: Perdew & Zunger, Phys. Rev. B 23, 5048 (1981).
|
||
"""
|
||
|
||
import numpy as np
|
||
|
||
# PZ81 correlation constants (spin-unpolarized)
|
||
_PZ = {
|
||
'gamma': -0.1423,
|
||
'beta1': 1.0529,
|
||
'beta2': 0.3334,
|
||
'A': 0.0311,
|
||
'B': -0.0480,
|
||
'C': 0.0020,
|
||
'D': -0.0116,
|
||
}
|
||
|
||
|
||
def lda_xc(rho: np.ndarray):
|
||
"""
|
||
LDA exchange-correlation energy density and potential.
|
||
|
||
Parameters
|
||
----------
|
||
rho : (N,) array Electron density (e/Bohr³), must be non-negative.
|
||
|
||
Returns
|
||
-------
|
||
exc : (N,) array XC energy density ε_xc(ρ) [Ha].
|
||
vxc : (N,) array XC potential v_xc(ρ) = δE_xc/δρ [Ha].
|
||
"""
|
||
rho = np.maximum(rho, 1e-20) # avoid log/division singularities
|
||
|
||
# ---------- Exchange ----------
|
||
rho_cbrt = rho ** (1.0 / 3.0)
|
||
CX = -(3.0 / 4.0) * (3.0 / np.pi) ** (1.0 / 3.0)
|
||
ex = CX * rho_cbrt
|
||
vx = (4.0 / 3.0) * ex # v_x = dε_x/dρ · ρ + ε_x ... but v_x = (4/3)ε_x
|
||
|
||
# ---------- Correlation (PZ81) ----------
|
||
rs = (3.0 / (4.0 * np.pi * rho)) ** (1.0 / 3.0)
|
||
|
||
ec = np.empty_like(rs)
|
||
vc = np.empty_like(rs)
|
||
|
||
hi = rs >= 1.0 # low-density regime
|
||
lo = ~hi # high-density regime
|
||
|
||
# --- High-density (r_s < 1) ---
|
||
if lo.any():
|
||
rs_lo = rs[lo]
|
||
A, B, C, D = _PZ['A'], _PZ['B'], _PZ['C'], _PZ['D']
|
||
lnrs = np.log(rs_lo)
|
||
ec[lo] = A * lnrs + B + C * rs_lo * lnrs + D * rs_lo
|
||
# v_c = ε_c − (r_s/3) dε_c/dr_s
|
||
dec_drs = A / rs_lo + C * (lnrs + 1) + D
|
||
vc[lo] = ec[lo] - (rs_lo / 3.0) * dec_drs
|
||
|
||
# --- Low-density (r_s ≥ 1) ---
|
||
if hi.any():
|
||
rs_hi = rs[hi]
|
||
gamma, beta1, beta2 = _PZ['gamma'], _PZ['beta1'], _PZ['beta2']
|
||
sqrt_rs = np.sqrt(rs_hi)
|
||
denom = 1.0 + beta1 * sqrt_rs + beta2 * rs_hi
|
||
ec[hi] = gamma / denom
|
||
dec_drs = -gamma * (0.5 * beta1 / sqrt_rs + beta2) / denom**2
|
||
vc[hi] = ec[hi] - (rs_hi / 3.0) * dec_drs
|
||
|
||
exc = ex + ec
|
||
vxc = vx + vc
|
||
return exc, vxc
|
||
|
||
|
||
def xc_energy(rho: np.ndarray, dV: float) -> float:
|
||
"""Integrate XC energy: E_xc = ∫ ε_xc(ρ) ρ dr."""
|
||
exc, _ = lda_xc(rho)
|
||
return float(np.sum(exc * rho) * dV)
|
||
|
||
|
||
def xc_potential(rho: np.ndarray) -> np.ndarray:
|
||
"""Return XC potential v_xc(r) on the grid."""
|
||
_, vxc = lda_xc(rho)
|
||
return vxc
|