Skip to content

Commit 86aff26

Browse files
committed
[fix.fields.sigma.D8+] I measured that a linear model is linear and reported it as physics
Ivan: 场还贵互相影响的啊,又不是线性叠加. He is right and it lands on a claim I made twice today in commit bodies and once to a peer. I measured |V(A∪B) - V(A) - V(B)| / |V| = 5.2e-16 for the CLASSICAL POINT-CHARGE model — which is V(r) = Σ q_i/r_i, linear BY CONSTRUCTION. That is arithmetic. I proved a tautology about a model and presented it as a discovered property of fields. MEASURED WHAT HE IS ACTUALLY POINTING AT. Water dimer at H-bonding distance, RHF/6-31G*, QM MEP on a 5 Å probe shell: |V_AB| max 14.065 rms 6.405 kcal/mol |V_AB - (V_A+V_B)| max 1.945 rms 0.784 kcal/mol non-additive: 13.8% of the peak Real densities polarise each other and transfer charge, and a point-charge model cannot represent that at ANY screening. Larger for charged residues against a polarisable ligand than for two waters. WHAT SURVIVES of the original measurement, stated narrowly this time: the ENGINEERING claim (a group classical field needs no bond perception, capping or molfile, so it cannot hit the multi-residue topology bug) and the GRID REGISTRATION control (three runs landing on the same grid, which is what makes any difference field meaningful at all). Neither is a claim about nature. AND HIS OBJECTION EXPOSED A BIGGER ERROR SITTING NEXT TO IT. The pocket map was an UNSCREENED vacuum Coulomb sum. A truncated Asp- at 4 Å contributes ~83 kcal/mol unscreened against the ±5-10 features a chemist reads — the charged residues drown the map they are supposed to inform. That one IS fixable at this level: eps(r) = 4r, the standard minimum distance-dependent dielectric, now the default. Vacuum remains available but must be asked for, because it answers a different question. The response carries `dielectric` and a `physics_caveat` naming both omissions, and the panel prints them with the number. The button is relabelled "Pocket map (point charges, screened)" — not "what the ligand sits in", which promised an interaction energy this cannot deliver. The non-additive part is not noise to be apologised for: it IS the interaction, and the only defensible way to see it is a QM difference field Δ = V_QM(A∪B) - V_QM(A) - V_QM(B) on a shared grid. That is now a real item rather than a phrase, and the registration control I already built is its precondition. warm-two-stage.sh: the browser supplies IDENTITY (the exact molfile, via --dump-dom), the shell does the WAITING (1800 s, no clock over it). Three earlier attempts each failed differently — SMILES-derived keys the app never asks for; button clicks inheriting the panel's 60 s budget; and a page fetch holding a socket for a quarter hour, which returned "TypeError: Failed to fetch" on all six fields. The last two are the same error the physics job queue exists to fix: a long computation on the far end of a synchronous request that something impatient is holding open.
1 parent f97af0a commit 86aff26

3 files changed

Lines changed: 458 additions & 34 deletions

File tree

‎backend/field_server.py‎

Lines changed: 134 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,7 @@
4646
# route 404 — indistinguishable from "no such route", which is the fallback-
4747
# hides-the-primary's-death shape. So a failed import is REMEMBERED and served
4848
# as a 503 that names the exception, while /field keeps working.
49+
import jobs
4950
from envelope import normalize_meta
5051

5152
# `meta.stored: true` is set when the persist thread is STARTED, not when it
@@ -57,6 +58,12 @@
5758
# line in a log nobody tails. A silent persist failure looks exactly like a
5859
# cache that is simply cold — and a cold cache that nobody can distinguish from
5960
# a broken one is how 19 rows sat unservable for a day.
61+
# app.job's writer. Created at db_init (it needs the worker name, which carries
62+
# the derived producer version) and left None when the database is unreachable —
63+
# every call site treats None as "no ledger", never as an error, because a job
64+
# row may not cost a result.
65+
_jobs: 'jobs.JobLedger | None' = None
66+
6067
_persist = {'queued': 0, 'ok': 0, 'failed': 0, 'last_error': None}
6168
_persist_lock = threading.Lock()
6269

