Repository navigation
Expand file tree
/
Copy pathstereo.py
More file actions
1222 lines (1116 loc) · 45.7 KB
/
Copy pathstereo.py
File metadata and controls
1222 lines (1116 loc) · 45.7 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# -*- coding: utf-8 -*-
#
# Copyright 2019-2024 Ramil Nugmanov <nougmanoff@protonmail.com>
# This file is part of chython.
#
# chython is free software; you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, see <https://www.gnu.org/licenses/>.
#
from collections import defaultdict
from functools import cached_property
from itertools import combinations, product
from logging import getLogger, INFO
from typing import Dict, Set, Tuple, Union, List, Optional, TYPE_CHECKING
from .morgan import _morgan
from ..exceptions import AtomNotFound, IsChiral, NotChiral
logger = getLogger('chython.stereo')
logger.setLevel(INFO)
if TYPE_CHECKING:
from chython import MoleculeContainer
# atomic number constants
H = 1
C = 6
# 1 2
# \ |
# \|
# n---3
# /
# /
# 0
_tetrahedron_translate = {(0, 1, 2): False, (1, 2, 0): False, (2, 0, 1): False,
(0, 2, 1): True, (1, 0, 2): True, (2, 1, 0): True,
(0, 3, 1): False, (3, 1, 0): False, (1, 0, 3): False,
(0, 1, 3): True, (1, 3, 0): True, (3, 0, 1): True,
(0, 2, 3): False, (2, 3, 0): False, (3, 0, 2): False,
(0, 3, 2): True, (3, 2, 0): True, (2, 0, 3): True,
(1, 3, 2): False, (3, 2, 1): False, (2, 1, 3): False,
(1, 2, 3): True, (2, 3, 1): True, (3, 1, 2): True}
# 2 1
# \ /
# n---m
# / \
# 0 3
_alkene_translate = {(0, 1): False, (1, 0): False, (0, 3): True, (3, 0): True,
(2, 3): False, (3, 2): False, (2, 1): True, (1, 2): True}
def _pyramid_sign(n, u, v, w):
#
# | n /
# | |\
# | | \
# | /| \
# | / u---v
# |/___\_/___
# w
#
nx, ny, nz = n
ux, uy, uz = u
vx, vy, vz = v
wx, wy, wz = w
q1x = ux - nx
q1y = uy - ny
q1z = uz - nz
q2x = vx - nx
q2y = vy - ny
q2z = vz - nz
q3x = wx - nx
q3y = wy - ny
q3z = wz - nz
vol = q1x * (q2y * q3z - q2z * q3y) + q1y * (q2z * q3x - q2x * q3z) + q1z * (q2x * q3y - q2y * q3x)
if vol > 0:
return 1
elif vol < 0:
return -1
return 0
def _cis_trans_sign(n, u, v, w):
# n w
# \ /
# u--v
# / \
# x x
nx, ny = n
ux, uy = u
vx, vy = v
wx, wy = w
q1x = ux - nx
q1y = uy - ny
q2x = vx - ux
q2y = vy - uy
q3x = wx - vx
q3y = wy - vy
# cross vectors
q1q2z = q1x * q2y - q1y * q2x
q2q3z = q2x * q3y - q2y * q3x
dot = q1q2z * q2q3z
if dot > 0:
return 1
elif dot < 0:
return -1
return 0
def _allene_sign(mark, u, v, w):
# n w
# | /
# u--v
ux, uy = u
vx, vy = v
wx, wy = w
q2x = vx - ux
q2y = vy - uy
q3x = wx - vx
q3y = wy - vy
# cross vectors
q2q3z = q2x * q3y - q2y * q3x
dot = -mark * q2q3z
if dot > 0:
return 1
elif dot < 0:
return -1
return 0
class MoleculeStereo:
__slots__ = ()
def clean_stereo(self: 'MoleculeContainer'):
"""
Remove stereo data.
"""
for _, a in self.atoms():
a._stereo = None
for *_, b in self.bonds():
b._stereo = None
self.flush_cache(keep_sssr=True, keep_components=True)
@cached_property
def tetrahedrons(self: 'MoleculeContainer') -> Tuple[int, ...]:
"""
Carbon sp3 atom numbers.
"""
tetra = []
for n, atom in self.atoms():
if atom == C and not atom.charge and not atom.is_radical:
env = self._bonds[n]
if all(b == 1 for b in env.values()):
if sum(int(b) for b in env.values()) > 4:
continue
tetra.append(n)
return tuple(tetra)
@cached_property
def cumulenes(self: 'MoleculeContainer') -> List[Tuple[int, ...]]:
"""
All double-bonds chains (e.g. alkenes, allenes, cumulenes).
"""
atoms = self._atoms
bonds = self._bonds
adj = defaultdict(set) # double bonds adjacency matrix
for n, atom in atoms.items():
if atom.is_forming_double_bonds:
adj_n = adj[n].add
for m, bond in bonds[n].items():
if bond == 2 and atoms[m].is_forming_double_bonds:
adj_n(m)
if not adj:
return []
terminals = [x for x, y in adj.items() if len(y) == 1] # list to keep atoms order!
cumulenes = []
while terminals:
n = terminals.pop(0)
m = adj[n].pop()
path = [n, m]
while m not in terminals:
if len(bonds[m]) > 2: # not cumulene. SO3, SO4- etc.
cumulenes.extend(zip(path, path[1:])) # keep single double bonds instead of cumulene chain.
break
adj_m = adj[m]
adj_m.discard(n)
n, m = m, adj_m.pop()
path.append(m)
else:
terminals.remove(m)
adj[m].pop()
cumulenes.append(tuple(path))
return cumulenes
@cached_property
def stereogenic_tetrahedrons(self: 'MoleculeContainer') -> Dict[int, Union[Tuple[int, int, int], Tuple[int, int, int, int]]]:
"""
Tetrahedrons which contains at least 3 non-hydrogen neighbors and corresponding neighbors order.
"""
# 2
# |
# 1--K--3
# |
# 4?
atoms = self._atoms
bonds = self._bonds
tetrahedrons = {}
for n in self.tetrahedrons:
if any(not atoms[x].is_forming_single_bonds for x in bonds[n]):
continue # skip metal-carbon complexes
env = tuple(x for x in bonds[n] if atoms[x] != H)
if len(env) in (3, 4):
tetrahedrons[n] = env
return tetrahedrons
@cached_property
def stereogenic_cumulenes(self: 'MoleculeContainer') -> Dict[Tuple[int, ...], Tuple[int, int, Optional[int], Optional[int]]]:
"""
Cumulenes which contains at least one non-hydrogen neighbor on both ends and corresponding neighbors order.
"""
# 5 4
# \ /
# 2---3
# / \
# 1 6
bonds = self._bonds
atoms = self._atoms
cumulenes = {}
for path in self.cumulenes:
nf = bonds[path[0]]
nl = bonds[path[-1]]
n1, m1 = path[1], path[-2]
if any(b == 3 or not atoms[m].is_forming_single_bonds and b != 8
for m, b in nf.items() if m != n1):
continue # skip X=C=C structures and metal-carbon complexes
if any(b == 3 or not atoms[m].is_forming_single_bonds and b != 8
for m, b in nl.items() if m != m1):
continue # skip X=C=C structures and metal-carbon complexes
nn = [x for x, b in nf.items() if x != n1 and atoms[x] != H and b != 8]
mn = [x for x, b in nl.items() if x != m1 and atoms[x] != H and b != 8]
if nn and mn:
sn = nn[1] if len(nn) == 2 else None
sm = mn[1] if len(mn) == 2 else None
cumulenes[path] = (nn[0], mn[0], sn, sm)
return cumulenes
@cached_property
def stereogenic_allenes(self) -> Dict[int, Tuple[int, int, Optional[int], Optional[int]]]:
"""
Allenes which contains at least one non-hydrogen neighbor on both ends and corresponding neighbors order.
"""
return {path[len(path) // 2]: env for path, env in self.stereogenic_cumulenes.items() if len(path) % 2}
@cached_property
def stereogenic_cis_trans(self) -> Dict[Tuple[int, int], Tuple[int, int, Optional[int], Optional[int]]]:
"""
Cis-trans bonds which contains at least one non-hydrogen neighbor on both ends and corresponding neighbors order.
"""
stereo = {}
for path, env in self.stereogenic_cumulenes.items():
if len(path) % 2:
continue
stereo[(path[0], path[-1])] = env
return stereo
@cached_property
def ring_tetrahedrons(self: 'MoleculeContainer') -> Dict[int, Union[Tuple[int, int], Tuple[int], Tuple]]:
"""
Tetrahedrons in rings, except ring-linkers. Values are non-ring atoms.
"""
out = {}
atoms_rings = self.atoms_rings
tetrahedrons = self.stereogenic_tetrahedrons
points = self.rings_linker_tetrahedrons
environment = self.not_special_connectivity
for n, r in atoms_rings.items():
if n in tetrahedrons and n not in points:
out[n] = tuple(environment[n].difference(atoms_rings))
return out
@cached_property
def rings_linker_tetrahedrons(self: 'MoleculeContainer') -> Dict[int, Tuple[int, int, int, int]]:
"""
A dictionary where the keys are tetrahedron atoms shared between two rings (not condensed rings) and the values
are tuples representing their neighbors in the first and second rings respectively.
"""
out = {}
tetrahedrons = self.stereogenic_tetrahedrons
for n, r in self.atoms_rings.items():
if n in tetrahedrons:
for nr, mr in combinations(r, 2):
if len(set(nr).intersection(mr)) == 1:
ni = nr.index(n)
mi = mr.index(n)
out[n] = (nr[ni - 1], nr[ni - len(nr) + 1], mr[mi - 1], mr[mi - len(mr) + 1])
break
return out
@cached_property
def ring_cumulenes_terminals(self: 'MoleculeContainer') -> Set[Tuple[int, int]]:
"""
Terminal atoms of inside ring cumulenes.
"""
out = set()
ar = self.atoms_rings
for n, *_, m in self.stereogenic_cumulenes:
if n in ar and m in ar and not set(ar[n]).isdisjoint(ar[m]):
out.add((n, m))
return out
@cached_property
def rings_linker_cumulenes_terminals(self: 'MoleculeContainer') -> Dict[Tuple[int, int], Tuple[int, int, int, int]]:
"""
Terminal atoms of cumulenes connecting two rings. Values are neighbors in first and second rings.
"""
out = {}
ar = self.atoms_rings
chord = self.ring_cumulenes_terminals
for (n, *_, m), (n1, m1, n2, m2) in self.stereogenic_cumulenes.items():
if n in ar and m in ar and (n, m) not in chord:
out[(n, m)] = (n1, n2, m1, m2)
return out
@cached_property
def ring_attached_cumulenes(self: 'MoleculeContainer') -> Dict[Tuple[int, int], Union[Tuple[int, int], Tuple[int]]]:
"""
Cumulenes attached to rings from one side. Values are out of ring neighbor atoms.
"""
ar = self.atoms_rings
out = {}
for (n, *_, m), (n1, m1, n2, m2) in self.stereogenic_cumulenes.items():
if n in ar:
if m in ar:
continue
if m2:
out[(n, m)] = (m1, m2)
else:
out[(n, m)] = (m1,)
elif m in ar:
if n2:
out[(n, m)] = (n1, n2)
else:
out[(n, m)] = (n1,)
return out
@property
def chiral_tetrahedrons(self) -> Set[int]:
"""
Chiral tetrahedrons except already labeled ones.
"""
return self.__chiral_centers[0]
@property
def chiral_cis_trans(self) -> Set[Tuple[int, int]]:
"""
Chiral cis-trans bonds except already labeled ones.
"""
return self.__chiral_centers[1]
@property
def chiral_allenes(self) -> Set[int]:
"""
Chiral allenes except already labeled ones.
"""
return self.__chiral_centers[2]
def add_wedge(self: 'MoleculeContainer', n: int, m: int, mark: int, *, clean_cache=True):
"""
Add stereo data by wedge notation of bonds. Use it for tetrahedrons of allenes.
:param n: number of atom from which wedge bond started
:param m: number of atom to which wedge bond coming
:param mark: up bond is 1, down is -1
"""
atoms = self._atoms
if n not in atoms or m not in atoms or m not in self._bonds[n]:
raise AtomNotFound
elif atoms[n].stereo is not None:
raise IsChiral
elif c := self._stereo_allenes_centers.get(n):
# allenes
if atoms[c].stereo is not None:
raise IsChiral
elif c not in self.chiral_allenes:
raise NotChiral
t1, t2 = self._stereo_allenes_terminals[c]
order = self.stereogenic_allenes[c]
if atoms[m] == H:
if t1 == n:
m1 = order[1]
else:
t1, t2 = t2, t1
m1 = order[0]
r = True
else:
w = order.index(m)
if w == 0:
m1 = order[1]
r = False
elif w == 1:
m1 = order[0]
t1, t2 = t2, t1
r = False
elif w == 2:
m1 = order[1]
r = True
else:
m1 = order[0]
t1, t2 = t2, t1
r = True
if s := _allene_sign(mark, atoms[t1].xy, atoms[t2].xy, atoms[m1].xy):
atoms[c]._stereo = s < 0 if r else s > 0
if clean_cache:
self.flush_cache(keep_sssr=True, keep_components=True)
# tetrahedrons
elif n in self.chiral_tetrahedrons:
th = self.stereogenic_tetrahedrons[n]
am = atoms[m]
if am == H:
order = []
for x in th:
ax = atoms[x]
order.append((ax.x, ax.y, 0))
s = _pyramid_sign((am.x, am.y, mark), *order)
else:
order = []
for x in th:
ax = atoms[x]
order.append((ax.x, ax.y, mark if x == m else 0))
if len(order) == 3:
if len(self._bonds[n]) == 4: # explicit hydrogen
x = next(x for x in self._bonds[n] if x not in th)
ax = atoms[x]
s = _pyramid_sign((ax.x, ax.y, 0), *order)
else:
an = atoms[n]
s = _pyramid_sign((an.x, an.y, 0), *order)
else:
s = _pyramid_sign(order[-1], *order[:3])
if s:
atoms[n]._stereo = s > 0
if clean_cache:
self.flush_cache(keep_components=True, keep_sssr=True)
else:
raise NotChiral
def calculate_cis_trans_from_2d(self: 'MoleculeContainer', *, clean_cache=True):
"""
Calculate cis-trans stereo bonds from given 2d coordinates. Unusable for SMILES and INCHI.
"""
atoms = self._atoms
flag = False
while self.chiral_cis_trans:
stereo = False
for nm in self.chiral_cis_trans:
n, m = nm
n1, m1, *_ = self.stereogenic_cis_trans[nm]
s = _cis_trans_sign(atoms[n1].xy, atoms[n].xy, atoms[m].xy, atoms[m1].xy)
if s:
stereo = True
i, j = self._stereo_cis_trans_centers[n]
self._bonds[i][j]._stereo = s > 0
if stereo:
flag = True
self.flush_stereo_cache()
else:
break
if flag and clean_cache:
self.flush_cache(keep_components=True, keep_sssr=True)
def add_atom_stereo(self: 'MoleculeContainer', n: int, env: Tuple[int, ...], mark: bool, *, clean_cache=True):
"""
Add stereo data for specified neighbors bypass. Use it for tetrahedrons or allenes.
:param n: number of tetrahedron atom or central atom of allene.
:param env: numbers of atoms with specified bypass
:param mark: clockwise or anti bypass.
See <https://www.daylight.com/dayhtml/doc/theory/theory.smiles.html> and <http://opensmiles.org/opensmiles.html>
"""
try:
atom = self._atoms[n]
except KeyError:
raise AtomNotFound
if atom.stereo is not None:
raise IsChiral
if not isinstance(mark, bool):
raise TypeError('stereo mark should be bool')
if n in self.chiral_tetrahedrons:
atom._stereo = self._translate_tetrahedron_sign(n, env, mark)
if clean_cache:
self.flush_cache(keep_components=True, keep_sssr=True)
elif n in self.chiral_allenes:
atom._stereo = self._translate_allene_sign(n, *env, mark)
if clean_cache:
self.flush_cache(keep_components=True, keep_sssr=True)
else: # only tetrahedrons supported
raise NotChiral
def add_cis_trans_stereo(self: 'MoleculeContainer', n: int, m: int, n1: int, n2: int, mark: bool, *,
clean_cache=True):
"""
Add stereo data to cis-trans double bonds (not allenes).
n1/n=m/n2
:param n: number of starting atom of double bonds chain (alkenes of cumulenes)
:param m: number of ending atom of double bonds chain (alkenes of cumulenes)
:param n1: number of neighboring atom of starting atom
:param n2: number of neighboring atom of ending atom
:param mark: cis or trans
See <https://www.daylight.com/dayhtml/doc/theory/theory.smiles.html> and <http://opensmiles.org/opensmiles.html
"""
atoms = self._atoms
if n not in atoms or m not in atoms or n1 not in atoms or n2 not in atoms:
raise AtomNotFound
if not isinstance(mark, bool):
raise TypeError('stereo mark should be bool')
if n not in self._stereo_cis_trans_counterpart or self._stereo_cis_trans_counterpart[n] != m:
raise NotChiral
i, j = self._stereo_cis_trans_centers[n]
if self._bonds[i][j].stereo is not None:
raise IsChiral
if (n, m) in self.chiral_cis_trans:
self._bonds[i][j]._stereo = self._translate_cis_trans_sign(n, m, n1, n2, mark)
if clean_cache:
self.flush_cache(keep_components=True, keep_sssr=True)
elif (m, n) in self.chiral_cis_trans:
self._bonds[i][j]._stereo = self._translate_cis_trans_sign(m, n, n2, n1, mark)
if clean_cache:
self.flush_cache(keep_components=True, keep_sssr=True)
else:
raise NotChiral
def flush_stereo_cache(self):
"""
Flush chiral morgan and chiral centers cache.
"""
self.__dict__.pop('_chiral_morgan', None)
self.__dict__.pop('_MoleculeStereo__chiral_centers', None)
def fix_stereo(self: 'MoleculeContainer'):
"""
Reset stereo marks.
"""
stereo_tetrahedrons = self.stereogenic_tetrahedrons
stereo_allenes = self.stereogenic_allenes
stereo_cis_trans = self._stereo_cis_trans_terminals
atoms_stereo = []
allenes_stereo = []
cis_trans_stereo = []
for n, a in self.atoms():
if a.stereo is None:
continue
elif n in stereo_tetrahedrons:
atoms_stereo.append((n, a, a.stereo))
elif n in stereo_allenes:
allenes_stereo.append((n, a, a.stereo))
a._stereo = None # flush stereo label
for n, m, b in self.bonds():
if b.stereo is None:
continue
elif ta := stereo_cis_trans.get(n):
cis_trans_stereo.append((ta, b, b.stereo))
b._stereo = None # flush stereo label
self.flush_stereo_cache()
old_stereo = len(atoms_stereo) + len(allenes_stereo) + len(cis_trans_stereo)
while old_stereo:
chiral_tetrahedrons = self.chiral_tetrahedrons
chiral_allenes = self.chiral_allenes
chiral_cis_trans = self.chiral_cis_trans
# filter out resolved
tmp = []
for n, a, s in atoms_stereo:
if n in chiral_tetrahedrons:
a._stereo = s # restore stereo
else:
tmp.append((n, a, s))
atoms_stereo = tmp
tmp = []
for n, a, s in allenes_stereo:
if n in chiral_allenes:
a._stereo = s # restore stereo
else:
tmp.append((n, a, s))
allenes_stereo = tmp
tmp = []
for ta, b, s in cis_trans_stereo:
if ta in chiral_cis_trans:
b._stereo = s
else:
tmp.append((ta, b, s))
cis_trans_stereo = tmp
fail_stereo = len(atoms_stereo) + len(allenes_stereo) + len(cis_trans_stereo)
if fail_stereo == old_stereo:
break
old_stereo = fail_stereo
self.flush_stereo_cache()
@cached_property
def _cis_trans_count(self) -> int:
return sum(b.stereo is not None for *_, b in self.bonds())
@cached_property
def _stereo_cis_trans_centers(self) -> Dict[int, Tuple[int, int]]:
"""
Cis-Trans terminal atoms to cis-trans key mapping. Key is central double bond in a cumulene chain.
"""
terminals = {}
for path in self.stereogenic_cumulenes:
if len(path) % 2:
continue
n, m = path[0], path[-1]
i = len(path) // 2
terminals[n] = terminals[m] = (path[i - 1], path[i])
return terminals
@cached_property
def _stereo_cis_trans_terminals(self) -> Dict[int, Tuple[int, int]]:
"""
Cis-Trans terminal and central atoms to terminal pair mapping.
"""
terminals = {}
for path in self.stereogenic_cumulenes:
if len(path) % 2:
continue
n, m = path[0], path[-1]
i = len(path) // 2
terminals[n] = terminals[m] = terminals[path[i]] = terminals[path[i - 1]] = (n, m)
return terminals
@cached_property
def _stereo_cis_trans_counterpart(self) -> Dict[int, int]:
"""
Cis-Trans terminal atoms counterparts
"""
counterpart = {}
for path in self.stereogenic_cumulenes:
if len(path) % 2:
continue
n, m = path[0], path[-1]
counterpart[n] = m
counterpart[m] = n
return counterpart
@cached_property
def _stereo_allenes_centers(self) -> Dict[int, int]:
"""
Allene terminal atom to center mapping
"""
terminals = {}
for c, (n, m) in self._stereo_allenes_terminals.items():
terminals[n] = terminals[m] = c
return terminals
@cached_property
def _stereo_allenes_terminals(self) -> Dict[int, Tuple[int, int]]:
"""
Allene center atom to terminals mapping
"""
return {path[len(path) // 2]: (path[0], path[-1]) for path in self.stereogenic_cumulenes if len(path) % 2}
def _translate_tetrahedron_sign(self: 'MoleculeContainer', n, env, s=None):
"""
Get sign of chiral tetrahedron atom for specified neighbors order
:param n: stereo atom
:param env: neighbors order
:param s: if None, use existing sign else translate given to molecule
"""
if s is None:
s = self._atoms[n].stereo
if s is None:
raise KeyError
order = self.stereogenic_tetrahedrons[n]
if len(order) == 3:
if len(env) == 4: # hydrogen atom passed to env
# hydrogen always last in order
try:
order = (*order, next(x for x in env if self._atoms[x] == H)) # see translate scheme
except StopIteration:
raise KeyError
elif len(env) != 3: # pyramid or tetrahedron expected
raise ValueError('invalid atoms list')
elif len(env) not in (3, 4): # pyramid or tetrahedron expected
raise ValueError('invalid atoms list')
translate = tuple(order.index(x) for x in env[:3])
if _tetrahedron_translate[translate]:
return not s
return s
def _translate_cis_trans_sign(self: 'MoleculeContainer', n, m, nn, nm, s=None):
"""
Get sign for specified opposite neighbors
:param n: first double bonded atom
:param m: last double bonded atom
:param nn: neighbor of first atom
:param nm: neighbor of last atom
:param s: if None, use existing sign else translate given to molecule
"""
try:
n0, n1, n2, n3 = self.stereogenic_cis_trans[(n, m)]
except KeyError:
n0, n1, n2, n3 = self.stereogenic_cis_trans[(m, n)]
n, m = m, n # in alkenes sign not order depended
nn, nm = nm, nn
if s is None:
i, j = self._stereo_cis_trans_centers[n]
s = self._bonds[i][j].stereo
if s is None:
raise KeyError
if nn == n0: # same start
t0 = 0
if nm == n1:
t1 = 1
elif nm == n3 or n3 is None and self._atoms[nm] == H:
t1 = 3
else:
raise KeyError
elif nn == n1:
t0 = 1
if nm == n0:
t1 = 0
elif nm == n2 or n2 is None and self._atoms[nm] == H:
t1 = 2
else:
raise KeyError
elif nn == n2 or n2 is None and self._atoms[nn] == H:
t0 = 2
if nm == n1:
t1 = 1
elif nm == n3 or n3 is None and self._atoms[nm] == H:
t1 = 3
else:
raise KeyError
elif nn == n3 or n3 is None and self._atoms[nn] == H:
t0 = 3
if nm == n0:
t1 = 0
elif nm == n2 or n2 is None and self._atoms[nm] == H:
t1 = 2
else:
raise KeyError
else:
raise KeyError
if _alkene_translate[(t0, t1)]:
return not s
return s
def _translate_allene_sign(self: 'MoleculeContainer', c, nn, nm, s=None):
"""
get sign for specified opposite neighbors
:param c: central double bonded atom
:param nn: neighbor of first double bonded atom
:param nm: neighbor of last double bonded atom
:param s: if None, use existing sign else translate given to molecule
"""
if s is None:
s = self._atoms[c].stereo
if s is None:
raise KeyError
n0, n1, n2, n3 = self.stereogenic_allenes[c]
if nn == n0: # same start
t0 = 0
if nm == n1:
t1 = 1
elif nm == n3 or n3 is None and self._atoms[nm] == H:
t1 = 3
else:
raise KeyError
elif nn == n1:
t0 = 1
if nm == n0:
t1 = 0
elif nm == n2 or n2 is None and self._atoms[nm] == H:
t1 = 2
else:
raise KeyError
elif nn == n2 or n2 is None and self._atoms[nn] == H:
t0 = 2
if nm == n1:
t1 = 1
elif nm == n3 or n3 is None and self._atoms[nm] == H:
t1 = 3
else:
raise KeyError
elif nn == n3 or n3 is None and self._atoms[nn] == H:
t0 = 3
if nm == n0:
t1 = 0
elif nm == n2 or n2 is None and self._atoms[nm] == H:
t1 = 2
else:
raise KeyError
else:
raise KeyError
if _alkene_translate[(t0, t1)]:
return not s
return s
@cached_property
def _wedge_map(self: Union['MoleculeContainer', 'MoleculeStereo']):
atoms = self._atoms
overlap = set()
space = []
solved = []
seen = set()
for n, env in self.stereogenic_allenes.items():
if atoms[n].stereo is None:
continue
term = self._stereo_allenes_terminals[n]
overlap.update(term) # don't allow incoming wedge to allenes terminals
orders = [(*env[:2], *term, n, True), (*env[1::-1], *term[::-1], n, True)]
if env[2]:
orders.append((env[2], env[1], *term, n, True))
if env[3]:
orders.append((env[3], env[0], *term[::-1], n, True))
space.append(orders)
for n, env in self.stereogenic_tetrahedrons.items():
if atoms[n].stereo is None:
continue
overlap.add(n) # don't allow incoming wedge to stereo tetrahedrons
order = list(env)
orders = [(*order, n, False)]
for _ in range(1, len(order)):
order = order.copy()
order.append(order.pop(0))
orders.append((*order, n, False))
space.append(orders)
while space:
ls = len(space)
unsolved = []
for orders in space:
good = []
if orders[0][-1]:
for x in orders:
x0 = x[0]
if x0 in seen or x0 not in overlap:
good.append(x)
seen.add(x[2])
if good:
solved.append(max(good, key=lambda x: (not atoms[x[0]].in_ring, atoms[x[0]].atomic_number)))
else:
unsolved.append(orders)
else:
for x in orders:
x0 = x[0]
if x0 in seen or x0 not in overlap:
good.append(x)
if good:
seen.add(x[-2])
solved.append(max(good, key=lambda x: (not atoms[x[0]].in_ring, atoms[x[0]].atomic_number)))
else:
unsolved.append(orders)
space = unsolved
if len(unsolved) == ls:
break
solved = [y for x in solved if (y := self.__wedge_sign(x))]
if not space:
return solved
for orders in product(*space):
used = set()
wedge = []
for order in orders:
if order[-1]: # allene
if (order[0], order[2]) in used:
break
used.add((order[2], order[0]))
wedge.append(self.__wedge_sign(order))
else: # TH
n = order[-2]
if (order[0], n) in used:
break
used.add((n, order[0]))
wedge.append(self.__wedge_sign(order))
else: # found
solved.extend(wedge)
return solved
logger.info('wedge stereo mapping failed')
return solved
def __wedge_sign(self: 'MoleculeContainer', order):
if order[-1]: # allene
s = self._translate_allene_sign(order[-2], order[0], order[1])
v = _allene_sign(1, self._atoms[order[2]].xy, self._atoms[order[3]].xy, self._atoms[order[1]].xy)
if not v:
logger.info(f'need 2d clean. allenes wedge stereo ambiguous for atom {order[-2]}')
if s:
return order[2], order[0], v
else:
return order[2], order[0], -v
else: # TH
n = order[-2]
s = self._translate_tetrahedron_sign(n, order[:-2])
# need recalculation if XY changed
ao0 = self._atoms[order[0]]
ao1 = self._atoms[order[1]]
ao2 = self._atoms[order[2]]
if len(order) == 5:
an = self._atoms[n]
v = _pyramid_sign((an.x, an.y, 0),
(ao0.x, ao0.y, 1),
(ao1.x, ao1.y, 0),
(ao2.x, ao2.y, 0))
else:
ao3 = self._atoms[order[3]]
v = _pyramid_sign((ao3.x, ao3.y, 0),
(ao0.x, ao0.y, 1),
(ao1.x, ao1.y, 0),
(ao2.x, ao2.y, 0))
if not v:
logger.info(f'need 2d clean. tetrahedron wedge stereo ambiguous for atom {n}')
if s:
return n, order[0], v
else:
return n, order[0], -v
@cached_property
def _chiral_morgan(self: Union['MoleculeContainer', 'MoleculeStereo']) -> Dict[int, int]:
stereo_atoms = {n for n, a in self.atoms() if a.stereo is not None}
stereo_bonds = {n for n, mb in self._bonds.items() if any(b.stereo is not None for m, b in mb.items())}
if not stereo_atoms and not stereo_bonds:
return self.atoms_order
morgan = self.atoms_order.copy()
atoms_stereo = stereo_atoms.intersection(self.tetrahedrons)
allenes_stereo = stereo_atoms - atoms_stereo
cis_trans_terminals = self._stereo_cis_trans_terminals
cis_trans_stereo = {cis_trans_terminals[n] for n in stereo_bonds}
while True:
# try iteratively differentiate stereo atoms.
morgan, atoms_stereo, cis_trans_stereo, allenes_stereo, atoms_groups, cis_trans_groups, allenes_groups = \
self.__differentiation(morgan, atoms_stereo, cis_trans_stereo, allenes_stereo)
if not atoms_groups and not cis_trans_groups and not allenes_groups:
break
# for some rings differentiation by morgan impossible. try randomly set new weights.
# sometimes this will lead to pseudo chiral centers and non-unique morgan.
for group in atoms_groups:
for n in group[:len(group) // 2]: # set new weight in half of group randomly.
morgan[n] = -morgan[n]
for group in cis_trans_groups:
for n, _ in group[:len(group) // 2]: # set new weight in half of group randomly.
morgan[n] = -morgan[n]
for group in allenes_groups:
for n in group[:len(group) // 2]: # set new weight in half of group randomly.
morgan[n] = -morgan[n]
morgan = _morgan(morgan, self.int_adjacency)
return morgan
@cached_property
def __chiral_centers(self: Union['MoleculeStereo', 'MoleculeContainer']):
atoms_rings = self.atoms_rings