sasa.bend

bend-prot@0.1.0.0 · 84 lines · 10 definitions

1import Base
2import bend-prot@0.1.0.0/geom.bend as G
3import bend-prot@0.1.0.0/protein.bend as P
4import bend-prot@0.1.0.0/topology.bend as T
5
6# --- SASA lite v0: Shrake-Rupley with 6 octahedral sample points.
7# Uniform radii: rad scales the sample sphere, rc2 is the squared occlusion
8# cutoff, area_pt is the surface area each exposed point stands for (4*pi*r^2/6).
9# The occluder set `ys` must contain every atom INCLUDING the source atom;
10# the source atom is skipped by serial (which callers must keep unique),
11# because its own center sits exactly `rad` from each of its sample points
12# and would otherwise self-occlude whenever rc2 > rad^2.
13
14def sasa_offsets() -> List<&2, G.Vec3>:
15  [G.V3{1.0, 0.0, 0.0}, G.V3{(0.0 - 1.0 : F32), 0.0, 0.0}, G.V3{0.0, 1.0, 0.0}, G.V3{0.0, (0.0 - 1.0 : F32), 0.0}, G.V3{0.0, 0.0, 1.0}, G.V3{0.0, 0.0, (0.0 - 1.0 : F32)}]
16
17def occ_go(same: Bool, +pos: G.Vec3, +pt: G.Vec3, +rc2: F32) -> Bool:
18  match same:
19    case True{}:
20      False{}
21    case False{}:
22      F32.is_lt(G.Vec3.dist2(pos, pt), rc2)
23
24def occ_head(h: P.Atom, +pt: G.Vec3, +rc2: F32, skip: U32) -> Bool:
25  match h:
26    case P.Atom{serial, elem, +pos}:
27      occ_go(U32.is_eq(serial, skip), pos, pt, rc2)
28
29# Any atom but `skip` within rc2 of the sample point occludes it.
30def has_other_close(
31  xs: List<&2, P.Atom>, +pt: G.Vec3, +rc2: F32, +skip: U32
32) -> Bool:
33  match xs:
34    case Nil{}:
35      False{}
36    case h <> t:
37      Bool.or(has_other_close(t, pt, rc2, skip), occ_head(h, pt, rc2, skip))
38
39def sasa_pt(
40  +c: G.Vec3, +o: G.Vec3, +rad: F32, +ys: List<&2, P.Atom>, +rc2: F32,
41  +skip: U32
42) -> Nat:
43  T.hit_nat_go(Bool.not(has_other_close(ys, G.Vec3.add(c, G.Vec3.scale(rad, o)), rc2, skip)))
44
45def sasa_cover(
46  +c: G.Vec3, +offs: List<&2, G.Vec3>, +rad: F32, +ys: List<&2, P.Atom>,
47  +rc2: F32, +skip: U32
48) -> Nat:
49  match offs:
50    case Nil{}:
51      0n
52    case o <> t:
53      Nat.add(sasa_cover(c, t, rad, ys, rc2, skip), sasa_pt(c, o, rad, ys, rc2, skip))
54
55def sasa_of_atom(
56  h: P.Atom, +offs: List<&2, G.Vec3>, +rad: F32, +ys: List<&2, P.Atom>,
57  +rc2: F32, +area_pt: F32
58) -> F32:
59  match h:
60    case P.Atom{+serial, elem, +pos}:
61      (area_pt * F32.from_nat(sasa_cover(pos, offs, rad, ys, rc2, serial)) : F32)
62
63def sasa_all(
64  xs: List<&2, P.Atom>, +offs: List<&2, G.Vec3>, +rad: F32,
65  +ys: List<&2, P.Atom>, +rc2: F32, +area_pt: F32
66) -> List<&2, F32>:
67  match xs:
68    case Nil{}:
69      Nil{}
70    case h <> t:
71      sasa_of_atom(h, offs, rad, ys, rc2, area_pt) <> sasa_all(t, offs, rad, ys, rc2, area_pt)
72
73def sasa_of(
74  +xs: List<&2, P.Atom>, +rad: F32, +rc2: F32, +area_pt: F32
75) -> List<&2, F32>:
76  sasa_all(xs, sasa_offsets(), rad, xs, rc2, area_pt)
77
78# Total area: sum of per-atom areas.
79def sasa_total(areas: List<&2, F32>) -> F32:
80  match areas:
81    case Nil{}:
82      0.0
83    case h <> t:
84      (h + sasa_total(t) : F32)