
   0 | module Chem.Aromaticity
   1 | 
   2 | import Chem.Atom
   3 | import Chem.AtomType
   4 | import Chem.Elem
   5 | import Chem.Isotope
   6 | import Chem.Types
   7 | import Data.Graph.Indexed
   8 | import Data.Graph.Indexed.Util.Bipartite
   9 | import Data.Graph.Indexed.Ring.Relevant
  10 | import Data.SortedMap
  11 | import Data.SortedSet
  12 | import Derive.Prelude
  13 | 
  14 | %default total
  15 | %language ElabReflection
  16 | 
  17 | --------------------------------------------------------------------------------
  18 | -- Types
  19 | --------------------------------------------------------------------------------
  20 | 
  21 | public export
  22 | record AromBond (a : Type) where
  23 |   constructor AB
  24 |   type : a
  25 |   arom : Bool
  26 | 
  27 | %runElab derive "AromBond" [Show,Eq,Ord]
  28 | 
  29 | ||| True, if the given node is connected to at least one aromatic edge.
  30 | export %inline
  31 | isAromatic : IGraph k (AromBond e) n -> Fin k -> Bool
  32 | isAromatic g = any arom . neighboursAsAL g
  33 | 
  34 | export %inline
  35 | Cast a BondOrder => Cast (AromBond a) BondOrder where
  36 |   cast = cast . type
  37 | 
  38 | --------------------------------------------------------------------------------
  39 | -- AromBonds
  40 | --------------------------------------------------------------------------------
  41 | 
  42 | ||| Bonds count similar to `AtomType.Bonds` but including a counter for
  43 | ||| aromatic bonds.
  44 | public export
  45 | record AromBonds where
  46 |   constructor ABS
  47 |   single : Nat
  48 |   arom   : Nat
  49 |   double : Nat
  50 |   triple : Nat
  51 | 
  52 | ||| Total number of bonds
  53 | export
  54 | numBonds : AromBonds -> Nat
  55 | numBonds (ABS s a d t) = s + a + d + t
  56 | 
  57 | ||| Total bond order
  58 | |||
  59 | ||| Note: This assumes that there is a reasonable configuration of aromatic
  60 | |||       bonds (2 or 3). Otherwise, aromatic bonds are just ignored.
  61 | export
  62 | totalBondOrder : AromBonds -> Nat
  63 | totalBondOrder (ABS s a d t) =
  64 |   let rem = s + 2 * d + 3 * t
  65 |    in case a of
  66 |         2 => rem + 3
  67 |         3 => rem + 4
  68 |         _ => rem
  69 | 
  70 | --------------------------------------------------------------------------------
  71 | -- Utilities
  72 | --------------------------------------------------------------------------------
  73 | 
  74 | 0 Edges : Nat -> Type
  75 | Edges k = SortedSet (Fin k, Fin k)
  76 | 
  77 | toPair : Fin k -> Fin k -> (Fin k, Fin k)
  78 | toPair x y = if x > y then (y,x) else (x,y)
  79 | 
  80 | %inline hasEdge : Fin k -> Fin k -> Edges k -> Bool
  81 | hasEdge x y = contains (toPair x y)
  82 | 
  83 | -- we use `n & 7 == 4` here, because the contribution of
  84 | -- every node is counted twice
  85 | hueckel : Nat -> Bool
  86 | hueckel n = prim__and_Integer (cast n) 7 == 4
  87 | 
  88 | --------------------------------------------------------------------------------
  89 | -- Aromaticity Perception
  90 | --------------------------------------------------------------------------------
  91 | 
  92 | parameters {0 e,n   : Type}
  93 |            (contrib : n -> Maybe Nat)
  94 | 
  95 |   ||| Perceives aromaticity for a mol graph using the given contribution
  96 |   ||| function.
  97 |   export
  98 |   aromatizeI : {k : _} -> IGraph k e n -> IGraph k (AromBond e) n
  99 |   aromatizeI {k} g =
 100 |     let bs := foldl (pairs [] 0) empty (map ecycle . cr $ computeCrAndMCB g)
 101 |      in IG $ mapWithIndex (toArom bs) g.graph
 102 | 
 103 |     where
 104 |       toArom : Edges k -> Fin k -> Adj k e n -> Adj k (AromBond e) n
 105 |       toArom ps m = {neighbours $= mapKV (\n,v => AB v $ hasEdge m n ps)}
 106 | 
 107 |       pairs : List (Fin k, Fin k) -> Nat -> Edges k -> ECycle k -> Edges k
 108 |       pairs xs k es all@(E x y _ :: ys) =
 109 |         case [| contrib (lab g x) + contrib (lab g y) |] of
 110 |           Nothing => es
 111 |           Just v  => pairs (toPair x y :: xs) (k+v) es ys
 112 |       pairs xs k es [] = if hueckel k then foldl (flip insert) es xs else es
 113 | 
 114 |   ||| Perceives aromaticity for a mol graph using the given contribution
 115 |   ||| function.
 116 |   export %inline
 117 |   aromatize : Graph e n -> Graph (AromBond e) n
 118 |   aromatize (G o g) = G o $ aromatizeI g
 119 | 
 120 | --------------------------------------------------------------------------------
 121 | -- Models
 122 | --------------------------------------------------------------------------------
 123 | 
 124 | atMap : SortedMap (Elem,Charge,Hybridization) Nat
 125 | atMap =
 126 |   SortedMap.fromList
 127 |     [ ((C,0,SP2), 1)
 128 |     , ((C,0,SP), 1)
 129 |     , ((C,1,Planar), 0)
 130 |     , ((C,(-1),SP), 1)
 131 |     , ((C,(-1),SP2), 1)
 132 |     , ((C,(-1),SP3), 2)
 133 |     , ((B,0,SP2), 0)
 134 |     , ((O,0,SP3), 2)
 135 |     , ((S,0,SP3), 2)
 136 |     , ((N,0,SP3), 2)
 137 |     , ((N,0,SP2), 1)
 138 |     , ((N,(-1),SP3), 2)
 139 |     , ((N,(-1),SP2), 1)
 140 |     , ((N,1,SP2), 1)
 141 |     , ((P,0,SP3), 2)
 142 |     , ((P,0,SP2), 1)
 143 |     , ((P,(-1),SP3), 2)
 144 |     , ((P,(-1),SP2), 1)
 145 |     , ((P,1,SP2), 1)
 146 |     , ((As,0,SP3), 2)
 147 |     , ((Se,0,SP3), 2)
 148 |     ]
 149 | 
 150 | export
 151 | atomType : Cast e Elem => Atom e Charge p r h AtomType c l -> Maybe Nat
 152 | atomType a = lookup (cast a.elem, a.charge, a.type.hybridization) atMap
 153 | 
 154 | --------------------------------------------------------------------------------
 155 | -- Kekulization
 156 | --------------------------------------------------------------------------------
 157 | 
 158 | public export
 159 | 0 KAtom : (e,c,p,r,h,ch,l : Type) -> Type
 160 | KAtom e c p r h ch l = Atom e c p r h AtomType ch l
 161 | 
 162 | parameters {0 b,e,c,p,r,h,ch,l : Type}
 163 |            {auto cb : Cast b BondOrder}
 164 |            (setdbl  : b -> b)
 165 | 
 166 |   available : IGraph k b (KAtom e c p r h ch l) -> Fin k -> Bool
 167 |   available g x =
 168 |    let A l es := adj g x
 169 |     in S (count ((Dbl ==) . cast) es) == l.type.double
 170 | 
 171 |   export
 172 |   kekulize : Graph b (KAtom e c p r h ch l) -> Graph b (KAtom e c p r h ch l)
 173 |   kekulize (G k g) =
 174 |    let es := map setdbl <$> matchEdgesWhere g (available g)
 175 |     in G k $ insEdges es g
 176 |