Εφαρμογή DFT: Θεμελιώδης Κατάσταση Ατόμου Be

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve_banded, eigh_tridiagonal

# ─────────────────────────────────────────────
# Radial grid
# ─────────────────────────────────────────────
rmax = 65.0
dr   = 0.015
r    = np.arange(dr/2, rmax + dr/2, dr)
N    = len(r)
vol  = 4 * np.pi * r**2

# ─────────────────────────────────────────────
# Atomic number: Beryllium Z=4, 4 electrons
# occupying 1s² and 2s²
# ─────────────────────────────────────────────
Z        = 4.0
N_elec   = 4      # total number of electrons
N_orb    = 2      # number of orbitals: 1s and 2s
occ      = [2, 2] # occupancy of each orbital

# ─────────────────────────────────────────────
# Kinetic energy (Laplacian) – same as before
# ─────────────────────────────────────────────
laplacian_diag = -2.0 * np.ones(N)     / dr**2
laplacian_off  =  1.0 * np.ones(N - 1) / dr**2

# ─────────────────────────────────────────────
# Poisson solver matrix – same as before
# ─────────────────────────────────────────────
poisson_diag   =  2.0 * np.ones(N)     / dr**2
poisson_off    = -1.0 * np.ones(N - 1) / dr**2
poisson_banded = np.vstack([
    np.hstack((0, poisson_off)),
    poisson_diag,
    np.hstack((poisson_off, 0))
])

# ─────────────────────────────────────────────
# Initial electron density for Be (4 electrons)
# ─────────────────────────────────────────────
n = (N_elec * Z**3 / np.pi) * np.exp(-2 * Z * r)

# ─────────────────────────────────────────────
# LDA (Perdew-Zunger 1981) parameters
# ─────────────────────────────────────────────
A      =  0.0311
B      = -0.048
C      =  0.002
D      = -0.0116
gamma  = -0.1423
beta1  =  1.0529
beta2  =  0.3334

def get_Vxc(n):
    """LDA exchange-correlation potential."""
    rs = (3.0 / (4.0 * np.pi * n))**(1.0/3.0)
    Vx = -(3.0 / np.pi * n)**(1.0/3.0)
    Vc = np.where(
        rs < 1,
        A*np.log(rs) + B + C*rs*np.log(rs) + D*rs
            - rs/3.0 * (A/rs + C*np.log(rs) + C + D),
        gamma*(1 + 7.0/6.0*beta1*np.sqrt(rs) + 4.0/3.0*beta2*rs)
            / (1 + beta1*np.sqrt(rs) + beta2*rs)**2
    )
    return Vx + Vc

def get_eps_xc(n):
    """LDA exchange-correlation energy density."""
    rs_ec = (3.0 / (4.0 * np.pi * n))**(1.0/3.0)
    eps_x = -0.75 * (3.0 / np.pi * n)**(1.0/3.0)
    eps_c = np.where(
        rs_ec < 1,
        A*np.log(rs_ec) + B + C*rs_ec*np.log(rs_ec) + D*rs_ec,
        gamma / (1 + beta1*np.sqrt(rs_ec) + beta2*rs_ec)
    )
    return eps_x, eps_c

# ─────────────────────────────────────────────
# SCF loop
# ─────────────────────────────────────────────
alpha    = 0.3
tol      = 1e-5
max_iter = 300

energy_history = []
E_ref_Be = -14.5730  # Reference LDA energy for Be (Ha)

print("SCF loop for Beryllium (Z=4, 1s² 2s²)\n")

