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,181 @@
+
+
+#include "inpdef.h"
+#include "ring_detection.h"
+#include "ichister.h"
+
+
+int find_atropisomeric_atoms_and_bonds(inp_ATOM* out_at,
+                                       int num_atoms,
+                                       RingSystems *ring_result,
+                                       ORIG_ATOM_DATA *orig_inp_data) {
+
+    int ret = 0;
+
+    if (out_at == NULL || num_atoms <= 0) {
+        return ret;
+    }
+
+    if (ring_result == NULL || orig_inp_data == NULL) {
+        return ret;
+    }
+
+    int fused_atom_partner[num_atoms];
+    for (int i = 0; i < num_atoms; i++) {
+        // out_at[i].fused_partner_atom_id = -1;
+        fused_atom_partner[i] = -1;
+    }
+    for (int 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;
+            }
+        }
+    }
+
+    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;
+                ret = 1;
+            } else if (score > 9) {
+                if (is_planar == 0) {
+                    is_atropisomer = 1;
+                    ret = 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->bAtropisomer = 1;
+
+                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]);
+        }
+    }
+    return ret;
+}

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

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

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,25 @@ 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->bAtropisomer = 0;
+    if (ip->Atropisomers) {
+        RingSystems *ring_result = find_rings(out_at, num_atoms);
+
+        // print_ring_result(ring_result);
+
+        int ret_ai = find_atropisomeric_atoms_and_bonds(out_at, num_atoms, ring_result, orig_inp_data);
+
+        //map values to orig_inp_data
+        if (orig_inp_data->bAtropisomer) {
+            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)

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

@@ -124,7 +124,8 @@ int set_common_options_by_parg(const char* pArg,
     int* pbMergeHash,
     int* pbHideInChI,
     int* pbMolecularInorganics,     /* @nnuk */
-    int* pbEnhancedStereochemistry);
+    int* pbEnhancedStereochemistry,
+    int* pbAtropisomers);
 
 
 /****************************************************************************
@@ -169,7 +170,8 @@ int set_common_options_by_parg(const char* pArg,
     int* pbMergeHash,
     int* pbHideInChI,
     int* pbMolecularInorganics,     /*@nnuk*/
