Skip to content

Commit da29404

Browse files
committed
[fix.fields.sigma.D8+] the crate was never an isovalue problem, it was a box problem
Third attempt, and the first one that is about the molecule. A clean-context computational chemist took the second one apart and the measurement confirms them; this is my own change from this morning being retracted, for the second time today, on better evidence. 1. iso = 0.05 MLP, hand-tuned on aspirin. Measured: clipped the grid wall on all four molecules tried, INCLUDING aspirin. Rendered as a gold crate. 2. iso = 97th percentile of |field|. Never clips. Also a function of `pad`, a RENDERING constant — measured on aspirin, pad 3->10 A moves it 25.5 -> 3.4 kcal/mol, a 7.6x swing with the physics IDENTICAL. I replaced an arbitrary constant with a HIDDEN one carrying more leverage and printed it as a measurement. It also selects the highest-|V| voxels, which are always interior: 84-94% of the enclosed volume sits inside the vdW surface, where nothing binds and where the softened Coulomb has the wrong sign anyway. 3. This. The contour is FIXED in physical units and the BOX adapts. The reframe is the finding: an adaptive RULER destroys comparison, an adaptive BOX costs nothing. When a picture clips, widen the frame — do not move the ruler and call the result measured. Verified after the change: aspirin / retinoic acid / benzene mep iso 10.00 pad 4.0 closes same three mlp iso 0.25 pad 4.0 closes acetate -1, choline +1 mep iso 10.00 pad 12.0 DOES NOT CLOSE That last row is the part worth keeping. A charged ligand's monopole does not decay fast enough to close inside any usable box, so the honest answer is not a cleverer isovalue — it is to say so. `contour_closes_in_box: false` and the wall value ride on the response, and the panel prints "open surface (field is 27.761 at the box edge)". Three-valued, as ever: closed, open, or refused — never a flat face drawn as if it were chemistry. Two more from the same review, both measured before I believed them: - MEP and QM MEP were contoured in their own bare units, kcal/mol against Ha/e, so toggling between two renderings of the SAME physical quantity moved the contour from 18 to 31.4 kcal/mol without saying so. One unit on screen now, converted at the isovalue by cubeScale. - MLP is one-signed for most drug-like molecules (aspirin 99.9% positive, ibuprofen 100%). Calling it "diverging" and drawing a negative lobe that does not exist is a claim about the molecule. `single_signed` is measured per response and the readout says "no negative lobe exists" instead of implying one is hiding below the contour. The isoLock/re-fit machinery from three commits ago is deleted. It was the right instinct — hold the ruler still — applied to a ruler that should never have been derived from the box in the first place. A fixed contour is comparable by construction and needs no lock. NEXT: the same review measured that /surface/mep_at has no cost gate and no basis whitelist, that the physics daemon's SCF has no watchdog at all (the 36-minute fix landed on one of two daemons), and that STO-3G inverts the HOMO ordering of a substituted-benzene series — nitrobenzene comes out more electron-rich than benzene. Those are next.
1 parent 8f32bac commit da29404

2 files changed

Lines changed: 127 additions & 96 deletions

File tree

‎backend/field_server.py‎

Lines changed: 72 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -127,12 +127,14 @@
127127
# this version is re-registered with different source — a forgotten bump is a
128128
# loud startup error, never a silently stale cache (design: migration 006).
129129
PRODUCER_SERVICE = 'dirac-fields'
130-
PRODUCER_VERSION = '1.8'
130+
PRODUCER_VERSION = '1.9'
131131
PRODUCER_NOTES = ('wall-clock deadline inside the SCF loop + measured cube-cost '
132132
'refusal; open-shell d/f metals refused without an explicit '
133133
'spin (group 12 exempt); salt stripping refuses to discard a '
134134
'coordination metal; iso_suggest = 97th percentile of |field| '
135-
'for mep/mlp; _eri dropped and the SCF cache bounded to 6 LRU; '
135+
'for mep/mlp RETIRED — it was a measurement of the padding (7.6x '
136+
'swing on aspirin); fixed physical contour + adaptive box instead; '
137+
'_eri dropped and the SCF cache bounded to 6 LRU; '
136138
'molblock refusals name the defect')
137139