for iteration in range(max_iter):

    # Hartree potential via Poisson equation
    rhs      = 4 * np.pi * r * n
    rhs[0]   = 0.0
    rhs[-1]  = N_elec   # boundary: U(rmax) → N_elec
    U        = solve_banded((1, 1), poisson_banded, rhs)
    VH       = U / r

    # XC and effective potential
    Vxc  = get_Vxc(n)
    Vext = -Z / r
    Veff = Vext + VH + Vxc

    # Solve KS equation – get lowest N_orb orbitals
    diag = -0.5 * laplacian_diag + Veff
    off  = -0.5 * laplacian_off
    eps, eigvec = eigh_tridiagonal(
        diag, off,
        select='i', select_range=(0, N_orb - 1)
    )

    # Normalize orbitals and build new density
    orbitals = []
    n_new    = np.zeros(N)
    for i in range(N_orb):
        u_i  = eigvec[:, i]
        u_i /= np.sqrt(np.trapz(u_i**2, r))
        orbitals.append(u_i)
        n_new += occ[i] * u_i**2 / (4 * np.pi * r**2)

    # Energies
    eps_x_vals, eps_c_vals = get_eps_xc(n)
    Ts   = sum(occ[i]*eps[i] for i in range(N_orb)) - np.trapz(Veff * n * vol, r)
    EH   = 0.5 * np.trapz(VH          * n * vol, r)
    Ex   =       np.trapz(eps_x_vals  * n * vol, r)
    Ec   =       np.trapz(eps_c_vals  * n * vol, r)
    Eext =       np.trapz(Vext        * n * vol, r)
    Etot = Ts + EH + Ex + Ec + Eext
    energy_history.append(Etot)

    # Convergence check
    delta_n = np.max(np.abs(n_new - n))
    print(f"Iter {iteration+1:4d} | Etot = {Etot:.6f} Ha | "
          f"ε_1s = {eps[0]:.4f} | ε_2s = {eps[1]:.4f} | Δn = {delta_n:.2e}")

    # Density mixing
    n = (1 - alpha) * n + alpha * n_new

    if delta_n < tol:
        print(f"\nConverged in {iteration+1} iterations")
        print(f"  Final total energy : Etot  = {Etot:.6f} Ha")
        print(f"  Reference (LDA)    : E_Be  = {E_ref_Be:.6f} Ha")
        print(f"  Difference         : ΔE    = {abs(Etot - E_ref_Be):.6f} Ha")
        print(f"  1s orbital energy  : ε_1s  = {eps[0]:.6f} Ha")
        print(f"  2s orbital energy  : ε_2s  = {eps[1]:.6f} Ha")
        break
else:
    print("Did not converge within max iterations!")

# ─────────────────────────────────────────────
# Plots
# ─────────────────────────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Plot 1: KS orbitals
u_1s = orbitals[0]
u_2s = orbitals[1]
axes[0].plot(r, u_1s, color='steelblue', linewidth=2, label='1s')
axes[0].plot(r, u_2s, color='darkorange', linewidth=2, label='2s')
axes[0].set_xlabel('r [Bohr]')
axes[0].set_ylabel('u(r) [a.u.]')
axes[0].set_title('KS Orbitals u(r)')
axes[0].set_xlim(0, 10)
axes[0].axhline(0, color='black', linewidth=0.5, linestyle='--')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Plot 2: Electron density
axes[1].plot(r, vol * n, color='tomato', linewidth=2)
axes[1].set_xlabel('r [Bohr]')
axes[1].set_ylabel('4πr²n(r) [a.u.]')
axes[1].set_title('Electron Density n(r)')
axes[1].set_xlim(0, 10)
axes[1].grid(True, alpha=0.3)

# Plot 3: Energy convergence
axes[2].plot(range(1, len(energy_history)+1), energy_history,
             color='purple', linewidth=2, marker='o', markersize=3)
axes[2].axhline(y=E_ref_Be, color='red', linestyle='--',
                linewidth=1.5, label=f'Reference: {E_ref_Be} Ha')
axes[2].set_xlabel('SCF Iteration')
axes[2].set_ylabel('Total Energy [Ha]')
axes[2].set_title('Total Energy Convergence')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.suptitle('Beryllium Atom (Z=4) DFT Results',
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

# ─────────────────────────────────────────────
# Diagnostics
# ─────────────────────────────────────────────
N_e_num = np.trapz(vol * n, r)
print(f"\nElectron Count Check: N_e = {N_e_num:.6f} (should be {N_elec})")
for i, u_i in enumerate(orbitals):
    norm = np.trapz(u_i**2, r)
    print(f"Normalization orbital {i+1}: {norm:.10f} (should be 1)")