Sarkas Yukawa MD — Plasma Transport from Molecular Dynamics

Paper: Stanton & Murillo, PRE 93, 043203 (2016) — transport coefficients
Reference: Daligault, PRE 86, 047401 (2012) — D* analytical fit
What we reproduce: Velocity-Verlet molecular dynamics of a Yukawa one-component plasma (OCP) in reduced units. We compute the velocity autocorrelation function (VACF) via Green-Kubo and extract the self-diffusion coefficient $D^*$. Small systems run live; production results (N=500+) are loaded from frozen JSON.

Live: small N=32 MD for demonstration. Production frozen data for quantitative comparison.
Rust parity: barracuda/src/md/ — GPU-accelerated MD via WGSL compute shaders.*

Physics

The Yukawa potential models screened Coulomb interactions:

$$V(r) = \frac{\Gamma}{r} e^{-\kappa r}$$

in reduced units ($a_{\text{ws}} = 1$, $\omega_p = 1$). The coupling parameter $\Gamma = q^2/(a_{\text{ws}} k_B T)$ controls the thermodynamic state:

• $\Gamma \ll 1$: weakly coupled, gas-like
• $\Gamma \sim 1-10$: moderately coupled liquid
• $\Gamma \gg 100$: strongly coupled, crystalline

The self-diffusion coefficient is obtained from the Green-Kubo relation:

$$D = \frac{1}{3} \int_0^\infty \langle \mathbf{v}(0) \cdot \mathbf{v}(t) \rangle dt$$

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import time

trapezoid = np.trapezoid if hasattr(np, 'trapezoid') else np.trapz

C_PASS = '#2ecc71'
C_FAIL = '#e74c3c'
C_INFO = '#3498db'
C_RUST = '#1abc9c'
C_PYTHON = '#f39c12'

def fcc_lattice(n_particles, box_side):
    n_cells = int(np.ceil((n_particles / 4)**(1.0/3.0)))
    basis = np.array([[0,0,0],[.5,.5,0],[.5,0,.5],[0,.5,.5]])
    positions = []
    cell_len = box_side / n_cells
    for ix in range(n_cells):
        for iy in range(n_cells):
            for iz in range(n_cells):
                for b in basis:
                    positions.append((np.array([ix,iy,iz])+b)*cell_len)
                    if len(positions) >= n_particles:
                        return np.array(positions[:n_particles])
    return np.array(positions[:n_particles])

def yukawa_forces(pos, n, box, kappa, rc):
    pe = 0.0
    forces = np.zeros_like(pos)
    rc2 = rc*rc
    for i in range(n):
        for j in range(i+1, n):
            d = pos[i] - pos[j]
            d -= box*np.round(d/box)
            r2 = np.dot(d, d)
            if r2 < rc2:
                r = np.sqrt(r2)
                exp_kr = np.exp(-kappa*r)
                pe += exp_kr/r
                f = exp_kr*(1/r2 + kappa/r)/r
                fv = d*f
                forces[i] += fv
                forces[j] -= fv
    return forces, pe

print("MD functions ready")

Small Live MD — N=32 Yukawa OCP

A small-system demonstration of velocity-Verlet MD with Yukawa forces. Too small for quantitative transport, but shows the algorithm.

N = 32
Gamma = 1.0
kappa = 1.0
dt = 0.005
n_equil = 200
n_prod = 500
rc = 5.0

box = (4*np.pi*N/3)**(1/3)
T_target = 1.5/Gamma

pos = fcc_lattice(N, box)
rng = np.random.default_rng(42)
vel = rng.normal(0, np.sqrt(T_target), (N, 3))
vel -= vel.mean(axis=0)

forces, pe = yukawa_forces(pos, N, box, kappa, rc)

t0 = time.perf_counter()
# Equilibration with Berendsen thermostat
for step in range(n_equil):
    vel += 0.5*dt*forces
    pos += dt*vel
    pos %= box
    forces, pe = yukawa_forces(pos, N, box, kappa, rc)
    vel += 0.5*dt*forces
    ke = 0.5*np.sum(vel**2)
    T_inst = 2*ke/(3*N)
    lam = np.sqrt(T_target/max(T_inst, 1e-10))
    vel *= 0.99 + 0.01*lam  # weak thermostat

# Production NVE
pe_hist, ke_hist, vacf_data = [], [], []
v0 = vel.copy()
for step in range(n_prod):
    vel += 0.5*dt*forces
    pos += dt*vel
    pos %= box
    forces, pe = yukawa_forces(pos, N, box, kappa, rc)
    vel += 0.5*dt*forces
    ke = 0.5*np.sum(vel**2)
    pe_hist.append(pe)
    ke_hist.append(ke)
    vacf = np.mean(np.sum(v0*vel, axis=1))
    vacf_data.append(vacf)