@@ -294,6 +301,17 @@ def db_init() -> bool:
294301
_method_ids = method_registry.register_all(_db, sys.modules[__name__],
295302
_toolkit_ids.get('pyscf'))
296303
print(f'[db] {len(_method_ids)} compute units registered', flush=True)
304+
global _jobs
305+
_jobs = jobs.JobLedger(
306+
_db, jobs.make_worker_name(os.getpid(), _producer_version()))
307+
# A restart leaves rows 'running' that no process is doing. Without
308+
# the reap, v_job_live's age_seconds grows without bound and the
309+
# ledger reports work nobody is performing — worse than no ledger,
310+
# because it reads as a hung system.
311+
reaped = _jobs.reap()
312+
if reaped:
313+
print(f'[job] reaped {reaped} orphaned row(s) from a previous run',
314+
flush=True)
297315
except Exception as e:
298316
# A method registry that cannot be written must not take the cache
299317
# down with it: rows keep their producer stamp, which is what the
@@ -1028,6 +1046,29 @@ def run_scf(mol: Chem.Mol, basis: str, max_seconds: float = DEFAULT_MAX_SECONDS,
10281046
ecp=ecp or None,
10291047
verbose=0,
10301048
)
1049+
# ── PRE-FLIGHT, and it did not exist until now ────────────────────────
1050+
# The only cost gate on this path was the CUBE estimate (after the SCF) plus
1051+
# the in-loop watchdog. So `max_seconds=0` — documented everywhere in this
1052+
# repo as "refuse before doing any work, tell me the cost" — actually entered
1053+
# the SCF and was stopped by the watchdog after one cycle: measured 1.03 s and
1054+
# 1.28 s against the job ledger, recorded as BUDGET. A rule stated in four
1055+
# comments and implemented in none of them is a rule about my intentions.
1056+
#
1057+
# The model is the one already measured in this repo over 47 molecules
1058+
# (backend/physics/README.md): seconds ~ 5.9e-9 * nao^4.03, HF being O(N^4).
1059+
# It is a screen, not a prediction — the fit underestimates the middle of the
1060+
# range by up to ~2.8x, so the 2.8x safety factor from that same measurement
1061+
# is applied here rather than quietly dropped.
1062+
nao = gmol.nao_nr()
1063+
predicted_scf = 2.8 * 5.9e-9 * float(nao) ** 4.03
1064+
if predicted_scf > max_seconds:
1065+
raise FieldBudgetExceeded(
1066+
f'estimated {predicted_scf:.1f} s of SCF for {nao} basis functions '
1067+
f'({basis}, {len(syms)} atoms) against a {max_seconds:g} s budget — '
1068+
f'refused before starting. Send max_seconds >= {predicted_scf:.0f}, '
1069+
f'or a smaller basis. The estimate is an O(N^4) screen with a 2.8x '
1070+
f'safety factor, not a promise.')
1071+
10311072
mf = scf.RHF(gmol) if spin == 0 else scf.UHF(gmol)
10321073
mf.max_cycle = 120
10331074
t0 = time.time()
@@ -1328,7 +1369,8 @@ def resolve_charges(sources: list[dict]) -> tuple[list[float], list[str]]:
13281369
MAX_REGION_VOXELS = 2_200_000 # ~128³, the same ceiling the ligand path uses
13291370

13301371