138140

@@ -595,32 +597,47 @@ def field_mep(mol: Chem.Mol, spacing=0.4, pad=4.0):
595597
raise ValueError('all Gasteiger charges are zero — classical MEP would '
596598
'be an empty picture; use mep_qm')
597599

598-
lo = coords.min(axis=0) - pad
599-
hi = coords.max(axis=0) + pad
600-
dims = np.maximum(np.ceil((hi - lo) / spacing).astype(int) + 1, 8)
601-
dims = np.minimum(dims, GRID_MAX_DIM) # hard cap: 128³ ≈ 2M voxels
602-
axes = np.diag((hi - lo) / (dims - 1))
603-
604-
xs = np.linspace(lo[0], hi[0], dims[0])
605-
ys = np.linspace(lo[1], hi[1], dims[1])
606-
zs = np.linspace(lo[2], hi[2], dims[2])
607-
gx, gy, gz = np.meshgrid(xs, ys, zs, indexing='ij')
608-
pts = np.stack([gx, gy, gz], axis=-1) # (nx,ny,nz,3)
609-
610-
v = np.zeros(tuple(dims), dtype=float)
611-
for q, c in zip(charges, coords):
612-
r = np.linalg.norm(pts - c, axis=-1)
613-
v += q / np.maximum(r, 0.5)
614-
v *= 332.06 # kcal/mol per (e²/Å)
600+
iso = FIXED_ISO['mep']
601+
602+
def evaluate(pad_a: float):
603+
lo = coords.min(axis=0) - pad_a
604+
hi = coords.max(axis=0) + pad_a
605+
dims = np.maximum(np.ceil((hi - lo) / spacing).astype(int) + 1, 8)
606+
dims = np.minimum(dims, GRID_MAX_DIM)
607+
xs = np.linspace(lo[0], hi[0], dims[0])
608+
ys = np.linspace(lo[1], hi[1], dims[1])
609+
zs = np.linspace(lo[2], hi[2], dims[2])
610+
gx, gy, gz = np.meshgrid(xs, ys, zs, indexing='ij')
611+
pts = np.stack([gx, gy, gz], axis=-1)
612+
vv = np.zeros(tuple(dims), dtype=float)
613+
for q, c in zip(charges, coords):
614+
r = np.linalg.norm(pts - c, axis=-1)
615+
vv += q / np.maximum(r, 0.5)
616+
return lo, hi, dims, vv * 332.06 # kcal/mol per (e²/Å)
617+
618+
# Grow the BOX, not the ruler, until the contour closes inside it.
619+
pad_used = pad
620+
lo, hi, dims, v = evaluate(pad_used)
621+
while wall_max(v) >= iso and pad_used < PAD_MAX:
622+
pad_used = min(pad_used + PAD_STEP, PAD_MAX)
623+
lo, hi, dims, v = evaluate(pad_used)
624+
clear = wall_max(v) < iso
615625

626+
axes = np.diag((hi - lo) / (dims - 1))
616627
cube = write_cube(lo, axes, dims, v, syms, coords,
617628
'kind=mep units=kcal/mol charges=gasteiger')
618629
meta = {
619630
'kind': 'mep', 'units': 'kcal/mol', 'charges': 'gasteiger',
620631
'net_charge': int(round(charges.sum())),
621632
'dims': dims.tolist(), **grid_spacing_meta(lo, hi, dims, spacing),
622633
'vmin': float(v.min()), 'vmax': float(v.max()),
623-
'iso_suggest': suggest_iso(v),
634+
'iso_fixed': iso,
635+
'pad_used_angstrom': round(pad_used, 1),
636+
'wall_max': round(wall_max(v), 3),
637+
# False means the surface still runs off the edge of the grid and is
638+
# being drawn as a flat face. Stated, not hidden: a charged ligand's
639+
# monopole does not decay fast enough to close inside any usable box.
640+
'contour_closes_in_box': bool(clear),
624641
}
625642
return cube, meta
626643

@@ -641,30 +658,49 @@ def field_mlp(mol: Chem.Mol, spacing=0.4, pad=4.0):
641658
raise ValueError('Crippen contributions undefined for this molecule')
642659
syms, coords = mol_atoms(mol)
643660