-    int* pbEnhancedStereochemistry
+    int* pbEnhancedStereochemistry,
+    int* pbAtropisomers
 )
 {
     int got = 0;
@@ -275,6 +277,11 @@ int set_common_options_by_parg(const char* pArg,
         *pbEnhancedStereochemistry = 1;
         got = 1;
     }
+    else if (!inchi_stricmp(pArg, "Atropisomers"))
+    {
+        *pbAtropisomers = 1;
+        got = 1;
+    }
 
 #ifndef USE_STDINCHI_API
     /* These options DO TURN OFF Std flag   */
@@ -664,6 +671,7 @@ int ReadCommandLineParms(int argc,
     int bUnchargedAcidTaut = (CHARGED_SALTS_ONLY == 0);
     int bMergeSaltTGroups = (DISCONNECT_SALTS == 1);
     int bEnhancedStereochemistry = 0;
+    int bAtropisomers = 0;
 #if ( MIN_SB_RING_SIZE > 0 )
     int nMinDbRinSize = MIN_SB_RING_SIZE, mdbr = 0;
 #endif
@@ -816,7 +824,7 @@ int ReadCommandLineParms(int argc,
                 &bLargeMolecules, &bPolymers,
                 &bFoldPolymerSRU, &bFrameShiftScheme,
                 &bStereoAtZz, &bNPZz,
-                &bNoWarnings, &bMergeHash, &bHideInChI, &bMolecularInorganics, &bEnhancedStereochemistry);
+                &bNoWarnings, &bMergeHash, &bHideInChI, &bMolecularInorganics, &bEnhancedStereochemistry, &bAtropisomers);
             if (got)
             {
                 ;
@@ -1264,7 +1272,7 @@ int ReadCommandLineParms(int argc,
                 &bLargeMolecules, &bPolymers,
                 &bFoldPolymerSRU, &bFrameShiftScheme,
                 &bStereoAtZz, &bNPZz,
-                &bNoWarnings, &bMergeHash, &bHideInChI, &bMolecularInorganics, &bEnhancedStereochemistry);
+                &bNoWarnings, &bMergeHash, &bHideInChI, &bMolecularInorganics, &bEnhancedStereochemistry, &bAtropisomers);
 
             if ( got )
             {
@@ -2143,6 +2151,8 @@ int ReadCommandLineParms(int argc,
 
     ip->bEnhancedStereo = bEnhancedStereochemistry;
 
+    ip->Atropisomers = bAtropisomers;
+
     return 0;
 }
 

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

@@ -1619,6 +1619,14 @@ int OutputINChI1( CANON_GLOBALS *pCG,
                     }
                 }
 
+                if (ip->Atropisomers)
+                {
+                    if (set_Atropisomer_t_m_layers(orig_inp_data, pINChI, pINChI_Aux)) {
+                        io.sDifSegs[io.nCurINChISegment][DIFS_t_SATOMS] = DIFV_OUTPUT_FILL_T;
+                        io.sDifSegs[io.nCurINChISegment][DIFS_m_SP3INV] = DIFV_OUTPUT_FILL_T;
+                    }
+                }
+
                 if (ip->bEnhancedStereo)
                 {
                     set_EnhancedStereo_t_m_layers(orig_inp_data, pINChI, pINChI_Aux);
@@ -1714,6 +1722,10 @@ int OutputINChI1( CANON_GLOBALS *pCG,
         {
             is_beta = 1;
         }
+        else if (ip->Atropisomers)
+        {
+            is_beta = 1;
+        }
 
         OutputINCHI_VersionAndKind(out_file, strbuf, bINChIOutputOptions, is_beta, pLF, pTAB);
     }

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

@@ -106,7 +106,7 @@ static int CompDble( const void *a1, const void *a2, void * );
 static int Get2DTetrahedralAmbiguity( CANON_GLOBALS *pCG, double at_coord[][3], int bAddExplicitNeighbor, int bFix2DstereoBorderCase, double vMinAngle );
 static double triple_prod_and_min_abs_sine2( double at_coord[][3], double central_at_coord[], int bAddedExplicitNeighbor,
                                             double *min_sine, int *bAmbiguous, double vMinSine );
-static int are_4at_in_one_plane( double at_coord[][3], double min_sine );
+// static int are_4at_in_one_plane( double at_coord[][3], double min_sine );
 static int bInpAtomHasRequirdNeigh( inp_ATOM *at, int cur_at, int RequirdNeighType, int NumDbleBonds, int bStereoAtZz );
 static int bIsSuitableHeteroInpAtom( inp_ATOM  *at );
 static int bIsOxide( inp_ATOM  *at, int cur_at );
@@ -2190,7 +2190,7 @@ int half_stereo_bond_parity( inp_ATOM *at,
     {
         for (k = 0; k < 3; k++)
         {
-            at_coord[j][k] = 0; 
+            at_coord[j][k] = 0;
         }
     }
 

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

@@ -56,9 +56,11 @@ extern "C" {
     int bCanAtomBeMiddleAllene( char *elname, S_CHAR charge, S_CHAR radical );
     int bAtomHasValence3( char *elname, S_CHAR charge, S_CHAR radical );
 
-    double dot_prod3(const double a[], const double b[]); 
+    double dot_prod3(const double a[], const double b[]);
     void* cross_prod3(const double a[], const double b[], double result[]);
 
+    int are_4at_in_one_plane( double at_coord[][3], double min_sine );
+
     struct tagCANON_GLOBALS;
     int set_stereo_parity( struct tagCANON_GLOBALS *pCG,
                            inp_ATOM* at,

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

@@ -184,6 +184,8 @@ typedef struct tagInputAtom
     S_CHAR sb_parity[MAX_NUM_STEREO_BONDS];
     AT_NUMB sn_orig_at_num[MAX_NUM_STEREO_BONDS]; /* orig. at number of sn_ord[] neighbors        */
 
+    AT_NUMB bAtropisomeric; /* flag indicating if the atom is part of an atropisomeric bond */
+
 #if (FIND_RING_SYSTEMS == 1)
     S_CHAR bCutVertex;
     AT_NUMB nRingSystem;
@@ -459,7 +461,9 @@ typedef struct tagOrigAtom
     OAD_Polymer *polymer;
     OAD_V3000 *v3000;
     int valid_polymer;
-    int n_zy; /* number of non-polymeric pseudoatoms (Zy)             */
+    int n_zy;               /* number of non-polymeric pseudoatoms (Zy)             */
+
+    int bAtropisomer;     /* flag indicating whether the structure is an atropisomer; it is set to 1 if the structure has been identified as an atropisomer during input processing, and 0 otherwise */
 
 } ORIG_ATOM_DATA;