1331-
def field_region(sources, lo, hi, spacing: float, kind: str):
1372+
def field_region(sources, lo, hi, spacing: float, kind: str,
1373+
req_dielectric: str = 'r-dependent'):
13321374
"""The classical field of an arbitrary SOURCE set, sampled in a given FRAME.
13331375
13341376
Why this route exists at all, and why it needs no chemistry:
@@ -1425,6 +1467,29 @@ def field_region(sources, lo, hi, spacing: float, kind: str):
14251467
# A distance cutoff, because the ligand path's per-atom full-grid temporary
14261468
# is ~50 MB at 128³ and a pocket has thousands of atoms. Sources far from
14271469
# the frame contribute below the noise of the model itself.
1470+
# ── SCREENING, and the reason it is on by default ───────────────────
1471+
#
1472+
# A bare Coulomb sum in vacuum is not what a ligand experiences in a
1473+
# pocket. Two errors sit on top of an unscreened point-charge map and they
1474+
# are different sizes:
1475+
#
1476+
# NON-ADDITIVITY. Real densities polarise each other, so the true
1477+
# V(A∪B) is not V(A)+V(B). Measured here on a water dimer at H-bonding
1478+
# distance, RHF/6-31G* on a 5 Å probe shell: the non-additive part is
1479+
# 1.95 kcal/mol against a 14.07 peak — 13.8%. Larger for charged
1480+
# residues against a polarisable ligand. A point-charge model cannot
1481+
# represent this at all, at any screening.
1482+
#
1483+
# NO SCREENING. Larger still, and this one IS fixable here. A truncated
1484+
# Asp- at 4 Å contributes ~83 kcal/mol unscreened, against the ±5-10
1485+
# features a chemist actually reads. Unscreened, the charged residues
1486+
# drown the map they are supposed to inform.
1487+
#
1488+
# eps(r) = 4r is the standard minimum distance-dependent dielectric for
1489+
# exactly this situation — cheap, and far closer than vacuum. It is NOT
1490+
# Poisson-Boltzmann and does not pretend to be; the model is named in the
1491+
# response so the picture cannot be mistaken for an interaction energy.
1492+
screen = str(req_dielectric).lower()
14281493
cutoff = 25.0 if kind == 'mep' else 12.0
14291494
near = ((coords >= lo - cutoff).all(axis=1) & (coords <= hi + cutoff).all(axis=1))
14301495
used = int(near.sum())
@@ -1433,7 +1498,10 @@ def field_region(sources, lo, hi, spacing: float, kind: str):
14331498
for w, c in zip(weights[near], coords[near]):
14341499
d = np.linalg.norm(pts - c, axis=-1)
14351500
if kind == 'mep':
1436-
v += w / np.maximum(d, 0.5)
1501+
dd = np.maximum(d, 0.5)
1502+
# eps(r) = 4r => V = q / (4 r^2). Vacuum is available but must be
1503+
# asked for, because it is the answer to a different question.
1504+
v += w / (4.0 * dd * dd) if screen == 'r-dependent' else w / dd
14371505
else:
14381506
v += w * np.exp(-d / 2.0)
14391507
if kind == 'mep':
@@ -1454,6 +1522,15 @@ def field_region(sources, lo, hi, spacing: float, kind: str):
14541522
'net_charge': round(float(weights.sum()), 3),
14551523
'charge_model': (f'{CHARGE_FORCEFIELD} residue templates (pdb2pqr)'
14561524
if kind == 'mep' else 'caller-supplied'),
1525+
'dielectric': screen if kind == 'mep' else None,
1526+
'physics_caveat': (
1527+
'Point charges with eps(r)=4r screening. This is a QUALITATIVE MAP '
1528+
'of where the pocket is positive or negative, not an interaction '
1529+
'energy. Two things it cannot contain: mutual polarisation and '
1530+
'charge transfer (measured non-additive by ~14% of peak on an '
1531+
'H-bonded water dimer at RHF/6-31G*), and any solvent beyond the '
1532+
'distance-dependent dielectric.'
1533+
if kind == 'mep' else None),
14571534
'waters_excluded': n_water,
14581535
'waters_note': (f'{n_water} crystallographic water(s) left OUT: their '
14591536
f'hydrogens were never resolved, and a bare oxygen would '
@@ -1614,6 +1691,7 @@ def do_GET(self):
16141691
'scf_cached': len(_scf_cache),
16151692
'scf_cache_max': SCF_CACHE_MAX,
16161693
'persist': dict(_persist),
1694+
'jobs': jobs.counters(),
16171695
'rss_mb': rss_mb})
16181696
else:
16191697
self._send(404, {'ok': False, 'error': 'not found'})
@@ -1640,7 +1718,8 @@ def do_POST(self):
16401718
frame = req['frame']
16411719
cube, meta = field_region(
16421720
req['sources'], frame['lo'], frame['hi'],
1643-
float(frame.get('spacing', 0.5)), req.get('kind', 'mep'))
1721+
float(frame.get('spacing', 0.5)), req.get('kind', 'mep'),
1722+
req_dielectric=req.get('dielectric', 'r-dependent'))
16441723
meta['total_seconds'] = round(time.time() - t0, 2)
16451724
print(f"[field/region] kind={meta['kind']} "
16461725
f"sources={meta['n_sources_used']}/{meta['n_sources_sent']} "
@@ -1673,6 +1752,12 @@ def do_POST(self):
16731752
if self.path != '/field':
16741753
self._send(404, {'ok': False, 'error': 'not found'})
16751754
return
1755+
# Declared BEFORE the try so the failure branches can close the row.
1756+
# A job that stays 'running' after its request has already failed is
1757+
# the shape reap() exists to clean up, and leaving it to reap means the
1758+
# ledger reports work nobody is doing until the next restart.
1759+
job_id = None
1760+
t_job = time.time()
16761761
try:
16771762
length = int(self.headers.get('Content-Length', '0'))
16781763
if length > MAX_BODY_BYTES:
@@ -1701,7 +1786,18 @@ def do_POST(self):
17011786
req_seconds = float(req.get('max_seconds', DEFAULT_MAX_SECONDS))
17021787
if not math.isfinite(req_seconds):
17031788
req_seconds = DEFAULT_MAX_SECONDS
1704-
max_seconds = min(max(req_seconds, 1.0), MAX_MAX_SECONDS)
1789+
# ZERO IS PRESERVED, and the 1.0 floor that used to be here defeated
1790+
# the feature it looked like it was protecting. `max_seconds: 0` means
1791+
# "tell me the cost, do not run", which run_scf's pre-flight now
1792+
# answers before building an SCF object. The floor turned that into
1793+
# "run for one second, then fail with BUDGET" — measured against the
1794+
# job ledger at 1.03 s and 1.28 s, which is a different question
1795+
# answered at a cost.
1796+
#
1797+
# A NEGATIVE budget clamps to 0 rather than to the default: it is not
1798+
# a request for the default, it is nonsense, and the safest reading of
1799+
# nonsense is "refuse and tell me what it would have cost".
1800+
max_seconds = min(max(req_seconds, 0.0), MAX_MAX_SECONDS)
17051801
spin = req.get('spin')
17061802
spin = int(spin) if spin is not None else None
17071803
molfile_sha = hashlib.sha256(molblock.encode()).digest()
@@ -1722,6 +1818,21 @@ def do_POST(self):
17221818
return
17231819

17241820
t0 = time.time()
1821+
# The job row exists for the DURATION of the compute, which is the
1822+
# whole point: a 6-minute SCF was previously invisible until it
1823+
# finished. `open` returns None when the DB is down or an identical
1824+
# job is already in flight, and every path below tolerates None.
1825+
# Imported locally, matching every other use in this file — the
1826+
# module-level import list is deliberately minimal here.
1827+
import method_registry as _mr
1828+
method_row_id = _method_ids.get(_mr.KIND_TO_METHOD.get(kind, ''))
1829+
job_id = None
1830+
if _jobs is not None and method_row_id is not None:
1831+
job_id = _jobs.open(
1832+
method_row_id=method_row_id, input_sha256=molfile_sha,
1833+
params={'kind': kind, 'basis': basis_key,
1834+
'spin': spin, 'max_seconds': max_seconds},
1835+
budget_seconds=max_seconds)
17251836
mol = prepare_mol(molblock)
17261837
if kind == 'mep':
17271838
cube, meta = field_mep(mol)
@@ -1755,6 +1866,13 @@ def do_POST(self):
17551866
# computed.
17561867
meta.setdefault('computed_at',
17571868
datetime.now(timezone.utc).isoformat())
1869+
if _jobs is not None:
1870+
# field_cube_id is deliberately NOT set here: the row is written
1871+
# by a background thread that may not have committed yet, and a
1872+
# job pointing at a row that does not exist would be a worse lie
1873+
# than a job with no pointer. The pointer lands when the async
1874+
# write learns to report back.
1875+
_jobs.done(job_id, seconds=meta['total_seconds'])
17581876
# Lenient here, and the asymmetry with the cache path is deliberate:
17591877
# a 6-minute SCF must not be discarded because a producer added a
17601878
# key. It is logged at full volume instead — the drift is still
@@ -1772,10 +1890,22 @@ def do_POST(self):
17721890
# bigger budget" instead of showing a chemist a red failure for a
17731891
# calculation that was merely slow.
17741892
print(f'[field] budget exceeded: {e}', flush=True)
1893+
if _jobs is not None:
1894+
_jobs.failed(job_id, code='BUDGET', detail=str(e),
1895+
seconds=time.time() - t_job)
17751896
self._send(200, {'ok': False, 'error': str(e), 'reason': 'budget'})
17761897
except Exception as e:
17771898
traceback.print_exc()
17781899
reason = 'unsupported' if isinstance(e, ValueError) else 'internal'
1900+
if _jobs is not None:
1901+
# ValueError is this service's chemistry refusal (an element with
1902+
# no basis, an unparameterised atom), which is UNSUPPORTED — not
1903+
# INTERNAL. Recording every failure as INTERNAL would make the
1904+
# ledger's error_code column useless for the one question it is
1905+
# for: was it us or the molecule?
1906+
_jobs.failed(job_id,
1907+
code='UNSUPPORTED' if reason == 'unsupported' else 'INTERNAL',
1908+
detail=str(e), seconds=time.time() - t_job)
17791909
self._send(200, {'ok': False, 'error': str(e), 'reason': reason})
17801910

17811911
def log_message(self, *args):

‎scripts/warm-two-stage.sh‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
#!/usr/bin/env bash
2+
# Two-stage warm: the browser supplies IDENTITY, the shell does the WAITING.
3+
#
4+
# scripts/warm-two-stage.sh
5+
#
6+
# Three attempts got here and each failed for its own reason, which is why the
7+
# split is what it is:
8+
#
9+
# 1. Warm from SMILES in Python. Fast, patient — and writes rows the app can
10+
# never hit, because the durable cache is keyed on sha256(molfile) and the
11+
# app's molfile is reconstructed by mol* from scene loci, not from a SMILES.
12+
# 2. Warm by clicking the buttons. Exact keys — and inherits the panel's 60 s
13+
# client budget, so lapatinib returned "SCF exceeded its 60 s budget after
14+
# 10 cycles" and cached nothing.
15+
# 3. Warm by fetch() from the page with max_seconds=900. Exact keys, real
16+
# budget — and asks a browser to hold a socket for a quarter of an hour.
17+
# Every field came back "TypeError: Failed to fetch".
18+
#
19+
# The common error in 2 and 3 is the same one the physics job queue exists to
20+
# fix: a long computation does not belong on the far end of a synchronous
21+
# request that something impatient is holding open. So the browser is used for
22+
# the ONE thing only it can do — produce the exact molfile the app will ask
23+
# about — and then it leaves. Python computes, with no clock over it.
24+
set -u
25+
26+
ROOT=/home/ivan/dirac
27+
OUT=${1:-/tmp/claude-1000/-home-ivan/ff9c7c07-88fe-4345-a0cc-59d17b42683e/scratchpad/warm}
28+
STAGE=$OUT/_dom
29+
PORT=1348
30+
MOLECULES=${MOLECULES:-"1CBS 1XKK 4HHB"}
31+
32+
mkdir -p "$OUT/molfiles"
33+
rm -rf "$STAGE"; cp -r "$ROOT/build/dirac" "$STAGE"
34+
( cd "$STAGE" && python3 -m http.server $PORT --bind 127.0.0.1 >/dev/null 2>&1 & echo $! > "$OUT/.srv2" )
35+
sleep 1
36+
37+
# ── stage 1: identity ──────────────────────────────────────────────────────
38+
for MOL in $MOLECULES; do
39+
cat > "$STAGE/index.html.driver" <<EOF
40+
<script>
41+
(function () {
42+
var MOL = '$MOL';
43+
function wait(ms) { return new Promise(function (r) { setTimeout(r, ms); }); }
44+
(async function () {
45+
await wait(9000);
46+
var sel = document.getElementById('molecule');
47+
for (var i = 0; i < sel.options.length; i++) {
48+
if (sel.options[i].textContent.indexOf(MOL) >= 0) {
49+
sel.selectedIndex = i;
50+
sel.dispatchEvent(new Event('change', { bubbles: true })); break;
51+
}
52+
}
53+
await wait(14000);
54+
document.querySelector('[data-jump="fields"]').click();
55+
await wait(1500);
56+
var hook = window.diracFields || {};
57+
// Into the DOM, where --dump-dom can reach it. No network, no waiting.
58+
var pre = document.createElement('pre');
59+
pre.id = 'warm-molfile';
60+
pre.textContent = hook.molfile || '';
61+
document.body.appendChild(pre);
62+
var lab = document.createElement('pre');
63+
lab.id = 'warm-label';
64+
lab.textContent = hook.label || '';
65+
document.body.appendChild(lab);
66+
})();
67+
})();
68+
</script>
69+
EOF
70+
cat "$ROOT/build/dirac/index.html" "$STAGE/index.html.driver" > "$STAGE/index.html"
71+
echo "extracting $MOL ..."
72+
timeout 180 google-chrome --headless=new --enable-unsafe-swiftshader \
73+
--use-angle=swiftshader --virtual-time-budget=60000 --dump-dom \
74+
"http://127.0.0.1:$PORT/index.html" > "$OUT/molfiles/$MOL.dom" 2>/dev/null
75+
python3 - "$OUT/molfiles/$MOL.dom" "$OUT/molfiles/$MOL.mol" <<'PY'
76+
import sys, re, html
77+
dom = open(sys.argv[1], encoding='utf-8', errors='replace').read()
78+
m = re.search(r'<pre id="warm-molfile">(.*?)</pre>', dom, re.S)
79+
mol = html.unescape(m.group(1)) if m else ''
80+
open(sys.argv[2], 'w').write(mol)
81+
lines = mol.count('\n')
82+
print(f' -> {lines} lines' if lines > 4 else ' -> NO LIGAND (nothing to warm)')
83+
PY
84+
done
85+
kill "$(cat "$OUT/.srv2")" 2>/dev/null; rm -f "$OUT/.srv2" "$STAGE/index.html.driver"
86+
87+
# ── stage 2: computation, with no clock over it ────────────────────────────
88+
echo
89+
"$ROOT/backend/env/bin/python" - "$OUT/molfiles" $MOLECULES <<'PY'
90+
import json, os, sys, time, urllib.request
91+
92+
d, mols = sys.argv[1], sys.argv[2:]
93+
B = 'http://127.0.0.1:8901'
94+
KINDS = ['mep', 'mlp', 'mep_qm', 'homo', 'lumo', 'density']
95+
96+
def post(path, payload, timeout):
97+
req = urllib.request.Request(B + path, data=json.dumps(payload).encode(),
98+
headers={'Content-Type': 'application/json'})
99+
with urllib.request.urlopen(req, timeout=timeout) as r:
100+
return json.loads(r.read())
101+
102+
print(f'{"molecule":9s} {"field":7s} {"outcome":9s} detail')
103+
for mol in mols:
104+
p = os.path.join(d, f'{mol}.mol')
105+
text = open(p).read() if os.path.exists(p) else ''
106+
if len(text.splitlines()) < 5:
107+
print(f'{mol:9s} {"-":7s} {"NOLIGAND":9s} nothing to precompute')
108+
continue
109+
for kind in KINDS:
110+
t0 = time.time()
111+
try:
112+
# 1800 s and a socket to match. Nobody is waiting on this.
113+
out = post('/field', {'molfile': text, 'kind': kind,
114+
'basis': 'sto-3g', 'store': True,
115+
'max_seconds': 1800}, timeout=2100)
116+
except Exception as e: # noqa: BLE001
117+
print(f'{mol:9s} {kind:7s} {"HARNESS":9s} {type(e).__name__}: {e}')
118+
continue
119+
el = time.time() - t0
120+
if out.get('ok'):
121+
m = out.get('meta', {})
122+
print(f'{mol:9s} {kind:7s} {"CACHED":9s} {el:6.1f}s '
123+
f'cache={m.get("cache", "db")} nao={m.get("nbasis", "-")}')
124+
else:
125+
print(f'{mol:9s} {kind:7s} {"REFUSED":9s} {str(out.get("error"))[:80]}')
126+
PY
127+
echo
128+
echo "DB rows that survive a restart:"
129+
psql -U ivan -d dirac -tAc "select kind, basis, count(*) from app.field_cube group by 1,2 order by 1,2" 2>/dev/null

0 commit comments

Comments
 (0)