1	import Base
2	import ./geom.bend as G
3	import ./protein.bend as P
4	import ./topology.bend as T
5	
6	# --- RMSD v0: centroid-free aligned RMSD (no rotation) + ensemble matrix.
7	# Lockstep pairs zip: a longer list's extra elements are ignored, so
8	# `rmsd` of unequal lists silently answers 0. Use `rmsd_checked` when the
9	# mismatch must be visible instead.
10	
11	def pair_head_go(+p1: G.Vec3, y: P.Atom, rest: F32 & Nat) -> F32 & Nat:
12	  match y:
13	    case P.Atom{s2, e2, +p2}:
14	      T.sumd2_extend(G.Vec3.dist2(p1, p2), rest)
15	
16	def pair_head(x: P.Atom, y: P.Atom, rest: F32 & Nat) -> F32 & Nat:
17	  match x:
18	    case P.Atom{s1, e1, +p1}:
19	      pair_head_go(p1, y, rest)
20	
21	def pair_sumd2(xs: List<&2, P.Atom>, ys: List<&2, P.Atom>) -> F32 & Nat:
22	  match xs ys:
23	    case Nil{} Nil{}:
24	      (0.0, 0n)
25	    case Nil{} qh <> qt:
26	      (0.0, 0n)
27	    case h <> t Nil{}:
28	      (0.0, 0n)
29	    case x <> xt y <> yt:
30	      pair_head(x, y, pair_sumd2(xt, yt))
31	
32	def rmsd(xs: List<&2, P.Atom>, ys: List<&2, P.Atom>) -> F32:
33	  T.rg_of(pair_sumd2(xs, ys))
34	
35	# Length-checked RMSD: None{} on unequal lengths, Some{...} otherwise.
36	def rmsd_checked_go(
37	  ok: Bool, +xs: List<&2, P.Atom>, +ys: List<&2, P.Atom>
38	) -> Maybe<&2, F32>:
39	  match ok:
40	    case True{}:
41	      Some{rmsd(xs, ys)}
42	    case False{}:
43	      None{}
44	
45	def rmsd_checked(
46	  +xs: List<&2, P.Atom>, +ys: List<&2, P.Atom>
47	) -> Maybe<&2, F32>:
48	  rmsd_checked_go(Nat.is_eq(List.length(&2, P.Atom, xs), List.length(&2, P.Atom, ys)), xs, ys)
49	
50	def rmsd_against(+q: List<&2, P.Atom>, h: List<&2, P.Atom>) -> F32:
51	  T.rg_of(pair_sumd2(q, h))
52	
53	def rmsd_to_all(
54	  +q: List<&2, P.Atom>, ys: List<&2, List<&2, P.Atom>>
55	) -> List<&2, F32>:
56	  match ys:
57	    case Nil{}:
58	      Nil{}
59	    case h <> t:
60	      rmsd_against(q, h) <> rmsd_to_all(q, t)
61	
62	def rmsd_matrix(
63	  +xs: List<&2, List<&2, P.Atom>>, +all: List<&2, List<&2, P.Atom>>
64	) -> List<&2, List<&2, F32>>:
65	  match xs:
66	    case Nil{}:
67	      Nil{}
68	    case h <> t:
69	      rmsd_to_all(h, all) <> rmsd_matrix(t, all)
70	
71	def rmsd_matrix_of(+ens: List<&2, List<&2, P.Atom>>) -> List<&2, List<&2, F32>>:
72	  rmsd_matrix(ens, ens)