Implementation of atropisomer recognition & encoding

INCHI-1-SRC/INCHI_API/libinchi/src/CMakeLists.txt

@@ -54,6 +54,8 @@ target_sources(libinchi PRIVATE
 	${P_IXA}/ixa_read_mol.c
 	${P_IXA}/ixa_status.c
 	${P_IXA}/ixa_status.h
+	${P_BASE}/atropisomers.h
+	${P_BASE}/atropisomers.c
 	${P_BASE}/bcf_s.h
 	${P_BASE}/bcf_s.c
 	${P_BASE}/extr_ct.h
@@ -127,6 +129,8 @@ target_sources(libinchi PRIVATE
 	${P_BASE}/permutation_util.c
 	${P_BASE}/readinch.c
 	${P_BASE}/readinch.h
+	${P_BASE}/ring_detection.c
+	${P_BASE}/ring_detection.h
 	${P_BASE}/runichi.c
 	${P_BASE}/runichi2.c
 	${P_BASE}/runichi3.c

INCHI-1-SRC/INCHI_BASE/src/atropisomers.c

@@ -0,0 +1,288 @@
+
+
+#include "inpdef.h"
+#include "ring_detection.h"
+#include "ichister.h"
+
+
+
+void find_atropisomeric_atoms_and_bonds_old(inp_ATOM* out_at,
+                                        int num_atoms,
+                                        RingSystems *ring_result,
+                                        ORIG_ATOM_DATA *orig_inp_data,
+                                        int *fused_atom_partner) {
+
+
+    for (int i = 0; i < num_atoms; i++) {
+
+        int atom_id1 = i;
+        const inp_ATOM atom_i = out_at[atom_id1];
+
+        int num_neighbors_i = atom_i.valence;
+
+        const AT_NUMB *neighbors = atom_i.neighbor;
+
+        if (num_neighbors_i == 3) {
+            for (int j = 0; j < num_neighbors_i; j++) {
+
+                int atom_id2 = neighbors[j];
+
+                if (atom_id1 >= atom_id2) {
+                    continue;
+                }
+
+                const inp_ATOM atom_j = out_at[atom_id2];
+
+                int num_neighbors_j = atom_j.valence;
+
+                if (num_neighbors_j == 3 &&
+                    atom_i.bond_stereo[j] == 0 &&
+                    atom_i.bond_type[j] == 1) {
+
+                    if (ring_result->atom_to_ring_mapping[atom_id1].ring_count > 0 ||
+                        ring_result->atom_to_ring_mapping[atom_id2].ring_count > 0) {
+
+                        int nof_wedge_bonds_i = 0;
+                        int has_double_bond_i = 0;
+                        for (int k = 0; k < num_neighbors_i; k++) {
+                            if (atom_i.bond_stereo[k] == 1 ||
+                                atom_i.bond_stereo[k] == 4 ||
+                                atom_i.bond_stereo[k] == 6) {
+                                nof_wedge_bonds_i++;
+                            }
+                            if (atom_i.bond_type[k] == 2) {
+                                has_double_bond_i = 1;
+                            }
+                        }
+                        int nof_wedge_bonds_j = 0;
+                        int has_double_bond_j = 0;
+                        for (int k = 0; k < num_neighbors_j; k++) {
+                            if (atom_j.bond_stereo[k] == 1 ||
+                                atom_j.bond_stereo[k] == 4 ||
+                                atom_j.bond_stereo[k] == 6) {
+                                nof_wedge_bonds_j++;
+                            }
+                            if (atom_j.bond_type[k] == 2) {
+                                has_double_bond_j = 1;
+                            }
+                        }
+
+                        if (nof_wedge_bonds_i > 0 ||
+                            nof_wedge_bonds_j > 0) {
+
+                            // if ((has_double_bond_i || has_double_bond_j) || 1==1) {
+                                if ((fused_atom_partner[atom_id1] != j &&
+                                     fused_atom_partner[atom_id2] != i) ||
+                                    (fused_atom_partner[atom_id1] == -1 ||
+                                     fused_atom_partner[atom_id2] == -1)) {
+
+                                    if (are_atoms_in_same_small_ring(out_at,
+                                                                     num_atoms,
+                                                                     ring_result,
+                                                                     atom_id1, atom_id2,
+                                                                     6) == 0) {
+                                        printf(">>> is atropisomer\n");
+                                        printf("infos: atom %d with atom %d; bond type %d; bond stereo %d\n",
+                                                atom_id1, atom_id2, atom_i.bond_type[j], atom_i.bond_stereo[j]);
+                                        // printf("atom type %d %d\n", atom_i.el_number, atom_j.el_number);
+                                        // printf("has double bond %d %d\n", has_double_bond_i, has_double_bond_j);
+                                        // printf("fused pivot atoms %d %d\n", fused_atom_partner[atom_id1], fused_atom_partner[atom_id2]);
+
+                                        //TODO set atoms with atropisomeric bonds in out_at
+                                        orig_inp_data->is_atropisomer = 1;
+                                    }
+                                }
+                            // }
+                        } else {
+                            // printf(">>> not atropisomer because no wedge bonds: atom %d with atom %d; bond type %d; bond stereo %d\n",
+                            //         atom_id1, atom_id2, atom_i.bond_type[j], atom_i.bond_stereo[j]);
+                        }
+
+                    }
+                }
+            }
+        }
+    }
+
+    for (int i = 0; i < ring_result->count; i++) {
+        const Ring *cur_ring = &ring_result->rings[i];
+        if (cur_ring->size >= 8) {
+            int count_single_bonds = 0;
+            int count_double_bonds = 0;
+            for (int j = 0; j < cur_ring->size; j++) {
+                int atom_id1 = cur_ring->atom_ids[j];
+                int atom_id2 = cur_ring->atom_ids[(j + 1) % cur_ring->size];
+                const inp_ATOM atom1 = out_at[atom_id1];
+                for (int k = 0; k < atom1.valence; k++) {
+                    if (atom1.neighbor[k] == atom_id2) {
+                        if (atom1.bond_type[k] == 1) {
+                            count_single_bonds++;
+                        } else if (atom1.bond_type[k] == 2) {
+                            count_double_bonds++;
+                        }
+                    }
+                }
+            }
+
+            if (cur_ring->is_fused_ring == 0) {
+                printf(">>> ring id %d size %d nof single bonds %d nof double bonds %d fused ring %d\n",
+                    cur_ring->id, cur_ring->size,
+                    count_single_bonds, count_double_bonds,
+                    cur_ring->is_fused_ring);
+            }
+        }
+    }
+}
+
+void find_atropisomeric_atoms_and_bonds(inp_ATOM* out_at,
+                                        int num_atoms,
+                                        RingSystems *ring_result,
+                                        ORIG_ATOM_DATA *orig_inp_data,
+                                        int *fused_atom_partner) {
+
+
+    for (int i = 0; i < num_atoms; i++) {
+
+        int atom_id1 = i;
+        const inp_ATOM atom_i = out_at[atom_id1];
+
+        int num_neighbors_i = atom_i.valence;
+
+        const AT_NUMB *neighbors = atom_i.neighbor;
+
+        double at_coord[4][3];
+        at_coord[0][0] = atom_i.x;
+        at_coord[0][1] = atom_i.y;
+        at_coord[0][2] = atom_i.z;
+
+        if (num_neighbors_i < 3) {
+            continue;
+        }
+
+        for (int j = 0; j < num_neighbors_i; j++) {
+            int score = 0;
+
+            int atom_id2 = neighbors[j];
+
+            if (atom_id1 >= atom_id2) {
+                continue;
+            }
+
+            const inp_ATOM atom_j = out_at[atom_id2];
+
+            int num_neighbors_j = atom_j.valence;
+
+            if (num_neighbors_j < 3) {
+                continue;
+            }
+
+            score += (num_neighbors_j == 3);
+
+            score += (atom_i.bond_stereo[j] == 0);
+            score -= ((atom_i.bond_stereo[j] == STEREO_SNGL_UP) ||
+                      (atom_i.bond_stereo[j] == STEREO_SNGL_EITHER) ||
+                      (atom_i.bond_stereo[j] == STEREO_SNGL_DOWN)) * 2;
+
+            score += (atom_i.bond_type[j] == 1);
+            score -= (atom_i.bond_type[j] != 1) * 2;
+
+            score += (ring_result->atom_to_ring_mapping[atom_id1].ring_count > 0) * 3;
+            score += (ring_result->atom_to_ring_mapping[atom_id2].ring_count > 0) * 3;
+
+
+            int nof_wedge_bonds_i = 0;
+            int has_double_bond_i = 0;
+            int coord_count = 1;
+            for (int k = 0; k < num_neighbors_i; k++) {
+                if (atom_i.bond_stereo[k] == 1 ||
+                    atom_i.bond_stereo[k] == 4 ||
+                    atom_i.bond_stereo[k] == 6) {
+                    nof_wedge_bonds_i++;
+                }
+                if (atom_i.bond_type[k] == 2) {
+                    has_double_bond_i = 1;
+                }
+                if (atom_id2 != atom_i.neighbor[k]) {
+                    if (coord_count <= 2) {
+                        at_coord[coord_count][0] = out_at[atom_i.neighbor[k]].x;
+                        at_coord[coord_count][1] = out_at[atom_i.neighbor[k]].y;
+                        at_coord[coord_count][2] = out_at[atom_i.neighbor[k]].z;
+                        coord_count++;
+                    }
+                }
+            }
+            int nof_wedge_bonds_j = 0;
+            int has_double_bond_j = 0;
+            for (int k = 0; k < num_neighbors_j; k++) {
+                if (atom_j.bond_stereo[k] == 1 ||
+                    atom_j.bond_stereo[k] == 4 ||
+                    atom_j.bond_stereo[k] == 6) {
+                    nof_wedge_bonds_j++;
+                }
+                if (atom_j.bond_type[k] == 2) {
+                    has_double_bond_j = 1;
+                }
+                if (atom_id1 != atom_j.neighbor[k]) {
+                    at_coord[3][0] = out_at[atom_j.neighbor[k]].x;
+                    at_coord[3][1] = out_at[atom_j.neighbor[k]].y;
+                    at_coord[3][2] = out_at[atom_j.neighbor[k]].z;
+                }
+            }
+
+            score += (nof_wedge_bonds_i > 0) * 4;
+            score += (nof_wedge_bonds_j > 0) * 4;
+
+            score -= (nof_wedge_bonds_i == 0) * 3;
+            score -= (nof_wedge_bonds_j == 0) * 3;
+
+            score += has_double_bond_i * 2;
+            score += has_double_bond_j * 2;
+
+            score += (fused_atom_partner[atom_id1] != j && fused_atom_partner[atom_id2] != i);
+            score += (fused_atom_partner[atom_id1] == -1 || fused_atom_partner[atom_id2] == -1);
+
+            int both_atoms_in_same_small_ring = are_atoms_in_same_small_ring(out_at,
+                                                                        num_atoms,
+                                                                        ring_result,
+                                                                        atom_id1, atom_id2,
+                                                                        6);
+
+            score += (both_atoms_in_same_small_ring == 0);
+
+            // printf("are atom in small ring? %d %d\n", atom_in_same_small_ring, (atom_in_same_small_ring == 1) * 10);
+            score -= (both_atoms_in_same_small_ring == 1) * 20;
+
+            int is_planar = are_4at_in_one_plane(at_coord, 0.03);
+
+            int is_atropisomer = 0;
+
+            if (score > 10) {
+                is_atropisomer = 1;
+            } else if (score > 9) {
+                if (is_planar == 0) {
+                    is_atropisomer = 1;
+                } else {
+                    // TODO are more tests/rules needed??
+                }
+            }
+
+            if (is_atropisomer == 1) {
+                printf(">>> FOUND atropisomer (higher score): atom id %2d atom id %2d  is planar %d  --> score %2d (%d)\n", atom_id1, atom_id2, is_planar, score, both_atoms_in_same_small_ring);
+                orig_inp_data->is_atropisomer = 1;
+
+                //TODO this rule is not enough for setting the flag
+                orig_inp_data->is_enantiomeric_atropisomer = ((nof_wedge_bonds_i > 0) && (nof_wedge_bonds_j > 0));
+
+                out_at[atom_id1].bAtropisomeric = 1;
+                out_at[atom_id2].bAtropisomeric = 1;
+            }
+
+            // printf(">>> is atropisomer\n");
+            // printf("infos: atom %d with atom %d; bond type %d; bond stereo %d\n",
+            //         atom_id1, atom_id2, atom_i.bond_type[j], atom_i.bond_stereo[j]);
+            // printf("atom type %d %d\n", atom_i.el_number, atom_j.el_number);
+            // printf("has double bond %d %d\n", has_double_bond_i, has_double_bond_j);
+            // printf("fused pivot atoms %d %d\n", fused_atom_partner[atom_id1], fused_atom_partner[atom_id2]);
+        }
+    }
+}

INCHI-1-SRC/INCHI_BASE/src/atropisomers.h

@@ -0,0 +1,9 @@
+
+#include "inpdef.h"
+#include "ring_detection.h"
+
+void find_atropisomeric_atoms_and_bonds(inp_ATOM* out_at,
+                                        int num_atoms,
+                                        RingSystems *ring_result,
+                                        ORIG_ATOM_DATA *orig_inp_data,
+                                        int *fused_atom_partner);

INCHI-1-SRC/INCHI_BASE/src/ichidrp.h

@@ -189,7 +189,8 @@ typedef struct tagInputParms {
     int             bMolecularInorganics;       /* (@nnuk : Nauman Ullah Khan) :: Varaible for Molecular Inorganics parameter */
     int             bMolecularInorganicsReconnectedInChI;  /* (@nnuk : Nauman Ullah Khan) :: Custom flag to indicate reconnected InChI requirement */
 
-    int             bEnhancedStereo;
+    int             bEnhancedStereo;        /* v. 1.0?+ enable enhanced stereochemistry                             */
+    int             Atropisomers;           /* v. 1.0?+ enable atropisomeric stereochemistry                        */
 
     /* */
     INCHI_MODE      bTautFlags;

INCHI-1-SRC/INCHI_BASE/src/ichimake.c

@@ -54,6 +54,9 @@
 
 #include "bcf_s.h"
 
+#include "atropisomers.h"
+#include "ring_detection.h"
+
 /*
     Local functions
 */
@@ -3901,6 +3904,39 @@ int  Create_INChI(CANON_GLOBALS* pCG,
     /*fix_odd_things( num_atoms, out_at );*/
 #if ( FIND_RING_SYSTEMS == 1 )
     MarkRingSystemsInp(out_at, num_atoms, 0);
+
+    orig_inp_data->is_atropisomer = 0;
+    orig_inp_data->is_enantiomeric_atropisomer = 0;
+    if (ip->Atropisomers) {
+        RingSystems *ring_result = find_rings(out_at, num_atoms);
+
+        // print_ring_result(ring_result);
+        int fused_atom_partner[num_atoms];
+        for (i = 0; i < num_atoms; i++) {
+            // out_at[i].fused_partner_atom_id = -1;
+            fused_atom_partner[i] = -1;
+        }
+        for (i = 0; i < num_atoms; i++) {
+            for (int j = i + 1; j < num_atoms; j++) {
+                if(is_fused_ring_pivot(ring_result, out_at, i, j)) {
+                    fused_atom_partner[i] = j;
+                    fused_atom_partner[j] = i;
+                }
+            }
+        }
+
+        find_atropisomeric_atoms_and_bonds(out_at, num_atoms, ring_result, orig_inp_data, fused_atom_partner);
+
+        //map values to orig_inp_data
+        if (orig_inp_data->is_atropisomer) {
+            for (i = 0; i < num_atoms; i++) {
+                orig_inp_data->at[i].bAtropisomeric = out_at[i].bAtropisomeric;
+            }
+        }
+
+        free_ring_system(ring_result);
+    }
+
 #endif
     /*  duplicate the preprocessed structure so that all supplied out_norm_data[]->at buffers are filled */
     if (out_at != out_norm_data[TAUT_YES]->at && out_norm_data[TAUT_YES]->at)