feat(M11.3): add QDO-calibrated LJ and ATM interactions

openwave/xperiments/m11_cat_ept_soliton_qdo/__init__.py

@@ -1,15 +1,11 @@
 """M11 CAT/EPT pointwise soliton--Liouville--QDO particle model."""
-from .pointwise_soliton_carrier_m111 import (
-    PointwiseSolitonConfig, PointwiseSolitonState,
-    construct_pointwise_soliton, run_pointwise_soliton_study,
-)
-from .liouville_soliton_tensor_m112 import (
-    LiouvilleTensorConfig, LiouvilleTensorState,
-    construct_liouville_tensor, run_liouville_tensor_study,
-)
+from .pointwise_soliton_carrier_m111 import PointwiseSolitonConfig, PointwiseSolitonState, construct_pointwise_soliton, run_pointwise_soliton_study
+from .liouville_soliton_tensor_m112 import LiouvilleTensorConfig, LiouvilleTensorState, construct_liouville_tensor, run_liouville_tensor_study
+from .qdo_lj_atm_interaction_m113 import QDOLJATMConfig, run_qdo_lj_atm_study
 __all__ = [
     "PointwiseSolitonConfig", "PointwiseSolitonState",
-    "LiouvilleTensorConfig", "LiouvilleTensorState",
+    "LiouvilleTensorConfig", "LiouvilleTensorState", "QDOLJATMConfig",
     "construct_pointwise_soliton", "run_pointwise_soliton_study",
     "construct_liouville_tensor", "run_liouville_tensor_study",
+    "run_qdo_lj_atm_study",
 ]

openwave/xperiments/m11_cat_ept_soliton_qdo/formal/qdo_lj_atm_interaction.v1.json

@@ -0,0 +1,21 @@
+{
+  "schema": "openwave.m11.qdo-lj-atm-interaction-ledger.v1",
+  "model_id": "M11",
+  "milestone": "M11.3",
+  "formal_repository": "jagg-ix/entropic-physlib-private",
+  "formal_branch": "entropic-physlib-linear-full",
+  "formal_head": "8bafa9ab93cbb39e85909fc3837bb4b6e0dec748",
+  "equations": [
+    {"id": "lj-split", "equation": "E_LJ=ljRepulsive-ljAttractive", "formal_theorem": "ljPotential_split"},
+    {"id": "lj-tail", "equation": "-ljAttractive=-C6_LJ/r^6", "formal_theorem": "lj_tail_dispersion_form"},
+    {"id": "qdo-c6", "equation": "C6_LJ=C6_QDO=3/4 alpha1^2 hbar omega", "formal_theorem": "lj_tail_is_qdo_C6"},
+    {"id": "qdo-c9", "equation": "alpha1 C6=4 C9", "formal_theorem": "lattice_couplings_qdo_related"}
+  ],
+  "formal_sources": [
+    {"path": "Physlib/QuantumMechanics/ComplexAction/Particles/LennardJonesAxilrodTeller.lean", "sha": "93f5e1b55d1fbfe44dc5b8106b7331d08bdb5159", "theorem": "ljPotential_split"},
+    {"path": "Physlib/QuantumMechanics/ComplexAction/Particles/LennardJonesAxilrodTeller.lean", "sha": "93f5e1b55d1fbfe44dc5b8106b7331d08bdb5159", "theorem": "lj_tail_dispersion_form"},
+    {"path": "Physlib/QuantumMechanics/ComplexAction/Particles/LennardJonesAxilrodTeller.lean", "sha": "93f5e1b55d1fbfe44dc5b8106b7331d08bdb5159", "theorem": "lj_tail_eq_vdw_dispersionEnergy"},
+    {"path": "Physlib/QuantumMechanics/ComplexAction/Particles/LennardJonesAxilrodTeller.lean", "sha": "93f5e1b55d1fbfe44dc5b8106b7331d08bdb5159", "theorem": "lj_tail_is_qdo_C6"},
+    {"path": "Physlib/QuantumMechanics/ComplexAction/Particles/LennardJonesAxilrodTeller.lean", "sha": "93f5e1b55d1fbfe44dc5b8106b7331d08bdb5159", "theorem": "lattice_couplings_qdo_related"}
+  ]
+}