644-
lo = coords.min(axis=0) - pad
645-
hi = coords.max(axis=0) + pad
646-
dims = np.maximum(np.ceil((hi - lo) / spacing).astype(int) + 1, 8)
647-
dims = np.minimum(dims, GRID_MAX_DIM)
648-
axes = np.diag((hi - lo) / (dims - 1))
649-
xs = np.linspace(lo[0], hi[0], dims[0])
650-
ys = np.linspace(lo[1], hi[1], dims[1])
651-
zs = np.linspace(lo[2], hi[2], dims[2])
652-
gx, gy, gz = np.meshgrid(xs, ys, zs, indexing='ij')
653-
pts = np.stack([gx, gy, gz], axis=-1)
654-
655-
v = np.zeros(tuple(dims), dtype=float)
656-
for fi, c in zip(f, coords):
657-
d = np.linalg.norm(pts - c, axis=-1)
658-
v += fi * np.exp(-d / 2.0)
661+
iso = FIXED_ISO['mlp']
662+
663+
def evaluate(pad_a: float):
664+
lo = coords.min(axis=0) - pad_a
665+
hi = coords.max(axis=0) + pad_a
666+
dims = np.maximum(np.ceil((hi - lo) / spacing).astype(int) + 1, 8)
667+
dims = np.minimum(dims, GRID_MAX_DIM)
668+
xs = np.linspace(lo[0], hi[0], dims[0])
669+
ys = np.linspace(lo[1], hi[1], dims[1])
670+
zs = np.linspace(lo[2], hi[2], dims[2])
671+
gx, gy, gz = np.meshgrid(xs, ys, zs, indexing='ij')
672+
pts = np.stack([gx, gy, gz], axis=-1)
673+
vv = np.zeros(tuple(dims), dtype=float)
674+
for fi, c in zip(f, coords):
675+
d = np.linalg.norm(pts - c, axis=-1)
676+
vv += fi * np.exp(-d / 2.0)
677+
return lo, hi, dims, vv
678+
679+
pad_used = pad
680+
lo, hi, dims, v = evaluate(pad_used)
681+
while wall_max(v) >= iso and pad_used < PAD_MAX:
682+
pad_used = min(pad_used + PAD_STEP, PAD_MAX)
683+
lo, hi, dims, v = evaluate(pad_used)
684+
clear = wall_max(v) < iso
659685

686+
axes = np.diag((hi - lo) / (dims - 1))
660687
cube = write_cube(lo, axes, dims, v, syms, coords,
661688
'kind=mlp units=logP kernel=fauchere_exp_d_over_2 charges=crippen')
662689
meta = {
663690
'kind': 'mlp', 'units': 'MLP (Crippen/Fauchère)', 'method': 'crippen',
664691
'total_logp': float(f.sum()),
665692
'dims': dims.tolist(), **grid_spacing_meta(lo, hi, dims, spacing),
666693
'vmin': float(v.min()), 'vmax': float(v.max()),
667-
'iso_suggest': suggest_iso(v),
694+
'iso_fixed': iso,
695+
'pad_used_angstrom': round(pad_used, 1),
696+
'wall_max': round(wall_max(v), 4),
697+
'contour_closes_in_box': bool(clear),
698+
# The Fauchere kernel decays over ~2 A, so every atom contributes
699+
# everywhere and the field is one-signed for most drug-like molecules.
700+
# Measured: aspirin is 99.9% positive, ibuprofen 100%. Calling it
701+
# "diverging" and drawing a negative lobe that does not exist is a
702+
# claim about the molecule; this states the truth instead.
703+
'single_signed': bool(v.min() >= 0 or v.max() <= 0),
668704
}
669705
return cube, meta
670706

‎src/app.frontend.facets.molstar-rdkit.editable/facets/field-wells/index.ts‎

Lines changed: 55 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -107,7 +107,16 @@ export type FieldKind = 'mep' | 'mep_qm' | 'homo' | 'lumo' | 'density' | 'mlp';
107107

