nx_chem_descriptors.nx source
↩ module page · 377 lines · 18098 B
1// nx_chem_descriptors.nx -- C2.7 milestone: drug-likeness descriptors.
2//
3// Counting primitives for Lipinski's Rule of Five + Veber-style drug-
4// likeness screening. All composable on the existing MolGraph + the
5// implicit-H valence inference (C2.3a).
6//
7// Lipinski Ro5 (1997, Adv. Drug Deliv. Rev. 23:3-25):
8// - MW <= 500 Da
9// - logP <= 5 (DEFERRED to C2.8 -- Crippen-Wildman atom-type table)
10// - HBD <= 5 (count of N-H, O-H, S-H bonds)
11// - HBA <= 10 (count of N, O atoms)
12// Pass if at most 1 violation (Lipinski's original wording).
13//
14// Veber 2002 (J. Med. Chem. 45:2615-2623):
15// - Rotatable bonds <= 10
16// - TPSA <= 140 Ų (DEFERRED to C2.8 -- Ertl-class atom-contribution table)
17//
18// Honest gaps (DEFERRED):
19// - logP (Crippen-Wildman 1999): C2.8 ~70-atom-type lookup table
20// - TPSA (Ertl-Rohde-Selzer 2000): C2.8 ~40-pattern lookup table
21// - SSSR (Smallest Set of Smallest Rings): C2.7.1; current approx
22// aromatic-ring count assumes 6-membered rings
23// - PAINS / structural-alert SMARTS filters: C2.9 (composes on
24// substructure search which is also pending)
25//
26// nx_safety_envelope:
27// intended_use: Lipinski/Veber drug-likeness screening for the
28// supplement-arc adulterant pre-filter; foundation
29// for nx_chem_label_claim evidence-checks
30// sil_target: SIL1
31// evidence: [C2.7 KAT in nx_chem_descriptors_test.nx with
32// reference molecules from Lipinski 1997 + Veber 2002]
33// verdict: BENCH-PENDING
34
35import "nx_chem.nx"
36import "nx_chem_molecule.nx"
37import "nx_chem_valence.nx"
38const K_MAGIC_500000: i64 = 500000
39
40// =================================================================
41// Count hydrogen bond donors: total H attached to N / O / S atoms.
42// Sums Atom.h_count for atoms with z in {7, 8, 16}. Caller MUST
43// run nx_chem_compute_implicit_h(m) FIRST so non-bracket atoms have
44// their h_count inferred.
45// =================================================================
46func nx_chem_count_hbd(m: *MolGraph) -> nx_int {
47 var total: nx_int = 0
48 var i: nx_int = 0
49 while i < m.n_atoms {
50 let a: *Atom = ((m.atoms as nx_int) + (i * NX_ATOM_BYTES)) as *Atom
51 if a.z == 7 {
52 if a.h_count > 0 { total = total + a.h_count }
53 }
54 if a.z == 8 {
55 if a.h_count > 0 { total = total + a.h_count }
56 }
57 if a.z == 16 {
58 if a.h_count > 0 { total = total + a.h_count }
59 }
60 i = i + 1
61 }
62 return total
63}
64
65// =================================================================
66// Count hydrogen bond acceptors: total atoms with z in {7, 8}.
67// Standard Lipinski definition (does not exclude HBD; HBD-overlap is
68// counted on both sides for Ro5 thresholds). Positively charged
69// ammonium (N+) is excluded since it can't accept H-bonds.
70// =================================================================
71func nx_chem_count_hba(m: *MolGraph) -> nx_int {
72 var total: nx_int = 0
73 var i: nx_int = 0
74 while i < m.n_atoms {
75 let a: *Atom = ((m.atoms as nx_int) + (i * NX_ATOM_BYTES)) as *Atom
76 if a.z == 7 {
77 if a.charge <= 0 { total = total + 1 }
78 }
79 if a.z == 8 { total = total + 1 }
80 i = i + 1
81 }
82 return total
83}
84
85// =================================================================
86// Ring perception: walk the bond graph as undirected DFS; for each
87// back-edge (a non-tree edge connecting two visited atoms), trace the
88// cycle path via parent pointers and mark all bonds + atoms on the
89// cycle as in_ring = 1.
90//
91// The SMILES parser only sets Bond.in_ring = 1 on the explicit ring-
92// closure bond; this pass marks the OTHER bonds + atoms in each ring.
93// Required before nx_chem_count_rotatable_bonds returns meaningful
94// results for ring-bearing molecules.
95//
96// Multi-ring + fused ring systems: each back-edge contributes its own
97// cycle; fused atoms appear in multiple cycles and get marked once
98// (in_ring is set to 1, no count).
99// =================================================================
100func nx_chem_descriptors_perceive_rings(m: *MolGraph) -> nx_int {
101 if m.n_atoms == 0 { return 0 }
102 let visited: *nx_int = (sys_mmap((m.n_atoms * 8) as i64)) as *nx_int
103 let parent: *nx_int = (sys_mmap((m.n_atoms * 8) as i64)) as *nx_int
104 let stack: *nx_int = (sys_mmap((m.n_atoms * 8) as i64)) as *nx_int
105 var pi: nx_int = 0
106 while pi < m.n_atoms {
107 parent[pi] = -1
108 pi = pi + 1
109 }
110 var start: nx_int = 0
111 while start < m.n_atoms {
112 if visited[start] == 0 {
113 var stack_n: nx_int = 0
114 stack[stack_n] = start
115 stack_n = stack_n + 1
116 visited[start] = 1
117 while stack_n > 0 {
118 stack_n = stack_n - 1
119 let u: nx_int = stack[stack_n]
120 var ei: nx_int = 0
121 while ei < m.n_bonds {
122 let b: *Bond = ((m.bonds as nx_int) + (ei * NX_BOND_BYTES)) as *Bond
123 var v: nx_int = -1
124 if b.a == u { v = b.b }
125 if b.b == u { v = b.a }
126 if v >= 0 {
127 if v != parent[u] {
128 if visited[v] == 1 {
129 // Back-edge: u-v forms a cycle through the DFS tree.
130 // Mark this back-edge bond + both endpoints.
131 b.in_ring = 1
132 let u_atom: *Atom = ((m.atoms as nx_int) + (u * NX_ATOM_BYTES)) as *Atom
133 u_atom.in_ring = 1
134 let v_atom: *Atom = ((m.atoms as nx_int) + (v * NX_ATOM_BYTES)) as *Atom
135 v_atom.in_ring = 1
136 // Only record each undirected back-edge once (when u > v).
137 if u > v {
138 // Find LCA(u, v) in the DFS tree.
139 // Step 1: walk u up to root, marking visited_u[].
140 let on_u_path: *nx_int = (sys_mmap((m.n_atoms * 8) as i64)) as *nx_int
141 var cu: nx_int = u
142 var done_u_walk: nx_int = 0
143 while done_u_walk == 0 {
144 if cu < 0 { done_u_walk = 1 }
145 else {
146 on_u_path[cu] = 1
147 cu = parent[cu]
148 }
149 }
150 // Step 2: walk v up; first atom on u's path is the LCA.
151 var cv: nx_int = v
152 var lca: nx_int = -1
153 var done_lca: nx_int = 0
154 while done_lca == 0 {
155 if cv < 0 { done_lca = 1 }
156 else {
157 if on_u_path[cv] == 1 { lca = cv; done_lca = 1 }
158 else { cv = parent[cv] }
159 }
160 }
161 // Step 3: mark bonds + atoms on path u -> LCA
162 var cur: nx_int = u
163 var done_u: nx_int = 0
164 while done_u == 0 {
165 if cur == lca { done_u = 1 }
166 else {
167 let next: nx_int = parent[cur]
168 if next < 0 { done_u = 1 }
169 else {
170 var ej: nx_int = 0
171 while ej < m.n_bonds {
172 let bj: *Bond = ((m.bonds as nx_int) + (ej * NX_BOND_BYTES)) as *Bond
173 var touches: nx_int = 0
174 if bj.a == cur {
175 if bj.b == next { touches = 1 }
176 }
177 if bj.b == cur {
178 if bj.a == next { touches = 1 }
179 }
180 if touches == 1 { bj.in_ring = 1 }
181 ej = ej + 1
182 }
183 let cur_atom: *Atom = ((m.atoms as nx_int) + (cur * NX_ATOM_BYTES)) as *Atom
184 cur_atom.in_ring = 1
185 let next_atom: *Atom = ((m.atoms as nx_int) + (next * NX_ATOM_BYTES)) as *Atom
186 next_atom.in_ring = 1
187 cur = next
188 }
189 }
190 }
191 // Step 4: mark bonds + atoms on path v -> LCA
192 var cur2: nx_int = v
193 var done_v: nx_int = 0
194 while done_v == 0 {
195 if cur2 == lca { done_v = 1 }
196 else {
197 let next2: nx_int = parent[cur2]
198 if next2 < 0 { done_v = 1 }
199 else {
200 var ek: nx_int = 0
201 while ek < m.n_bonds {
202 let bk: *Bond = ((m.bonds as nx_int) + (ek * NX_BOND_BYTES)) as *Bond
203 var touches2: nx_int = 0
204 if bk.a == cur2 {
205 if bk.b == next2 { touches2 = 1 }
206 }
207 if bk.b == cur2 {
208 if bk.a == next2 { touches2 = 1 }
209 }
210 if touches2 == 1 { bk.in_ring = 1 }
211 ek = ek + 1
212 }
213 let cur2_atom: *Atom = ((m.atoms as nx_int) + (cur2 * NX_ATOM_BYTES)) as *Atom
214 cur2_atom.in_ring = 1
215 let next2_atom: *Atom = ((m.atoms as nx_int) + (next2 * NX_ATOM_BYTES)) as *Atom
216 next2_atom.in_ring = 1
217 cur2 = next2
218 }
219 }
220 }
221 }
222 }
223 else {
224 // Tree edge
225 parent[v] = u
226 visited[v] = 1
227 stack[stack_n] = v
228 stack_n = stack_n + 1
229 }
230 }
231 }
232 ei = ei + 1
233 }
234 }
235 }
236 start = start + 1
237 }
238 return 0
239}
240
241// =================================================================
242// Heavy-atom degree of an atom (count of bonds to non-H neighbors).
243// Required for rotatable-bond detection (terminal heavy atoms have
244// degree 1; rotatable bonds must connect degree>=2 heavy atoms).
245// =================================================================
246func nx_chem_descriptors_heavy_degree(m: *MolGraph, atom_idx: nx_int) -> nx_int {
247 var deg: nx_int = 0
248 var bi: nx_int = 0
249 while bi < m.n_bonds {
250 let b: *Bond = ((m.bonds as nx_int) + (bi * NX_BOND_BYTES)) as *Bond
251 if b.a == atom_idx {
252 let other_idx: nx_int = b.b
253 let other_a: *Atom = ((m.atoms as nx_int) + (other_idx * NX_ATOM_BYTES)) as *Atom
254 if other_a.z != 1 { deg = deg + 1 }
255 }
256 if b.b == atom_idx {
257 let other_idx2: nx_int = b.a
258 let other_a2: *Atom = ((m.atoms as nx_int) + (other_idx2 * NX_ATOM_BYTES)) as *Atom
259 if other_a2.z != 1 { deg = deg + 1 }
260 }
261 bi = bi + 1
262 }
263 return deg
264}
265
266// =================================================================
267// Count rotatable bonds (Veber 2002 definition):
268// - bond order = single (NOT double/triple/aromatic)
269// - bond is NOT in a ring (in_ring == 0)
270// - both endpoints are heavy atoms (z != 1)
271// - both endpoints have heavy-degree >= 2 (not terminal)
272// - bond is NOT to a multiply-bonded sp/sp2 neighbor (skipped in MVP;
273// proper test requires hybridization assignment which is C2.8)
274//
275// MVP excludes terminal C-H rotors (trivially non-rotatable) and ring
276// bonds; doesn't yet exclude triple-bond-adjacent rotors per the
277// stricter Veber definition.
278// =================================================================
279func nx_chem_count_rotatable_bonds(m: *MolGraph) -> nx_int {
280 var total: nx_int = 0
281 var bi: nx_int = 0
282 while bi < m.n_bonds {
283 let b: *Bond = ((m.bonds as nx_int) + (bi * NX_BOND_BYTES)) as *Bond
284 if b.order == NX_BOND_SINGLE {
285 if b.in_ring == 0 {
286 let a_atom: *Atom = ((m.atoms as nx_int) + (b.a * NX_ATOM_BYTES)) as *Atom
287 let b_atom: *Atom = ((m.atoms as nx_int) + (b.b * NX_ATOM_BYTES)) as *Atom
288 if a_atom.z != 1 {
289 if b_atom.z != 1 {
290 let deg_a: nx_int = nx_chem_descriptors_heavy_degree(m, b.a)
291 let deg_b: nx_int = nx_chem_descriptors_heavy_degree(m, b.b)
292 if deg_a >= 2 {
293 if deg_b >= 2 { total = total + 1 }
294 }
295 }
296 }
297 }
298 }
299 bi = bi + 1
300 }
301 return total
302}
303
304// =================================================================
305// Count atoms with any aromaticity bit set. This is NOT a ring
306// count -- aromatic atoms can belong to fused-ring systems where
307// dividing by 6 over-counts. For naive "aromatic-ring approx",
308// see nx_chem_approx_aromatic_rings below.
309// =================================================================
310func nx_chem_count_aromatic_atoms(m: *MolGraph) -> nx_int {
311 var total: nx_int = 0
312 var i: nx_int = 0
313 while i < m.n_atoms {
314 let a: *Atom = ((m.atoms as nx_int) + (i * NX_ATOM_BYTES)) as *Atom
315 if a.aromaticity != 0 { total = total + 1 }
316 i = i + 1
317 }
318 return total
319}
320
321// =================================================================
322// Approximate aromatic ring count assuming benzene-class 6-membered
323// rings (most common in drug-like molecules). For naphthalene
324// (fused bicyclic 10 atoms) this approximation returns 10/6 = 1
325// (rounded down); proper SSSR is C2.7.1.
326// =================================================================
327func nx_chem_approx_aromatic_rings(m: *MolGraph) -> nx_int {
328 let arom_atoms: nx_int = nx_chem_count_aromatic_atoms(m)
329 return arom_atoms / 6
330}
331
332// =================================================================
333// Count heavy atoms (z >= 2; excludes H and wildcard).
334// Convenience wrapper composing on existing graph walk.
335// =================================================================
336func nx_chem_count_heavy_atoms(m: *MolGraph) -> nx_int {
337 var total: nx_int = 0
338 var i: nx_int = 0
339 while i < m.n_atoms {
340 let a: *Atom = ((m.atoms as nx_int) + (i * NX_ATOM_BYTES)) as *Atom
341 if a.z >= 2 { total = total + 1 }
342 i = i + 1
343 }
344 return total
345}
346
347// =================================================================
348// Partial Lipinski Rule of Five: MW + HBD + HBA components.
349// Returns 1 if molecule passes (at most 1 violation), 0 otherwise.
350// LogP component DEFERRED to C2.8 (Crippen-Wildman lookup table).
351//
352// Caller MUST run nx_chem_compute_implicit_h(m) first.
353// Takes the Element table for MW (typically nx_chem_periodic_table_118).
354// =================================================================
355func nx_chem_lipinski_ro5_partial(m: *MolGraph, table: *Element) -> nx_int {
356 let mw_q3: nx_int = nx_chem_mol_weight_q3(m, table)
357 let hbd: nx_int = nx_chem_count_hbd(m)
358 let hba: nx_int = nx_chem_count_hba(m)
359 var violations: nx_int = 0
360 if mw_q3 > K_MAGIC_500000 { violations = violations + 1 } // 500 Da in Q3 milli-AMU
361 if hbd > 5 { violations = violations + 1 }
362 if hba > 10 { violations = violations + 1 }
363 // Lipinski accepts at most 1 violation
364 if violations <= 1 { return 1 }
365 return 0
366}
367
368// =================================================================
369// Partial Veber rule: rotatable bonds <= 10.
370// TPSA component DEFERRED to C2.8 (Ertl atom-contribution table).
371// Returns 1 if passes, 0 otherwise.
372// =================================================================
373func nx_chem_veber_partial(m: *MolGraph) -> nx_int {
374 let rot: nx_int = nx_chem_count_rotatable_bonds(m)
375 if rot <= 10 { return 1 }
376 return 0
377}