openwave/xperiments/m11_cat_ept_soliton_qdo/qdo_lj_atm_interaction_m113.py

@@ -0,0 +1,231 @@
+"""M11.3 QDO-calibrated Lennard-Jones and Axilrod-Teller interactions."""
+from __future__ import annotations
+
+from dataclasses import asdict, dataclass
+from hashlib import sha256
+import itertools
+import json
+import math
+from typing import Any, Mapping
+
+import numpy as np
+
+MILESTONE = "M11.3"
+SCHEMA = "openwave.m11.qdo-lj-atm-interaction.v1"
+FORMAL_HEAD = "8bafa9ab93cbb39e85909fc3837bb4b6e0dec748"
+FORMAL_SOURCE_PATH = (
+    "Physlib/QuantumMechanics/ComplexAction/Particles/"
+    "LennardJonesAxilrodTeller.lean"
+)
+FORMAL_SOURCE_SHA = "93f5e1b55d1fbfe44dc5b8106b7331d08bdb5159"
+FORMAL_SOURCES = tuple(
+    {
+        "path": FORMAL_SOURCE_PATH,
+        "sha": FORMAL_SOURCE_SHA,
+        "theorem": theorem,
+    }
+    for theorem in (
+        "ljPotential_split",
+        "lj_tail_dispersion_form",
+        "lj_tail_eq_vdw_dispersionEnergy",
+        "lj_tail_is_qdo_C6",
+        "lattice_couplings_qdo_related",
+        "atm_equilateral",
+        "atm_linear",
+    )
+)
+
+
+def _canonical_json(value: Mapping[str, Any]) -> str:
+    return json.dumps(value, sort_keys=True, separators=(",", ":"), default=str)
+
+
+@dataclass(frozen=True)
+class QDOLJATMConfig:
+    alpha1: float = 1.0
+    hbar_omega: float = 1.0
+    repulsive_exponent: int = 12
+    equilibrium_distance: float = 1.0
+    epsilon: float = 0.375
+    core_radius: float = 0.08
+
+    def validate(self) -> None:
+        if self.alpha1 <= 0 or self.hbar_omega <= 0:
+            raise ValueError("positive QDO response parameters are required")
+        if self.repulsive_exponent <= 6:
+            raise ValueError("repulsive exponent must exceed six")
+        if self.equilibrium_distance <= 0 or self.epsilon <= 0:
+            raise ValueError("positive LJ parameters are required")
+        if self.core_radius < 0:
+            raise ValueError("nonnegative core radius required")
+
+
+def qdo_c6(alpha1: float, hbar_omega: float) -> float:
+    return 0.75 * alpha1**2 * hbar_omega
+
+
+def qdo_c9(alpha1: float, hbar_omega: float) -> float:
+    return 3.0 / 16.0 * alpha1**3 * hbar_omega
+
+
+def lj_c6(epsilon: float, equilibrium_distance: float, exponent: int) -> float:
+    return epsilon * exponent * equilibrium_distance**6 / (exponent - 6.0)
+
+
+def lj_repulsive(epsilon: float, re: float, r: float, n: int, m: int = 6) -> float:
+    return epsilon * n * m / (n - m) * (1.0 / n) * (re / r) ** n
+
+
+def lj_attractive(epsilon: float, re: float, r: float, n: int, m: int = 6) -> float:
+    return epsilon * n * m / (n - m) * (1.0 / m) * (re / r) ** m
+
+
+def lj_potential(epsilon: float, re: float, r: float, n: int, m: int = 6) -> float:
+    return lj_repulsive(epsilon, re, r, n, m) - lj_attractive(epsilon, re, r, n, m)
+
+
+def dispersion_energy(coefficient: float, r: float) -> float:
+    return -0.75 * coefficient / r**6
+
+
+def atm_angular_factor(cosines: tuple[float, float, float]) -> float:
+    c1, c2, c3 = cosines
+    return 1.0 + 3.0 * c1 * c2 * c3
+
+
+def atm_energy(
+    coupling: float,
+    cosines: tuple[float, float, float],
+    distances: tuple[float, float, float],
+) -> float:
+    return coupling * atm_angular_factor(cosines) / np.prod(distances) ** 3
+
+
+def triangle_cosines(points: np.ndarray) -> tuple[float, float, float]:
+    if points.shape != (3, points.shape[1]):
+        raise ValueError("three points required")
+    cosines: list[float] = []
+    for center in range(3):
+        left = points[(center + 1) % 3] - points[center]
+        right = points[(center + 2) % 3] - points[center]
+        cosine = float(np.dot(left, right) / (np.linalg.norm(left) * np.linalg.norm(right)))
+        cosines.append(float(np.clip(cosine, -1.0, 1.0)))
+    return tuple(cosines)  # type: ignore[return-value]
+
+
+def trimer_atm_energy(points: np.ndarray, coupling: float) -> float:
+    distances = tuple(
+        float(np.linalg.norm(points[i] - points[j]))
+        for i, j in ((0, 1), (0, 2), (1, 2))
+    )
+    return atm_energy(coupling, triangle_cosines(points), distances)
+
+
+def density_smeared_dispersion(
+    x: np.ndarray,
+    density_a: np.ndarray,
+    density_b: np.ndarray,
+    separation: float,
+    c6: float,
+    core_radius: float,
+) -> float:
+    dx = float(x[1] - x[0])
+    shifted = x[:, None] - (x[None, :] + separation)
+    kernel = 1.0 / (shifted**2 + core_radius**2) ** 3
+    return float(-c6 * np.sum(density_a[:, None] * density_b[None, :] * kernel) * dx**2)
+
+
+def canonical_payload(config: QDOLJATMConfig | None = None) -> dict[str, Any]:
+    cfg = QDOLJATMConfig() if config is None else config
+    return {
+        "schema": SCHEMA,
+        "model_id": "M11",
+        "milestone": MILESTONE,
+        "configuration": asdict(cfg),
+        "study_api": (
+            "openwave.xperiments.m11_cat_ept_soliton_qdo."
+            "qdo_lj_atm_interaction_m113:run_qdo_lj_atm_study"
+        ),
+        "formal_authority": {
+            "repository": "jagg-ix/entropic-physlib-private",
+            "branch": "entropic-physlib-linear-full",
+            "head": FORMAL_HEAD,
+            "sources": list(FORMAL_SOURCES),
+        },
+    }
+
+
+def fingerprint(payload: Mapping[str, Any] | None = None) -> str:
+    selected = canonical_payload() if payload is None else dict(payload)
+    return sha256(_canonical_json(selected).encode()).hexdigest()
+
+
+def run_qdo_lj_atm_study(config: QDOLJATMConfig | None = None) -> dict[str, Any]:
+    cfg = QDOLJATMConfig() if config is None else config
+    cfg.validate()
+    c6_qdo = qdo_c6(cfg.alpha1, cfg.hbar_omega)
+    c9_qdo = qdo_c9(cfg.alpha1, cfg.hbar_omega)
+    c6_lj = lj_c6(cfg.epsilon, cfg.equilibrium_distance, cfg.repulsive_exponent)
+    distances = np.geomspace(1.3, 8.0, 32)
+    split_errors = []
+    tail_errors = []
+    vdw_errors = []
+    tails = []
+    for r in distances:
+        repulsive = lj_repulsive(cfg.epsilon, cfg.equilibrium_distance, float(r), cfg.repulsive_exponent)
+        attractive = lj_attractive(cfg.epsilon, cfg.equilibrium_distance, float(r), cfg.repulsive_exponent)
+        total = lj_potential(cfg.epsilon, cfg.equilibrium_distance, float(r), cfg.repulsive_exponent)
+        split_errors.append(abs(total - (repulsive - attractive)))
+        tail = -attractive
+        tails.append(abs(tail))
+        tail_errors.append(abs(tail + c6_lj / float(r) ** 6))
+        vdw_errors.append(abs(tail - dispersion_energy(4.0 / 3.0 * c6_lj, float(r))))
+    slope = float(np.polyfit(np.log(distances[-12:]), np.log(np.asarray(tails[-12:])), 1)[0])
+    equilateral = np.asarray([[0.0, 0.0], [1.0, 0.0], [0.5, math.sqrt(3.0) / 2.0]])
+    linear = np.asarray([[-1.0, 0.0], [0.0, 0.0], [1.0, 0.0]])
+    permutations = [trimer_atm_energy(equilateral[list(order)], c9_qdo) for order in itertools.permutations(range(3))]
+    eq_factor = atm_angular_factor(triangle_cosines(equilateral))
+    linear_factor = atm_angular_factor(triangle_cosines(linear))
+    scale = 1.7
+    atm_scale_ratio = trimer_atm_energy(equilateral * scale, c9_qdo) / trimer_atm_energy(equilateral, c9_qdo)
+    diagnostics = {
+        "c6_lj": c6_lj,
+        "c6_qdo": c6_qdo,
+        "c6_match_error": abs(c6_lj - c6_qdo),
+        "c9_qdo": c9_qdo,
+        "threebody_invariant_error": abs(cfg.alpha1 * c6_qdo - 4.0 * c9_qdo),
+        "maximum_split_error": max(split_errors),
+        "maximum_tail_error": max(tail_errors),
+        "maximum_vdw_error": max(vdw_errors),
+        "far_field_log_slope": slope,
+        "equilateral_factor": eq_factor,
+        "linear_factor": linear_factor,
+        "atm_permutation_error": float(max(permutations) - min(permutations)),
+        "atm_scaling_error": abs(atm_scale_ratio - scale**-9),
+        "well_depth_error": abs(
+            lj_potential(cfg.epsilon, cfg.equilibrium_distance, cfg.equilibrium_distance, cfg.repulsive_exponent)
+            + cfg.epsilon
+        ),
+    }
+    acceptance = {
+        "lj_split_exact": diagnostics["maximum_split_error"] < 5.0e-13,
+        "tail_is_minus_c6_over_r6": diagnostics["maximum_tail_error"] < 5.0e-13,
+        "tail_equals_vdw_energy": diagnostics["maximum_vdw_error"] < 5.0e-13,
+        "lj_c6_is_qdo_c6": diagnostics["c6_match_error"] < 5.0e-13,
+        "qdo_two_three_body_invariant": diagnostics["threebody_invariant_error"] < 5.0e-13,
+        "tail_slope_is_minus_six": abs(slope + 6.0) < 5.0e-12,
+        "atm_equilateral_repulsive": abs(eq_factor - 11.0 / 8.0) < 5.0e-13,
+        "atm_linear_attractive": abs(linear_factor + 2.0) < 5.0e-13,
+        "atm_permutation_symmetric": diagnostics["atm_permutation_error"] < 5.0e-13,
+        "atm_scales_as_r_minus_nine": diagnostics["atm_scaling_error"] < 5.0e-13,
+        "well_depth_is_minus_epsilon": diagnostics["well_depth_error"] < 5.0e-13,
+    }
+    payload = canonical_payload(cfg)
+    return {
+        **payload,
+        "task": "M11.3",
+        "diagnostics": diagnostics,
+        "acceptance": acceptance,
+        "fingerprint": fingerprint(payload),
+        "passed": all(acceptance.values()),
+    }

openwave/xperiments/m11_cat_ept_soliton_qdo/research/zil/m11_qdo_lj_atm_m113.zc

@@ -0,0 +1,8 @@
+model:m11.3#extends@model:m11.2.
+model:m11.3#pair-interaction@potential:lennard-jones.
+potential:lennard-jones#splits@lean:LennardJonesDispersion.ljPotential_split.
+potential:lennard-jones#tail@equation:minus-c6-over-r6.
+equation:minus-c6-over-r6#identified-with@lean:LennardJonesDispersion.lj_tail_eq_vdw_dispersionEnergy.
+qdo:c6#equals@coefficient:lj-c6.
+qdo:c6#related-to@qdo:c9.
+qdo:c9#controls@potential:axilrod-teller-muto.

tests/test_m11_qdo_lj_atm_m113.py

@@ -0,0 +1,6 @@
+from openwave.xperiments.m11_cat_ept_soliton_qdo.qdo_lj_atm_interaction_m113 import run_qdo_lj_atm_study
+
+def test_qdo_lj_atm_campaign():
+    result=run_qdo_lj_atm_study()
+    assert result['passed'], result
+    assert abs(result['diagnostics']['far_field_log_slope']+6)<1e-10