108108
interface KindSpec {
109109
label: string;
110+
/** The FIXED contour, in `unit`. Identical across molecules by design. */
110111
iso: number;
112+
/**
113+
* display unit -> cube unit. The QM potential's cube is in Ha/e while the
114+
* classical one is in kcal/mol, and both were being contoured with their
115+
* own bare number — so toggling between two renderings of the SAME
116+
* physical quantity silently moved the contour from 18 to 31.4 kcal/mol.
117+
* One unit on screen, converted here.
118+
*/
119+
cubeScale: number;
111120
diverging: boolean;
112121
unit: string;
113122
/** Theme token names; the numbers beside them are the fallback only. */
@@ -125,14 +134,14 @@ interface KindSpec {
125134
// only chroma was cut. Divergent pairs verified ≥0.10 ΔE apart by
126135
// design/check_palette.py — run it after touching any value here.
127136
const Kinds: Record<FieldKind, KindSpec> = {
128-
mep: { label: 'Electrostatic well', iso: 8, diverging: true, unit: 'kcal/mol', posToken: '--viz-mep-pos', negToken: '--viz-mep-neg', posColor: 0x6788bc, negColor: 0xbd777b, quantum: false },
129-
mep_qm: { label: 'QM potential well', iso: 0.05, diverging: true, unit: 'Ha/e', posToken: '--viz-mep-pos', negToken: '--viz-mep-neg', posColor: 0x6788bc, negColor: 0xbd777b, quantum: true },
130-
homo: { label: 'HOMO', iso: 0.04, diverging: true, unit: 'amp', posToken: '--viz-orb-pos', negToken: '--viz-orb-neg', posColor: 0x7fc7a5, negColor: 0xa397d3, quantum: true },
131-
lumo: { label: 'LUMO', iso: 0.04, diverging: true, unit: 'amp', posToken: '--viz-orb-pos', negToken: '--viz-orb-neg', posColor: 0x7fc7a5, negColor: 0xa397d3, quantum: true },
132-
density: { label: 'e⁻ density', iso: 0.05, diverging: false, unit: 'e/Bohr³', posToken: '--viz-density', negToken: '--viz-density', posColor: 0xd8aa75, negColor: 0xd8aa75, quantum: true },
137+
mep: { label: 'Electrostatic well', iso: 10, cubeScale: 1, diverging: true, unit: 'kcal/mol', posToken: '--viz-mep-pos', negToken: '--viz-mep-neg', posColor: 0x6788bc, negColor: 0xbd777b, quantum: false },
138+
mep_qm: { label: 'QM potential well', iso: 10, cubeScale: 1 / 627.5094740631, diverging: true, unit: 'kcal/mol', posToken: '--viz-mep-pos', negToken: '--viz-mep-neg', posColor: 0x6788bc, negColor: 0xbd777b, quantum: true },
139+
homo: { label: 'HOMO', iso: 0.04, cubeScale: 1, diverging: true, unit: 'amp', posToken: '--viz-orb-pos', negToken: '--viz-orb-neg', posColor: 0x7fc7a5, negColor: 0xa397d3, quantum: true },
140+
lumo: { label: 'LUMO', iso: 0.04, cubeScale: 1, diverging: true, unit: 'amp', posToken: '--viz-orb-pos', negToken: '--viz-orb-neg', posColor: 0x7fc7a5, negColor: 0xa397d3, quantum: true },
141+
density: { label: 'e⁻ density', iso: 0.05, cubeScale: 1, diverging: false, unit: 'e/Bohr³', posToken: '--viz-density', negToken: '--viz-density', posColor: 0xd8aa75, negColor: 0xd8aa75, quantum: true },
133142
// Default iso must sit BELOW the hydrophilic side's typical |min| (~0.06
134143
// on aspirin) or the cyan lobes never exist and the field looks all-grease.
135-
mlp: { label: 'Lipophilicity', iso: 0.05, diverging: true, unit: 'MLP', posToken: '--viz-mlp-pos', negToken: '--viz-mlp-neg', posColor: 0xd5b979, negColor: 0x74ccdd, quantum: false },
144+
mlp: { label: 'Lipophilicity', iso: 0.25, cubeScale: 1, diverging: true, unit: 'MLP', posToken: '--viz-mlp-pos', negToken: '--viz-mlp-neg', posColor: 0xd5b979, negColor: 0x74ccdd, quantum: false },
136145
};
137146

138147
interface FieldMeta {
@@ -149,15 +158,22 @@ interface FieldMeta {
149158
units?: string;
150159
natoms?: number;
151160
nbasis?: number;
152-
/** Measured from the cube for MEP/MLP, whose scale is molecule-dependent. */
153-
iso_suggest?: number;
161+
/** The fixed contour the backend grew the box around. */
162+
iso_fixed?: number;
163+
pad_used_angstrom?: number;
164+
wall_max?: number;
165+
/** False = the surface runs off the grid and is drawn as a flat face. */
166+
contour_closes_in_box?: boolean;
167+
/** MLP is one-signed for most drug-like molecules; measured, not assumed. */
168+
single_signed?: boolean;
154169
}
155170

156171
let plugin: PluginContext | null = null;
157172
let molfile: string | null = null;
158173
let ligandLabel: string | null = null;
159174
let activeKind: FieldKind | null = null;
160175
let activeVolume: Volume | null = null;
176+
let lastMeta: FieldMeta | null = null;
161177
let busy = false;
162178

163179
/**
@@ -386,61 +402,51 @@ function isoMultiplier(): number {
386402
* The isovalue the slider's 1.0x position means.
387403
*
388404
* MEP and MLP are sums of atomic contributions, so their scale belongs to the
389-
* molecule and the backend measures it from the cube (`iso_suggest`, the ~97th
390-
* percentile of |field|). A constant there produced a surface that enclosed
391-
* the whole padded grid and clipped against the box wall — on retinoic acid
392-
* the lipophilicity field rendered as a gold crate with the ligand inside it.
393-
*
394-
* Orbitals and density keep their constants on purpose: a wavefunction is
395-
* normalised, which is exactly why the conventional values work for them.
405+
* molecule. Both earlier answers were wrong in opposite directions: a
406+
* hand-tuned constant clipped the grid wall on every molecule tried, and the
407+
* per-molecule percentile that replaced it was a function of the PADDING —
408+
* 7.6x swing on aspirin from pad alone. The contour is fixed now and the
409+
* BACKEND grows the box until it closes.
396410
*/
397-
let isoBase: number | null = null;
398-
399411
/**
400-
* The isovalue, LOCKED across molecules once it has been fitted.
412+
* The contour, in the field's own display units.
401413
*
402-
* Fitting per molecule is automatic gain control: two analogues render equally
403-
* hot even when one carries three times the potential, and the panel is built
404-
* to invite exactly that comparison. A fixed constant is not the alternative —
405-
* measured, the old one clipped the grid wall on every molecule tried, which
406-
* is how the lipophilicity field came to render as a box.
414+
* It is a CONSTANT, and that is the point. Two earlier attempts moved it:
415+
* a hand-tuned value that clipped the grid wall on every molecule tried, then
416+
* a per-molecule percentile that never clipped and was a function of the
417+
* PADDING — measured, 7.6x swing on aspirin from pad alone, with the physics
418+
* identical. A ruler derived from the box is a measurement of the box.
407419
*
408-
* So: MEASURE it once, then HOLD it. The first molecule of a session fits the
409-
* scale from its own field; every molecule after that is drawn on the same
410-
* absolute surface and is therefore comparable. Re-fit is a deliberate act
411-
* with a button, and the readout always says which regime it is in — a scale
412-
* that changed silently underneath a comparison is the failure being fixed.
420+
* The crate was a BOX problem, so the backend grows the box until the contour
421+
* closes inside it and reports whether it succeeded. The ruler holds still, so
422+
* the same colour means the same kcal/mol on every molecule in a series.
413423
*/
414-
const isoLock = new Map<FieldKind, { value: number, fittedOn: string }>();
415-
416424
function currentIso(): number {
417425
if (!activeKind) return 0;
418-
return (isoBase ?? Kinds[activeKind].iso) * isoMultiplier();
426+
return Kinds[activeKind].iso * isoMultiplier();
427+
}
428+
429+
/** The same contour expressed in the cube's own units. */
430+
function currentIsoCube(): number {
431+
if (!activeKind) return 0;
432+
return currentIso() * Kinds[activeKind].cubeScale;
419433
}
420434

421435
function updateIsoReadout() {
422436
const el = byId('field-iso-readout');
423437
if (!el) return;
424438
if (!activeKind) { el.textContent = ''; return; }
425439
const sign = Kinds[activeKind].diverging ? '±' : '';
426-
const lock = isoLock.get(activeKind);
427-
let provenance = '';
428-
if (lock && lock.fittedOn !== (ligandLabel ?? '')) {
429-
provenance = ` · locked (fitted on ${lock.fittedOn})`;
430-
} else if (lock) {
431-
provenance = ' · fitted to this molecule';
440+
let note = '';
441+
// Whether the surface actually closes is a fact about this picture, and it
442+
// is the one the crate bug turned on. Stated, not inferred from the shape.
443+
if (lastMeta && lastMeta.contour_closes_in_box === false) {
444+
note = ` · open surface (field is ${lastMeta.wall_max ?? '?'} at the box edge)`;
445+
} else if (lastMeta && lastMeta.single_signed) {
446+
note = ' · field is single-signed — no negative lobe exists';
432447
}
433448
el.textContent =
434-
`${sign}${currentIso().toPrecision(3)} ${Kinds[activeKind].unit}${provenance}`;
435-
const refit = byId<HTMLButtonElement>('field-refit');
436-
if (refit) refit.hidden = !lock || lock.fittedOn === (ligandLabel ?? '');
437-
}
438-
439-
/** Drop the held scale so the next render fits this molecule instead. */
440-
function refitIso() {
441-
if (activeKind) isoLock.delete(activeKind);
442-
const kind = activeKind;
443-
if (kind) void requestField(kind);
449+
`${sign}${currentIso().toPrecision(3)} ${Kinds[activeKind].unit}${note}`;
444450
}
445451

446452
function setButtonsEnabled() {
@@ -488,7 +494,7 @@ function reprParams(kind: FieldKind, sign: 1 | -1, shell: number) {
488494
return createVolumeRepresentationParams(plugin!, activeVolume ?? undefined, {
489495
type: 'isosurface',
490496
typeParams: {
491-
isoValue: Volume.IsoValue.absolute(sign * currentIso() * s.fraction),
497+
isoValue: Volume.IsoValue.absolute(sign * currentIsoCube() * s.fraction),
492498
visuals: s.visuals as ('solid' | 'wireframe')[],
493499
alpha: s.alpha,
494500
xrayShaded: s.visuals.includes('solid'),
@@ -518,17 +524,7 @@ async function renderCube(cubeText: string, kind: FieldKind, meta?: FieldMeta) {
518524
// Adopt the backend's measured isovalue when it sent one; fall back to the
519525
// kind's constant otherwise. Set BEFORE the representations are built, or
520526
// the first frame is drawn at the old scale and corrected a tick later.
521-
// Fit once, then hold. A held scale is what makes two molecules comparable;
522-
// re-fitting each time is the gain control this exists to avoid.
523-
const held = isoLock.get(kind);
524-
if (held) {
525-
isoBase = held.value;
526-
} else if (typeof meta?.iso_suggest === 'number') {
527-
isoBase = meta.iso_suggest;
528-
isoLock.set(kind, { value: meta.iso_suggest, fittedOn: ligandLabel ?? '' });
529-
} else {
530-
isoBase = null; // orbital/density: a normalised quantity keeps its constant
531-
}
527+
lastMeta = meta ?? null;
532528
if (!plugin) return;
533529
if (plugin.state.data.cells.has(REF_DATA)) {
534530
await plugin.build().delete(REF_DATA).commit();
@@ -736,7 +732,6 @@ export function initFieldWellsPanel(p: PluginContext) {
736732
// panel said "Solving…" for 36 minutes with every control disabled. The
737733
// daemon now stops itself too, but a bound the user cannot reach is not a
738734
// control, it is a promise.
739-
byId('field-refit')?.addEventListener('click', () => refitIso());
740735
byId('field-cancel')?.addEventListener('click', () => {
741736
abortInFlight();
742737
busy = false;

0 commit comments

Comments
 (0)