elapsed = time.perf_counter() - t0
print(f"MD: N={N}, Gamma={Gamma}, kappa={kappa}")
print(f"  {n_equil} equil + {n_prod} prod steps in {elapsed*1000:.0f} ms")
print(f"  <T> = {2*np.mean(ke_hist)/(3*N):.4f} (target {T_target:.4f})")

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Energy conservation
total_e = np.array(pe_hist) + np.array(ke_hist)
axes[0].plot(pe_hist, label='PE', color=C_INFO, alpha=0.7)
axes[0].plot(ke_hist, label='KE', color=C_FAIL, alpha=0.7)
axes[0].plot(total_e, label='Total', color=C_PASS, linewidth=2)
axes[0].set_xlabel('Step')
axes[0].set_ylabel('Energy (reduced)')
axes[0].set_title(f'Energy Conservation (N={N})')
axes[0].legend(fontsize=9)

# VACF
t_arr = np.arange(len(vacf_data))*dt
vacf_norm = np.array(vacf_data)/max(vacf_data[0], 1e-10)
axes[1].plot(t_arr, vacf_norm, color=C_INFO, linewidth=2)
axes[1].axhline(y=0, color='gray', ls='--', alpha=0.5)
axes[1].set_xlabel('Time ($\\omega_p^{-1}$)')
axes[1].set_ylabel('VACF / VACF(0)')
axes[1].set_title('Velocity Autocorrelation')

# D* from Green-Kubo
D_gk = trapezoid(vacf_data, t_arr)/3.0
axes[2].text(0.5, 0.7, f'$D^* \\approx {D_gk:.4f}$', transform=axes[2].transAxes,
             fontsize=16, ha='center', va='center',
             bbox=dict(boxstyle='round', facecolor=C_INFO, alpha=0.3))
axes[2].text(0.5, 0.4, f'$\\Gamma = {Gamma}$, $\\kappa = {kappa}$',
             transform=axes[2].transAxes, fontsize=14, ha='center')
axes[2].text(0.5, 0.2, f'N={N} (small-system demo)',
             transform=axes[2].transAxes, fontsize=10, ha='center', color='gray')
axes[2].set_title('Self-Diffusion Coefficient')
axes[2].axis('off')

fig.suptitle('Yukawa OCP Molecular Dynamics', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('/tmp/hotspring_paper_03_md.png', dpi=150, bbox_inches='tight')
plt.show()

Daligault Analytical Fit — D*(Gamma, kappa)

The Daligault (2012) model interpolates between weak-coupling (Landau-Spitzer) and strong-coupling (Einstein) limits for the self-diffusion coefficient.

def d_star_daligault(gamma, kappa):
    gamma_eff = gamma * np.exp(-kappa)
    cl = max(np.log(1+1/max(gamma_eff, 0.01)), 0.1)
    dw = 3*np.sqrt(np.pi)/4 / (gamma**2.5 * cl)
    a = 0.0094 + 0.018*kappa - 0.0025*kappa**2
    alpha = 1.09 + 0.12*kappa - 0.019*kappa**2
    ds = a * gamma**(-alpha)
    gamma_x = 10*np.exp(0.5*kappa)
    f = 1/(1+(gamma/gamma_x)**2)
    return dw*f + ds*(1-f)

gammas = np.logspace(-1, 3, 200)
kappas = [0.0, 1.0, 2.0, 3.0]

fig, ax = plt.subplots(figsize=(10, 6))
colors_k = [C_INFO, C_PASS, C_PYTHON, C_FAIL]
for kap, color in zip(kappas, colors_k):
    d_vals = [d_star_daligault(g, kap) for g in gammas]
    ax.loglog(gammas, d_vals, color=color, linewidth=2, label=f'$\\kappa={kap}$')

ax.set_xlabel('$\\Gamma$')
ax.set_ylabel('$D^*$')
ax.set_title('Daligault (2012) — Self-Diffusion Coefficient')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('/tmp/hotspring_paper_03_daligault.png', dpi=150, bbox_inches='tight')
plt.show()

Rust Parity

The full MD engine is in barracuda/src/md/ with GPU-accelerated force computation via WGSL compute shaders. Production runs (N=500-2000) demonstrate quantitative agreement with published transport coefficients.

Provenance

• Papers: Stanton & Murillo PRE 93 (2016), Daligault PRE 86 (2012)
• Control: control/sarkas/simulations/transport-study/scripts/
• Validation: barracuda/src/md/, validate_yukawa_md
• Next: Production sweeps via GPU primal composition, Sarkas comparison

hotSpring — ecoPrimals · AGPL-3.0 · primals.eco