"""Reproduce #1004: public terminal RHF to raw MP2 and two independent replays.""" import json from hashlib import sha256 from itertools import product from pathlib import Path import numpy as np from scipy.linalg import eigh try: from . import bipole_finite_operator as op except ImportError: # direct script invocation import bipole_finite_operator as op from vibeqc import run_pbc_bipole_rhf from vibeqc.dlpno.translation_driver import run_translation_dlpno_mp2 MESH = (3, 1, 1) LIMIT = 64 << 20 def bloch(r): """Analytic Gaussian transition integrals. Reverse the *second AO pair* relative to the native reader, hence (q,l,r)=(i-a,a,j), ordered (i,a,b,j). No native integral, pair matrix or amplitude enters this energy sum. """ c, e = map(np.array, (r.mo_coeffs, r.mo_energies)) n = len(c) tiles = {} err = forbidden = 0.0 kmax = r.resolved_domain.reciprocal_cutoff_bohr_inv for i, a, j in product(range(n), repeat=3): b = (i + j - a) % n q = op.centered((i - a) % n, n) _, ps, w = op.reciprocal(q, kmax) z = 0j for u, v, x, y in product(range(2), repeat=4): parts, _ = op.parts((u, v, y, x), q, op.centered(a, n), op.centered(j, n), ps, w) z += c[i, u, 0].conj() * c[a, v, 1] * c[j, x, 0].conj() * c[b, y, 1] * sum(parts) / n tiles[i, a, j, b] = z h = r.terminal_transitions err = max(err, abs(z - h.transition(i, a, j, b, max_bytes=LIMIT))) for wrong in set(range(n)) - {b}: forbidden = max(forbidden, abs(h.transition(i, a, j, wrong, max_bytes=LIMIT))) op.near(err, 0.0, 1e-10) op.near(forbidden, 0.0, 1e-12) by_q = np.zeros(n) uncrossed = 0.0 for (i, a, j, b), z in tiles.items(): den = e[i, 0] + e[j, 0] - e[a, 1] - e[b, 1] assert den < 0 # Both ERIs already carry 1/N. Only this final extensive sum is # divided by N. q here labels the excitation a-i, not native i-a. by_q[(a - i) % n] += (z.conjugate() * (2 * z - tiles[i, b, j, a])).real / (n * den) uncrossed += abs(z) ** 2 / (n * den) return ( by_q, err, forbidden, uncrossed, [[i, a, j, b, z.real, z.imag] for (i, a, j, b), z in tiles.items()], ) def explicit_torus(r): """Six AO BvK cell: accumulate erfc image integrals modulo three, use a *single supercell* reciprocal lattice and metric, then diagonalize F C=S C e. Neither this construction nor its canonical MP2 loop reads a momentum tile. """ n = 3 lengths = op.LENGTH * [n, 1, 1] volume = np.prod(lengths) kmax = r.resolved_domain.reciprocal_cutoff_bohr_inv bound = np.ceil(kmax * lengths / (2 * np.pi)).astype(int) labels = np.asarray(list(product(*(range(-m, m + 1) for m in bound)))) ps = 2 * np.pi * labels / lengths p2 = np.sum(ps**2, axis=1) keep = (p2 > 0) & (p2 <= kmax**2) ps, p2 = ps[keep], p2[keep] assert len(ps) == 226 + 2 * 218 w = 4 * np.pi * np.exp(-p2 / (4 * op.OMEGA**2)) / (volume * p2) ft = np.zeros((6, 6, len(ps)), complex) s = np.zeros((6, 6)) sr = np.zeros((6,) * 4) for t, u in product(range(n), range(2)): for v in range(2): for g in op.pair_cells(u, v): overlap, center = op.pair(u, v, [g, 0, 0]) row, col = 2 * t + u, 2 * ((t + g) % n) + v s[row, col] += overlap ft[row, col] += overlap * np.exp( -1j * (ps @ (center + [t * op.LENGTH[0], 0, 0])) - p2 / (8 * op.ALPHA) ) for x, y in product(range(2), repeat=2): for (g, p, h), val in op.images(u, v, x, y): sr[ 2 * t + u, 2 * ((t + g) % n) + v, 2 * ((t + p) % n) + x, 2 * ((t + h) % n) + y, ] += val lr = np.einsum("uvg,g,xyg->uvxy", ft.conj(), w, ft) op.near(lr.imag, 0.0, 1e-12) eri = sr + lr.real - np.pi / (volume * op.OMEGA**2) * np.einsum("uv,xy->uvxy", s, s) def fold(mat): # Use the actual Cartesian k vectors and AO-cell displacements. cells = np.arange(3) * op.LENGTH[0] phase = np.exp(1j * (cells[:, None] - cells[None, :])[:, :, None] * r.kpoints_cart[:, 0]) a = np.einsum("rtk,k,kab->ratb", phase, r.kpoint_weights, mat).reshape(6, 6) op.near(a.imag, 0.0, 1e-12) return a.real op.near(s, fold(r.overlap), 1e-12) f = fold(r.fock) eps, c = eigh(f, s) op.near(eps, np.sort(np.ravel(r.mo_energies)), 1e-10) op.near(2 * c[:, :3] @ c[:, :3].T, fold(r.accepted_density_k), 1e-10) den = ( eps[:3, None, None, None] - eps[None, 3:, None, None] + eps[None, None, :3, None] - eps[None, None, None, 3:] ) def canonical(o, v): tiles = np.einsum("ui,va,xj,yb,uvxy->iajb", o.conj(), v, o.conj(), v, eri, optimize=True) return float(np.sum((tiles.conj() * (2 * tiles - tiles.swapaxes(1, 3))).real / den) / n) raw = canonical(c[:, :3], c[:, 3:]) # The +/- third levels coincide to roundoff. Rotate each whole equal-energy # subspace with a complex unitary; neither projector nor energy may change. rot = np.array([[1, 1j], [1j, 1]]) / np.sqrt(2.0) o, v = c[:, :3].astype(complex), c[:, 3:].astype(complex) for z, ee in ((o, eps[:3]), (v, eps[3:])): left = np.argmin(np.diff(ee)) assert np.diff(ee)[left] < 1e-12 z[:, left : left + 2] = z[:, left : left + 2] @ rot op.near(o @ o.conj().T, c[:, :3] @ c[:, :3].T, 1e-12) op.near(v @ v.conj().T, c[:, 3:] @ c[:, 3:].T, 1e-12) rotation_error = canonical(o, v) - raw op.near(rotation_error, 0.0, 1e-13) return raw, rotation_error def run(): system, basis, mesh, options, kw = op.rhf_fixture(3) r = run_pbc_bipole_rhf(system, basis, mesh, options, terminal_transition_bytes=LIMIT, **kw) result = run_translation_dlpno_mp2( r, MESH, factor_handle=r.terminal_transitions, max_factor_bytes=LIMIT, n_frozen=0, tcut_mkn=0, tcut_pairs=0, tcut_pno=0, conv_tol_energy=1e-13, conv_tol_residual=1e-12, ) sol = result.solution assert ( sol.converged and result.pair_build.same_source and set(result.pair_build.pno_ranks) == {3} ) op.near(sol.max_residual, 0.0, 1e-12) op.near(sol.trace[-1].energy_change, 0.0, 1e-13) q, selected, forbidden, uncrossed, tiles = bloch(r) torus, rotation = explicit_torus(r) raw = sol.e_pair_per_cell op.near(raw, q.sum(), 1e-10) op.near(raw, torus, 1e-10) omitted = q[:2].sum() assert abs(omitted - raw) > 1e-10 assert abs(uncrossed - raw) > 1e-10 orbitals = r.terminal_transitions.orbitals def complex_array(values): return np.stack((values.real, values.imag), axis=-1).tolist() return { "source": r.terminal_transitions.identity, "solver": result.identity, "source_record": r.terminal_transitions.provenance, "terminal_arrays_re_im": { key: complex_array(getattr(orbitals, key)) for key in ( "fock", "overlap", "coefficients", "density", "occupied_projector", "virtual_projector", ) }, "bloch_i_a_j_b_re_im": tiles, "using_direct_for_exchange_minus_raw": uncrossed - raw, "raw": raw, "ordered_bloch": q.sum(), "explicit_supercell": torus, "raw_minus_bloch": raw - q.sum(), "raw_minus_supercell": raw - torus, "q_zero_plus_minus": q.tolist(), "omit_minus": omitted, "omit_minus_minus_raw": omitted - raw, "max_selected_error": selected, "forbidden": forbidden, "rotated_minus_supercell": rotation, "iterations": sol.n_iter, "max_residual": sol.max_residual, "final_energy_change": sol.trace[-1].energy_change, "hf": r.energy, "orbital_energies": np.asarray(r.mo_energies).tolist(), "source_reservation_bytes": r.terminal_transitions.required_bytes, "full_pair_ranks": result.pair_build.pno_ranks, "diagnostic_sha256": sha256(Path(__file__).read_bytes()).hexdigest(), } if __name__ == "__main__": print(json.dumps(run()), flush=True)