#!/usr/bin/env python3 """Wetting of trimethylsilylated amorphous silica from interfacial thermodynamics. Self-contained driver. Recomputes every number reported in SR-2026-001. No contact-angle measurement is used as input anywhere in this file. python3 driver.py # writes results.json """ import json, numpy as np from scipy.spatial import cKDTree from scipy.optimize import brentq, curve_fit # ------------------------------------------------------------------ constants kB, hbar, NA = 1.380649e-23, 1.054571817e-34, 6.02214076e23 T = 298.15 kT = kB * T xi1 = 2 * np.pi * kT / hbar # first Matsubara frequency, rad/s D0 = 1.65e-10 # universal contact separation, m GAM_W = 72.8e-3 # water surface tension, J/m^2 (L) # molecular geometry (bond lengths, vdW radii) -- all (L) d_SiC, d_SiO = 1.87, 1.63 # Angstrom tet = np.deg2rad(109.47) r_CH3, r_Ow, r_Si = 2.00, 1.52, 2.10 # Angstrom d_HB = 2.85 # H-bond O...O, Angstrom h_Si = 1.50 # grafted Si above the anchor plane SIG_OH0 = 0.046 # fully hydroxylated silica, /Angstrom^2 N_HB_PER_OH = 1.5 # donor + acceptor on an isolated silanol (A) N_HB_SILOX = 1.0 # siloxane bridge: acceptor only (A) Z_HB, L_HB = 3.6, 2.83e-10 # liquid-water H-bond count and length (L) TRAPZ = np.trapezoid if hasattr(np, "trapezoid") else np.trapz # ------------------------------------------- dielectric response, imaginary axis MAT = { "water": dict(e0=78.4, eIR=4.35, n=1.333, wMW=1.20e11, wIR=1.00e14, wUV=1.905e16), "tms": dict(e0=2.05, eIR=2.05, n=1.400, wMW=None, wIR=2.0e14, wUV=1.87e16), "alkane": dict(e0=2.02, eIR=2.02, n=1.421, wMW=None, wIR=2.0e14, wUV=1.87e16), "silica": dict(e0=3.82, eIR=3.82, n=1.448, wMW=None, wIR=2.07e14, wUV=2.03e16), "vac": dict(e0=1.0, eIR=1.0, n=1.000, wMW=None, wIR=1.0e14, wUV=1.0e16), } def eps(xi, p): """eps(i xi): one Debye (microwave) + one IR + one UV Lorentz oscillator.""" out = np.ones_like(np.atleast_1d(xi), dtype=float) if p["wMW"] is not None and p["e0"] > p["eIR"]: out = out + (p["e0"] - p["eIR"]) / (1.0 + xi / p["wMW"]) out = out + (p["eIR"] - p["n"] ** 2) / (1.0 + (xi / p["wIR"]) ** 2) return out + (p["n"] ** 2 - 1.0) / (1.0 + (xi / p["wUV"]) ** 2) def hamaker(m1, m2, m3="vac", N=8000, smax=8): """A_132 by the Matsubara sum with multiple reflections to order smax.""" xn = np.arange(0, N) * xi1 e1, e2, e3 = eps(xn, MAT[m1]), eps(xn, MAT[m2]), eps(xn, MAT[m3]) x = ((e1 - e3) / (e1 + e3)) * ((e2 - e3) / (e2 + e3)) ssum = sum(x ** s / s ** 3 for s in range(1, smax + 1)) w = np.ones(N); w[0] = 0.5 # prime on the n=0 term return 1.5 * kT * np.sum(w * ssum) W_of_A = lambda A: A / (12 * np.pi * D0 ** 2) # work of adhesion gam_of_A = lambda A: A / (24 * np.pi * D0 ** 2) # surface free energy # ------------------------------------------ H-bond energy from water's own gamma A_ww = hamaker("water", "water") GAM_W_LW = gam_of_A(A_ww) GAM_W_AB = GAM_W - GAM_W_LW RHO_N = 997 / 0.018 * NA # water molecules / m^3 N_CROSS = RHO_N * Z_HB / 2 * L_HB * 0.5 # H-bonds crossing unit plane EPS_HB = 2 * GAM_W_AB / N_CROSS # effective free energy per bond # -------------------------------- group polarizabilities by Lorentz-Lorenz (L n_D) ALKANES = { # name: (n_D20, rho kg/m3, M g/mol, nCH3, nCH2) "pentane": (1.3575, 626.2, 72.15, 2, 3), "hexane": (1.3749, 659.5, 86.18, 2, 4), "heptane": (1.3878, 679.5, 100.20, 2, 5), "octane": (1.3974, 702.5, 114.23, 2, 6), "decane": (1.4102, 730.0, 142.28, 2, 8), "dodecane": (1.4216, 748.7, 170.33, 2, 10), } _rows, _y = [], [] for _n, _rho, _M, _a, _b in ALKANES.values(): _N = _rho / (_M * 1e-3) * NA _rows.append([_a, _b]) _y.append(3 * (_n ** 2 - 1) / (4 * np.pi * _N * (_n ** 2 + 2)) * 1e30) A_CH3, A_CH2 = np.linalg.lstsq(np.array(_rows, float), np.array(_y), rcond=None)[0] N_H, RHO_H, M_H = 1.3770, 764.0, 162.38 # hexamethyldisiloxane (L) _NH = RHO_H / (M_H * 1e-3) * NA A_HMDSO = 3 * (N_H ** 2 - 1) / (4 * np.pi * _NH * (N_H ** 2 + 2)) * 1e30 A_SI = (A_HMDSO - 6 * A_CH3 - 0.80) / 2 # alpha(O,siloxane) ~ 0.80 A^3 (L) RHO_A_HMDSO = A_HMDSO * _NH * 1e-30 # A^3 per A^3, bulk reference W_UV_A, W_IR_A = 1.87e16, 2.0e14 DE_IR_H = 2.05 - N_H ** 2 NF = 600 XN = np.arange(NF) * xi1 E_W, E_V, E_S = eps(XN, MAT["water"]), np.ones(NF), eps(XN, MAT["silica"]) def eps_layer(rho, xnn): """Clausius-Mossotti dielectric response of grafted material at density rho.""" cm = np.clip((4 * np.pi / 3) * np.atleast_1d(rho)[:, None] / (1 + (xnn[None, :] / W_UV_A) ** 2), 0, 0.95) return ((1 + 2 * cm) / (1 - cm) + DE_IR_H * (np.atleast_1d(rho)[:, None] / RHO_A_HMDSO) / (1 + (xnn[None, :] / W_IR_A) ** 2)) def lifshitz_stack(eL, dL, e_sub, e_med, e_out, D, xnn, xmax=50.0, nx=1200): """Work of adhesion for [substrate | layers | gap D | outer], recursive reflection coefficients through the stack. eL ordered TOP-most layer first.""" x = np.linspace(1e-6, xmax, nx) k = x[None, :] / (2 * D) Nl = eL.shape[0] Tm = np.repeat((eL[-1][:, None] - e_sub[:, None]) / (eL[-1][:, None] + e_sub[:, None]), nx, axis=1) for j in range(Nl - 2, -1, -1): ph = np.exp(-2 * k * dL[j + 1]) rj = ((eL[j][:, None] - eL[j + 1][:, None]) / (eL[j][:, None] + eL[j + 1][:, None])) Tm = (rj + Tm * ph) / (1 + rj * Tm * ph) ph = np.exp(-2 * k * dL[0]) r0 = (e_med[:, None] - eL[0][:, None]) / (e_med[:, None] + eL[0][:, None]) DL = (r0 + Tm * ph) / (1 + r0 * Tm * ph) DR = ((e_out - e_med) / (e_out + e_med))[:, None] integ = x[None, :] * np.log(np.clip(1 - DL * DR * np.exp(-x[None, :]), 1e-300, None)) Fn = (kT / (8 * np.pi * D ** 2)) * TRAPZ(integ, x, axis=1) Fn[0] *= 0.5 return np.sum(Fn) W_LW_BARE = W_of_A(hamaker("water", "silica")) T_LAYER = 4.40 # methyl-layer thickness, A def W_dense_film(sig_tms, t=T_LAYER): """Water on a grafted film of areal density sig_tms (/A^2) over silica.""" S = sig_tms * (3 * A_CH3 + A_SI) # areal polarizability, A^3/A^2 return lifshitz_stack(eps_layer(S / t, XN), np.array([t * 1e-10]), E_S, E_V, E_W, D0, XN) theta_from = lambda W: float(np.degrees(np.arccos(np.clip(W / GAM_W - 1, -1, 1)))) # ---------------------------------------------------------- grafting geometry L = 80.0 # periodic box edge, A pbcv = lambda d, Lb=L: d - Lb * np.round(d / Lb) tile3 = lambda P, Lb=L: np.vstack([P + np.array([i * Lb, j * Lb, 0]) for i in (-1, 0, 1) for j in (-1, 0, 1)]) def tripod(origin, axis, phi): """Three methyl carbons of -Si(CH3)3 grafted at origin along axis.""" ax = np.asarray(axis, float); ax = ax / np.linalg.norm(ax) tmp = np.array([0, 0, 1.0]) if abs(ax[2]) < 0.9 else np.array([1.0, 0, 0]) e1 = np.cross(ax, tmp); e1 /= np.linalg.norm(e1) e2 = np.cross(ax, e1) Si = np.asarray(origin, float) + h_Si * ax C = np.array([Si + d_SiC * (np.cos(np.pi - tet) * ax + np.sin(np.pi - tet) * (np.cos(phi + 2 * np.pi * k / 3) * e1 + np.sin(phi + 2 * np.pi * k / 3) * e2)) for k in range(3)]) return Si, C def densest_lattice_sigma(): """Upper bound: densest triangular lattice of TMS tripods on a flat plane.""" def min_CC(d, phi): _, C0 = tripod([0, 0, 0], [0, 0, 1], phi) R = [d * np.array([i + 0.5 * j, j * np.sqrt(3) / 2, 0]) for i in (-1, 0, 1) for j in (-1, 0, 1) if (i, j) != (0, 0)] return min(np.linalg.norm(C0[:, None, :] - (r + C0)[None, :, :], axis=-1).min() for r in R) best = lambda d: max(min_CC(d, p) for p in np.linspace(0, 2 * np.pi / 3, 60)) d_min = brentq(lambda d: best(d) - 2 * r_CH3, 4.0, 12.0, xtol=1e-3) return 2 / (np.sqrt(3) * d_min ** 2), d_min def silanol_field(sig_nm2, seed, Lb=L, dmin=3.0): """Amorphous silanol positions: hard-disc RSA at areal density sig_nm2.""" rg = np.random.default_rng(seed) tgt = int(sig_nm2 / 100 * Lb * Lb) pts = np.zeros((0, 2)) for _ in range(400000): if len(pts) >= tgt: break p = rg.uniform(0, Lb, 2) if len(pts): d = pbcv(pts - p, Lb) if np.sqrt((d ** 2).sum(1)).min() < dmin: continue pts = np.vstack([pts, p]) return pts def graft_flat(sitf, seed, Lb=L, n_try=80): """RSA grafting of tripods onto a flat silanol field. Si-O-Si hinge allows up to 1 A of lateral offset; thermal tilt sigma = 15 deg.""" rg = np.random.default_rng(seed) gC, gSi, gr = [], [], np.zeros(len(sitf), bool) for idx in rg.permutation(len(sitf)): tree = cKDTree(tile3(np.vstack(gC), Lb)) if gC else None for _ in range(n_try): off, oa = rg.uniform(0, 1.0), rg.uniform(0, 2 * np.pi) xy = sitf[idx] + off * np.array([np.cos(oa), np.sin(oa)]) ax = np.array([0, 0, 1.0]) + rg.normal(0, np.deg2rad(15), 3) * [1, 1, 0] _, C = tripod([xy[0], xy[1], 0.0], ax, rg.uniform(0, 2 * np.pi)) Cw = C.copy(); Cw[:, :2] %= Lb if tree is not None and min(len(tree.query_ball_point(c, 2 * r_CH3)) for c in Cw) > 0: continue gC.append(Cw) gSi.append([xy[0] % Lb, xy[1] % Lb, h_Si]) gr[idx] = True break return (gr, np.vstack(gC) if gC else np.zeros((0, 3)), np.array(gSi) if gSi else np.zeros((0, 3))) # H-bond cone directions: upper hemisphere, fixed seed for reproducibility _u = np.random.default_rng(3).normal(size=(3000, 3)) _u /= np.linalg.norm(_u, axis=1)[:, None] U3 = _u[_u[:, 2] > 0] A_BARE = float(((d_HB * U3)[:, 2] > 1.0).mean()) # bare-surface reference def accessibility(sites, idx, C_, Si_, Lb=L, h=None): """nu: fraction of H-bond cone directions reaching bulk water, relative to the bare surface. A direction is blocked if the water O would overlap a methyl or a grafted Si.""" tC = cKDTree(tile3(C_, Lb)) if len(C_) else None tS = cKDTree(tile3(Si_, Lb)) if len(Si_) else None out = [] for i in idx: base = sites[i] if sites.shape[1] == 3 else np.array([sites[i][0], sites[i][1], 0.0]) W = base + d_HB * U3 z0 = base[2] ok = W[:, 2] > z0 + 1.0 Wq = W.copy(); Wq[:, :2] %= Lb if tC is not None: ok &= (tC.query(Wq)[0] > r_CH3 + r_Ow) if tS is not None: ok &= (tS.query(Wq)[0] > r_Si + r_Ow) out.append(ok.mean() / A_BARE) return np.array(out) def methyl_coverage(C_, seed, Lb=L, n=25000): """Projected areal fraction covered by methyl vdW discs.""" if not len(C_): return 0.0 g = np.random.default_rng(seed).uniform(0, Lb, (n, 2)) d = pbcv(g[:, None, :] - C_[None, :, :2], Lb) return float((np.sqrt((d ** 2).sum(-1)) < r_CH3).any(1).mean()) # --------------------------------------------------------- corrugated surfaces def rough_surface(rms, corr_len, seed, Lb=L, n=128): """Gaussian random rough surface, Gaussian height autocorrelation.""" rg = np.random.default_rng(seed) kx = 2 * np.pi * np.fft.fftfreq(n, d=Lb / n) KX, KY = np.meshgrid(kx, kx, indexing="ij") K = np.sqrt(KX ** 2 + KY ** 2) h = np.real(np.fft.ifft2(np.fft.fft2(rg.normal(size=(n, n))) * np.exp(-(K * corr_len / 2.0) ** 2))) return h * (rms / h.std()) def h_and_normal(h, x, y, Lb=L): """Bilinear height and outward unit normal at (x, y).""" n = h.shape[0]; dx = Lb / n fx, fy = (x % Lb) / dx, (y % Lb) / dx i0, j0 = np.floor(fx).astype(int) % n, np.floor(fy).astype(int) % n i1, j1 = (i0 + 1) % n, (j0 + 1) % n tx, ty = fx - np.floor(fx), fy - np.floor(fy) h00, h10, h01, h11 = h[i0, j0], h[i1, j0], h[i0, j1], h[i1, j1] z = h00 * (1 - tx) * (1 - ty) + h10 * tx * (1 - ty) + h01 * (1 - tx) * ty + h11 * tx * ty dhdx = ((h10 - h00) * (1 - ty) + (h11 - h01) * ty) / dx dhdy = ((h01 - h00) * (1 - tx) + (h11 - h10) * tx) / dx nrm = np.stack([-dhdx, -dhdy, np.ones_like(z)], axis=-1) return z, nrm / np.linalg.norm(nrm, axis=-1, keepdims=True) def graft_rough(rms, corr_len, seed, sig_site=4.6, Lb=L, n_try=90): """Silanols on a corrugated surface; tripods graft along the LOCAL normal.""" h = rough_surface(rms, corr_len, seed, Lb) if rms > 0 else np.zeros((128, 128)) rg = np.random.default_rng(seed + 77) tgt = int(sig_site / 100 * Lb * Lb) P = np.zeros((0, 3)) for _ in range(400000): if len(P) >= tgt: break p2 = rg.uniform(0, Lb, 2) z, _ = h_and_normal(h, p2[0], p2[1], Lb) p = np.array([p2[0], p2[1], float(z)]) if len(P): d = P - p; d[:, :2] = pbcv(d[:, :2], Lb) if np.sqrt((d ** 2).sum(1)).min() < 3.0: continue P = np.vstack([P, p]) gC, gSi, gr = [], [], np.zeros(len(P), bool) for idx in rg.permutation(len(P)): tr = cKDTree(tile3(np.vstack(gC), Lb)) if gC else None for _ in range(n_try): off, oa = rg.uniform(0, 1.0), rg.uniform(0, 2 * np.pi) xy = P[idx][:2] + off * np.array([np.cos(oa), np.sin(oa)]) zz, nn = h_and_normal(h, xy[0], xy[1], Lb) ax = nn + rg.normal(0, 0.22, 3); ax /= np.linalg.norm(ax) Si, C = tripod([xy[0], xy[1], float(zz)], ax, rg.uniform(0, 2 * np.pi)) Cw = C.copy(); Cw[:, :2] %= Lb if tr is not None and min(len(tr.query_ball_point(c, 2 * r_CH3)) for c in Cw) > 0: continue gC.append(Cw) Siw = Si.copy(); Siw[:2] %= Lb gSi.append(Siw) gr[idx] = True break gx, gy = np.gradient(h, Lb / 128, Lb / 128) area_ratio = float(np.mean(np.sqrt(1 + gx ** 2 + gy ** 2))) return (h, P, gr, np.vstack(gC) if gC else np.zeros((0, 3)), np.array(gSi) if gSi else np.zeros((0, 3)), area_ratio) def calcine(sitf, d_cond, seed, Lb=L): """Random sequential pairwise condensation 2 SiOH -> Si-O-Si + H2O.""" rg = np.random.default_rng(seed) alive = np.ones(len(sitf), bool) nbr = 0 for i in rg.permutation(len(sitf)): if not alive[i]: continue d = pbcv(sitf - sitf[i], Lb) r = np.sqrt((d ** 2).sum(1)) cand = np.where(alive & (r < d_cond) & (np.arange(len(sitf)) != i))[0] if not len(cand): continue j = cand[rg.integers(len(cand))] alive[i] = alive[j] = False nbr += 1 return alive, nbr def wetting(s_tms, s_res, nu, fcov, W_dense, s_silox=0.0, beta=0.0): """Assemble the work of adhesion and invert Young-Dupre.""" W_lw = fcov * W_dense + (1 - fcov) * W_LW_BARE W_ab = ((s_res * nu * N_HB_PER_OH + s_silox * nu * N_HB_SILOX * beta) * 1e20 * EPS_HB) return theta_from(W_lw + W_ab), W_lw, W_ab # ============================================================== main def main(): R = {} # -- validation of the Lifshitz machinery against known Hamaker constants BENCH = {("water", "water", "vac"): 3.70, ("alkane", "alkane", "vac"): 5.00, ("silica", "silica", "vac"): 6.50, ("tms", "tms", "water"): 0.45, ("silica", "silica", "water"): 0.85} R["hamaker"] = [{"pair": "|".join(k), "model_1e20J": hamaker(*k) * 1e20, "reference_1e20J": v, "rel_err": hamaker(*k) * 1e20 / v - 1} for k, v in BENCH.items()] # -- Matsubara truncation and quadrature convergence R["conv_matsubara"] = [{"N": N, "A_1e20J": hamaker("water", "silica", N=N) * 1e20} for N in [50, 100, 200, 500, 1000, 2000, 4000, 8000]] R["conv_reflection"] = [{"smax": s, "A_1e20J": hamaker("silica", "silica", "water", smax=s) * 1e20} for s in [1, 2, 4, 8, 16]] e_one = eps(XN, MAT["silica"])[None, :] R["conv_quad"] = [{"nx": nx, "xmax": xm, "W_mJm2": lifshitz_stack(e_one, np.array([1e-5]), E_S, E_V, E_W, D0, XN, xmax=xm, nx=nx) * 1e3} for nx, xm in [(200, 30.), (400, 40.), (800, 50.), (1200, 50.), (2400, 60.), (4800, 80.)]] R["stack_vs_analytic"] = { "silica_stack_mJm2": lifshitz_stack(e_one, np.array([1e-5]), E_S, E_V, E_W, D0, XN) * 1e3, "silica_analytic_mJm2": W_LW_BARE * 1e3, "tms_stack_mJm2": lifshitz_stack(eps(XN, MAT["tms"])[None, :], np.array([1e-5]), eps(XN, MAT["tms"]), E_V, E_W, D0, XN) * 1e3, "tms_analytic_mJm2": W_of_A(hamaker("water", "tms")) * 1e3} R["calibration"] = { "A_ww_1e20J": A_ww * 1e20, "gam_w_LW_mJm2": GAM_W_LW * 1e3, "gam_w_AB_mJm2": GAM_W_AB * 1e3, "n_cross_m2": N_CROSS, "eps_HB_J": EPS_HB, "eps_HB_kJmol": EPS_HB * NA / 1e3, "alpha_CH3_A3": A_CH3, "alpha_CH2_A3": A_CH2, "alpha_Si_A3": A_SI, "alpha_HMDSO_A3": A_HMDSO, "rho_a_HMDSO": RHO_A_HMDSO, "W_LW_bare_mJm2": W_LW_BARE * 1e3} # -- packing bounds on a flat plane sig_max, d_min = densest_lattice_sigma() rsa = [] for sd in (1, 2, 3): F = silanol_field(4.6, sd) gr, C_, Si_ = graft_flat(F, 100 + sd) rsa.append(gr.sum() / (L * L)) R["packing"] = {"sigma_max_nm2": sig_max * 100, "d_min_A": d_min, "sigma_rsa_nm2": float(np.mean(rsa)) * 100, "sigma_rsa_sd_nm2": float(np.std(rsa)) * 100, "n_replicates": len(rsa)} # -- reference configuration: fully hydroxylated, flat, jammed F = silanol_field(4.6, 1) gr, C_, Si_ = graft_flat(F, 101) nu_all = accessibility(F, np.where(~gr)[0], C_, Si_) fcov = methyl_coverage(C_, 5) W_dense = W_dense_film(gr.sum() / (L * L)) th, W_lw, W_ab = wetting((gr.sum()) / (L * L), (~gr).sum() / (L * L), float(nu_all.mean()), fcov, W_dense) R["reference"] = { "sigma_TMS_nm2": gr.sum() / (L * L) * 100, "sigma_OH_res_nm2": (~gr).sum() / (L * L) * 100, "f_OH_floor": float((~gr).sum() / (L * L) / SIG_OH0), "nu_mean": float(nu_all.mean()), "buried_frac": float((nu_all < 0.02).mean()), "sigma_OH_eff_nm2": (~gr).sum() / (L * L) * float(nu_all.mean()) * 100, "f_CH3": fcov, "W_LW_mJm2": W_lw * 1e3, "W_AB_mJm2": W_ab * 1e3, "W_adh_mJm2": (W_lw + W_ab) * 1e3, "theta_deg": th, "W_LW_dense_mJm2": W_dense * 1e3, "theta_if_accessible": theta_from( W_lw + (~gr).sum() / (L * L) * 1.0 * N_HB_PER_OH * EPS_HB * 1e20)} # -- sensitivity of the reference angle sens = {"baseline": th} for lbl, D in [("D0=1.60A", 1.60e-10), ("D0=1.70A", 1.70e-10)]: A_wt = hamaker("water", "tms") sens[lbl] = theta_from(A_wt / (12 * np.pi * D ** 2) * 0.90 + W_ab) for t in (4.0, 5.0): sens[f"t={t}A"] = theta_from(W_dense_film(gr.sum() / (L * L), t) + W_ab) sens["nu x3"] = theta_from(W_lw + 3 * W_ab) sens["sigma_TMS -15%"] = theta_from(0.9 * W_dense + 0.1 * W_LW_BARE + 3 * W_ab) R["sensitivity"] = sens # -- theta versus residual silanol fraction (flat, 3 configurations) cfgs = [] for sd in (1, 2, 3): Fs = silanol_field(4.6, sd) cfgs.append((Fs,) + graft_flat(Fs, 100 + sd)) sweep = [] for frac in np.linspace(0, 1, 13): per = [] for k, (Fs, g0, Ca, Sa) in enumerate(cfgs): gidx = np.where(g0)[0] rg2 = np.random.default_rng(100 + 7 * k) keep = rg2.permutation(len(gidx))[:int(round(frac * len(gidx)))] m = np.zeros(len(gidx), bool); m[keep] = True Cs = Ca.reshape(-1, 3, 3)[m].reshape(-1, 3) if m.sum() else np.zeros((0, 3)) Ss = Sa[m] free = np.setdiff1d(np.arange(len(Fs)), gidx[m]) nu = accessibility(Fs, free, Cs, Ss) per.append((m.sum() / (L * L), len(free) / (L * L), float(nu.mean()) if len(nu) else 0.0, methyl_coverage(Cs, 200 + k))) s_t, s_o, nu_m, fc = np.mean(per, axis=0) nu_sd = float(np.std([p[2] for p in per])) thh, wl, wa = wetting(s_t, s_o, nu_m, fc, W_dense_film(s_t) if s_t > 0 else W_LW_BARE) sweep.append({"f_OH": s_o / SIG_OH0, "sigma_OH_nm2": s_o * 100, "sigma_TMS_nm2": s_t * 100, "f_CH3": fc, "nu": nu_m, "nu_sd": nu_sd, "W_LW_mJm2": wl * 1e3, "W_AB_mJm2": wa * 1e3, "theta_deg": thh}) R["silanol_sweep"] = sweep # -- closed form: random-overlap coverage and cone-screening exponent sT = np.array([s["sigma_TMS_nm2"] for s in sweep]) / 100 fcv = np.array([s["f_CH3"] for s in sweep]) nuv = np.array([s["nu"] for s in sweep]) A_eff = curve_fit(lambda s, a: 1 - np.exp(-a * s), sT, fcv, p0=[40])[0][0] ok = fcv < 0.98 m_exp = curve_fit(lambda f, m: (1 - f) ** m, fcv[ok], nuv[ok], p0=[2.0])[0][0] def theta_closed(f_OH): s_t = SIG_OH0 * (1 - f_OH) fc = 1 - np.exp(-A_eff * s_t) nu = (1 - fc) ** m_exp return wetting(s_t, SIG_OH0 * f_OH, nu, fc, W_dense_film(max(s_t, 1e-9)))[0] f_floor = R["reference"]["f_OH_floor"] th_c = np.array([theta_closed(s["f_OH"]) for s in sweep]) th_m = np.array([s["theta_deg"] for s in sweep]) R["closed_form"] = { "A_eff_A2": float(A_eff), "methyl_shadow_A2": 3 * np.pi * r_CH3 ** 2, "m_exponent": float(m_exp), "rms_deg": float(np.sqrt(np.mean((th_c - th_m) ** 2))), "max_dev_deg": float(np.abs(th_c - th_m).max()), "f_OH_at_90deg": float(brentq(lambda f: theta_closed(f) - 90, 0.30, 0.999)), "f_OH_at_110deg": float(brentq(lambda f: theta_closed(f) - 110, 0.30, 0.999))} # -- floor: site-limited or sterically limited? R["floor_scan"] = [] for s0 in [4.6, 4.0, 3.4, 2.8, 2.2, 1.6, 1.0, 0.6]: per = [] for sd in (11, 12): Fq = silanol_field(s0, sd) gq, _, _ = graft_flat(Fq, 100 + sd) per.append((len(Fq) / (L * L), gq.sum() / (L * L), gq.mean())) m = np.mean(per, axis=0) R["floor_scan"].append({"sigma_sites_nm2": m[0] * 100, "sigma_TMS_nm2": m[1] * 100, "frac_capped": m[2]}) # -- calcination SISI = 2 * d_SiO * np.sin(np.deg2rad(144 / 2)) R["siloxane_bridge_A"] = SISI R["calcine_only"], R["calcine_then_silylate"] = [], [] for dc in [3.2, 3.5, 4.0, 4.5, 5.0, 6.0]: per = [] for sd in (21, 22, 23): Fq = silanol_field(4.6, sd) al, nbr = calcine(Fq, dc, 200 + sd) per.append((al.sum() / (L * L), nbr / (L * L))) s_oh, s_br = np.mean(per, axis=0) row = {"d_cond_A": dc, "sigma_OH_nm2": s_oh * 100, "sigma_silox_nm2": s_br * 100} for b in (0.0, 0.2, 0.4): row[f"theta_beta{b}"] = wetting(0.0, s_oh, 1.0, 0.0, W_LW_BARE, s_silox=s_br, beta=b)[0] R["calcine_only"].append(row) per2 = [] for sd in (21, 22): Fq = silanol_field(4.6, sd) al, nbr = calcine(Fq, dc, 200 + sd) Fs = Fq[al] if len(Fs) < 3: continue gq, Cq, Sq = graft_flat(Fs, 300 + sd) nu = accessibility(Fs, np.where(~gq)[0], Cq, Sq) per2.append((gq.sum() / (L * L), (~gq).sum() / (L * L), float(nu.mean()) if len(nu) else 0.0, methyl_coverage(Cq, 400 + sd), nbr / (L * L))) s_t, s_r, nu_m, fc, s_b = np.mean(per2, axis=0) thh, wl, wa = wetting(s_t, s_r, nu_m, fc, W_dense_film(max(s_t, 1e-9)), s_silox=s_b, beta=0.2) R["calcine_then_silylate"].append( {"d_cond_A": dc, "sigma_OH_in_nm2": s_oh * 100, "sigma_TMS_nm2": s_t * 100, "sigma_OH_res_nm2": s_r * 100, "nu": nu_m, "f_CH3": fc, "W_LW_mJm2": wl * 1e3, "theta_deg": thh}) # -- amorphous corrugation R["rough"] = [] for rms, corr in [(0.0, 0.0), (0.5, 8.0), (1.0, 8.0), (1.5, 8.0), (2.0, 8.0), (2.5, 8.0), (1.5, 4.0), (1.5, 15.0)]: per = [] for sd in (41, 42): h, P, gq, Cq, Sq, ar = graft_rough(rms, corr, sd) nu = accessibility(P, np.where(~gq)[0], Cq, Sq) per.append((len(P) / (L * L), gq.sum() / (L * L), gq.mean(), ar, (~gq).sum() / (L * L), float(nu.mean()) if len(nu) else 0.0, methyl_coverage(Cq, 600 + sd))) m = np.mean(per, axis=0) thh, wl, wa = wetting(m[1], m[4], m[5], m[6], W_dense_film(m[1])) R["rough"].append({"rms_A": rms, "corr_len_A": corr, "sigma_sites_nm2": m[0] * 100, "sigma_TMS_nm2": m[1] * 100, "frac_capped": m[2], "area_ratio": m[3], "sigma_TMS_per_true_area_nm2": m[1] * 100 / m[3], "nu": m[5], "f_CH3": m[6], "W_LW_mJm2": wl * 1e3, "W_AB_mJm2": wa * 1e3, "theta_intrinsic_deg": thh, "theta_wenzel_deg": theta_from( GAM_W * (m[3] * (np.cos(np.radians(thh)) + 1) - 1) + GAM_W) if False else float(np.degrees(np.arccos(np.clip( m[3] * np.cos(np.radians(thh)), -1, 1))))}) per_true = np.array([r["sigma_TMS_per_true_area_nm2"] for r in R["rough"]]) R["rough_invariance"] = {"mean_nm2": float(per_true.mean()), "sd_nm2": float(per_true.std()), "cv_percent": float(100 * per_true.std() / per_true.mean())} # -- silane-excluded micropores r_tri = d_SiC * np.sin(tet / 2) + r_CH3 ref = R["reference"] def theta_pore(f_pore): frac = 1.0 - f_pore s_t = ref["sigma_TMS_nm2"] / 100 * frac s_r = ref["sigma_OH_res_nm2"] / 100 * frac fc = ref["f_CH3"] * frac W_lw = fc * W_dense + (1 - fc) * W_LW_BARE W_ab = ((s_r * ref["nu_mean"] + SIG_OH0 * f_pore * 1.0) * 1e20 * N_HB_PER_OH * EPS_HB) return theta_from(W_lw + W_ab), W_lw, W_ab R["micropore"] = {"tripod_radius_A": r_tri, "silane_min_pore_A": 2 * r_tri, "water_min_pore_A": 2 * r_Ow, "f_pore_at_90deg": float(brentq(lambda f: theta_pore(f)[0] - 90, 0.001, 0.6)), "curve": [{"f_pore": f, "theta_deg": theta_pore(f)[0], "W_LW_mJm2": theta_pore(f)[1] * 1e3, "W_AB_mJm2": theta_pore(f)[2] * 1e3} for f in [0.0, 0.02, 0.05, 0.10, 0.15, 0.20, 0.30, 0.40]]} with open("results.json", "w") as fh: json.dump(R, fh, indent=1, default=float) print(f"theta(reference) = {R['reference']['theta_deg']:.1f} deg") print(f"f_OH floor = {R['reference']['f_OH_floor']:.3f}") print(f"nu = {R['reference']['nu_mean']:.4f}") print(f"per-true-area CV = {R['rough_invariance']['cv_percent']:.1f}%") return R if __name__ == "__main__": main()