1 /*
2 * International Chemical Identifier (InChI)
3 * Version 1
4 * Software version 1.04
5 * September 9, 2011
6 *
7 * The InChI library and programs are free software developed under the
8 * auspices of the International Union of Pure and Applied Chemistry (IUPAC).
9 * Originally developed at NIST. Modifications and additions by IUPAC
10 * and the InChI Trust.
11 *
12 * IUPAC/InChI-Trust Licence for the International Chemical Identifier (InChI)
13 * Software version 1.0.
14 * Copyright (C) IUPAC and InChI Trust Limited
15 *
16 * This library is free software; you can redistribute it and/or modify it under the
17 * terms of the IUPAC/InChI Trust Licence for the International Chemical Identifier
18 * (InChI) Software version 1.0; either version 1.0 of the License, or
19 * (at your option) any later version.
20 *
21 * This library is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
24 * See the IUPAC/InChI Trust Licence for the International Chemical Identifier (InChI)
25 * Software version 1.0 for more details.
26 *
27 * You should have received a copy of the IUPAC/InChI Trust Licence for the
28 * International Chemical Identifier (InChI) Software version 1.0 along with
29 * this library; if not, write to:
30 *
31 * The InChI Trust
32 * c/o FIZ CHEMIE Berlin
33 * Franklinstrasse 11
34 * 10587 Berlin
35 * GERMANY
36 *
37 */
38
39
40 #include <stdio.h>
41 #include <stdlib.h>
42 #include <math.h>
43 #include <string.h>
44
45 #include "mode.h"
46
47 #include "ichierr.h"
48 #include "inpdef.h"
49 #include "extr_ct.h"
50 #include "ichister.h"
51 #include "ichiring.h"
52 #include "ichi.h"
53
54 #include "ichicomp.h"
55 #include "util.h"
56
57 #define ZTYPE_DOWN (-1) /* should be equal to -ZTYPE_UP */
58 #define ZTYPE_NONE 0
59 #define ZTYPE_UP 1 /* should be equal to -ZTYPE_DOWN */
60 #define ZTYPE_3D 3
61 #define ZTYPE_EITHER 9999
62
63 /* criteria for ill-defined */
64 #define MIN_ANGLE 0.10 /* 5.73 degrees */
65 #define MIN_SINE 0.03 /* min edge/plane angle in case the tetrahedra has significantly different edge length */
66 #define MIN_ANGLE_DBOND 0.087156 /* 5 degrees = max angle considered as too small for unambiguous double bond stereo */
67 #define MIN_SINE_OUTSIDE 0.06 /* min edge/plane angle to determine whether the central atom is outside of the tetrahedra */
68 #define MIN_SINE_SQUARE 0.125 /* min edge/plane angle in case the tetrahedra is somewhat close to a parallelogram */
69 #define MIN_SINE_EDGE 0.167 /* min sine/(min.edge) ratio to avoid undefined in case of long edges */
70 #define MIN_LEN_STRAIGHT 1.900 /* min length of two normalized to 1 bonds in a straight line */
71 #define MAX_SINE 0.70710678118654752440084436210485 /* 1/sqrt(2)=sin(pi/4) */
72 #define MIN_BOND_LEN 0.000001
73 #define ZERO_LENGTH MIN_BOND_LEN
74 #define ZERO_FLOAT 1.0e-12
75 #define BOND_PARITY_UNDEFINED 64
76 #if ( STEREO_CENTER_BONDS_NORM == 1 )
77 #define MPY_SINE 1.00 /* was 3.0 */
78 #define MAX_EDGE_RATIO 2.50 /* max max/min edge ratio for a tetrahedra close to a parallelogram */
79 #else
80 #define MPY_SINE 3.00
81 #define MAX_EDGE_RATIO 6.00 /* max max/min edge ratio for a tetrahedra close to a parallelogram */
82 #endif
83 /* local prototypes */
84 static int save_a_stereo_bond( int z_prod, int result_action,
85 int at1, int ord1, AT_NUMB *stereo_bond_neighbor1, S_CHAR *stereo_bond_ord1, S_CHAR *stereo_bond_z_prod1, S_CHAR *stereo_bond_parity1,
86 int at2, int ord2, AT_NUMB *stereo_bond_neighbor2, S_CHAR *stereo_bond_ord2, S_CHAR *stereo_bond_z_prod2, S_CHAR *stereo_bond_parity2 );
87 static double get_z_coord( inp_ATOM* at, int cur_atom, int neigh_no, int *nType,int bPointedEdgeStereo );
88 static double len3( const double c[] );
89 static double len2( const double c[] );
90 static double* diff3( const double a[], const double b[], double result[] );
91 static double* add3( const double a[], const double b[], double result[] );
92 static double* mult3( const double a[], double b, double result[] );
93 static double* copy3( const double a[], double result[] );
94 static double* change_sign3( const double a[], double result[] );
95 static double dot_prod3( const double a[], const double b[] );
96 static int dot_prodchar3( const S_CHAR a[], const S_CHAR b[] );
97 static double* cross_prod3( const double a[], const double b[], double result[] );
98 static double triple_prod( double a[], double b[], double c[], double *sine_value );
99 static double triple_prod_and_min_abs_sine(double at_coord[][3], double *min_sine);
100 static int are_3_vect_in_one_plane( double at_coord[][3], double min_sine);
101 static int triple_prod_char( inp_ATOM *at, int at_1, int i_next_at_1, S_CHAR *z_dir1,
102 int at_2, int i_next_at_2, S_CHAR *z_dir2 );
103
104 static int CompDble( const void *a1, const void *a2 );
105 static int Get2DTetrahedralAmbiguity( double at_coord[][3], int bAddExplicitNeighbor, int bFix2DstereoBorderCase );
106 static double triple_prod_and_min_abs_sine2(double at_coord[][3], double central_at_coord[], int bAddedExplicitNeighbor, double *min_sine, int *bAmbiguous);
107 static int are_4at_in_one_plane( double at_coord[][3], double min_sine);
108 static int bInpAtomHasRequirdNeigh ( inp_ATOM *at, int cur_at, int RequirdNeighType, int NumDbleBonds );
109 static int bIsSuitableHeteroInpAtom( inp_ATOM *at );
110 static int bIsOxide( inp_ATOM *at, int cur_at );
111 static int half_stereo_bond_parity( inp_ATOM *at, int cur_at, inp_ATOM *at_removed_H, int num_removed_H, S_CHAR *z_dir,
112 int bPointedEdgeStereo, int vABParityUnknown );
113 static int get_allowed_stereo_bond_type( int bond_type );
114 static int can_be_a_stereo_bond_with_isotopic_H( inp_ATOM *at, int cur_at, INCHI_MODE nMode );
115 static int half_stereo_bond_action( int nParity, int bUnknown, int bIsotopic, int vABParityUnknown );
116 static int set_stereo_bonds_parity( sp_ATOM *out_at, inp_ATOM *at, int at_1, inp_ATOM *at_removed_H, int num_removed_H,
117 INCHI_MODE nMode, QUEUE *q, AT_RANK *nAtomLevel,
118 S_CHAR *cSource, AT_RANK min_sb_ring_size,
119 int bPointedEdgeStereo, int vABParityUnknown );
120 static int can_be_a_stereo_atom_with_isotopic_H( inp_ATOM *at, int cur_at, int bPointedEdgeStereo );
121 static int set_stereo_atom_parity( sp_ATOM *out_at, inp_ATOM *at, int cur_at, inp_ATOM *at_removed_H, int num_removed_H,
122 int bPointedEdgeStereo, int vABParityUnknown );
123 /*
124 int set_stereo_parity( inp_ATOM* at, sp_ATOM* at_output, int num_at, int num_removed_H,
125 int *nMaxNumStereoAtoms, int *nMaxNumStereoBonds, INCHI_MODE nMode, int bPointedEdgeStereo, vABParityUnknown );
126 int get_opposite_sb_atom( inp_ATOM *at, int cur_atom, int icur2nxt, int *pnxt_atom, int *pinxt2cur, int *pinxt_sb_parity_ord );
127 */
128 int ReconcileCmlIncidentBondParities( inp_ATOM *at, int cur_atom, int prev_atom, S_CHAR *visited, int bDisconnected );
129 int comp_AT_NUMB( const void* a1, const void* a2);
130 int GetHalfStereobond0DParity( inp_ATOM *at, int cur_at, AT_NUMB nSbNeighOrigAtNumb[], int nNumExplictAttachments, int bond_parity, int nFlag );
131 int GetStereocenter0DParity( inp_ATOM *at, int cur_at, int j1, AT_NUMB nSbNeighOrigAtNumb[], int nFlag );
132 int GetSbNeighOrigAtNumb( inp_ATOM *at, int cur_at, inp_ATOM *at_removed_H, int num_removed_H, AT_NUMB nSbNeighOrigAtNumb[]);
133 int FixSb0DParities( inp_ATOM *at, /* inp_ATOM *at_removed_H, int num_removed_H,*/ int chain_length,
134 int at_1, int i_next_at_1, S_CHAR z_dir1[],
135 int at_2, int i_next_at_2, S_CHAR z_dir2[],
136 int *pparity1, int *pparity2 );
137
138 /******************************************************************/
139
140
141 static double *pDoubleForSort;
142
143 /**********************************************************************************/
comp_AT_NUMB(const void * a1,const void * a2)144 int comp_AT_NUMB( const void* a1, const void* a2)
145 {
146 return (int)*(const AT_NUMB*)a1 - (int)*(const AT_NUMB*)a2;
147 }
148 /******************************************************************/
get_z_coord(inp_ATOM * at,int cur_atom,int neigh_no,int * nType,int bPointedEdgeStereo)149 double get_z_coord( inp_ATOM* at, int cur_atom, int neigh_no, int *nType, int bPointedEdgeStereo )
150 {
151 int stereo_value = at[cur_atom].bond_stereo[neigh_no];
152 int stereo_type = abs( stereo_value );
153 int neigh = (int)at[cur_atom].neighbor[neigh_no];
154 double z = at[neigh].z - at[cur_atom].z;
155 int bFlat;
156
157 if ( (bFlat = (fabs(z) < ZERO_LENGTH)) ) {
158 int i;
159 for ( i = 0; i < at[cur_atom].valence; i ++ ) {
160 if ( fabs(at[cur_atom].z - at[(int)at[cur_atom].neighbor[i]].z) > ZERO_LENGTH ) {
161 bFlat = 0;
162 break;
163 }
164 }
165 }
166
167 if ( bFlat ) {
168 if ( !bPointedEdgeStereo || bPointedEdgeStereo * stereo_value >= 0 ) {
169 /* bPointedEdgeStereo > 0: define stereo from pointed end of the stereo bond only */
170 /* bPointedEdgeStereo < 0: define stereo from wide end of the stereo bond only (case of removed H) */
171 switch( stereo_type ) {
172 /* 1=Up (solid triangle), 6=Down (Dashed triangle), 4=Either (zigzag triangle) */
173 case 0: /* No stereo */
174 *nType = ZTYPE_NONE;
175 break;
176 case STEREO_SNGL_UP: /* 1= Up */
177 *nType = ZTYPE_UP;
178 break;
179 case STEREO_SNGL_EITHER: /* 4 = Either */
180 *nType = ZTYPE_EITHER;
181 break;
182 case STEREO_SNGL_DOWN: /* 6 = Down */
183 *nType = ZTYPE_DOWN;
184 break;
185 default:
186 *nType = ZTYPE_NONE; /* ignore unexpected values */
187 }
188 if ( stereo_value < 0 && (*nType == ZTYPE_DOWN || *nType == ZTYPE_UP) )
189 *nType = -*nType;
190 } else {
191 *nType = ZTYPE_NONE; /* no stereo */
192 }
193 } else
194 if ( stereo_type == STEREO_SNGL_EITHER &&
195 ( !bPointedEdgeStereo || bPointedEdgeStereo * stereo_value >= 0 ) ) {
196 *nType = ZTYPE_EITHER;
197 } else {
198 *nType = ZTYPE_3D;
199 }
200 return z;
201 }
202 /******************************************************************/
len3(const double c[])203 double len3( const double c[] )
204 {
205 return sqrt( c[0]*c[0] + c[1]*c[1] + c[2]*c[2] );
206 }
207 /******************************************************************/
len2(const double c[])208 double len2( const double c[] )
209 {
210 return sqrt( c[0]*c[0] + c[1]*c[1] );
211 }
212 /******************************************************************/
diff3(const double a[],const double b[],double result[])213 double* diff3( const double a[], const double b[], double result[] )
214 {
215
216 result[0] = a[0] - b[0];
217 result[1] = a[1] - b[1];
218 result[2] = a[2] - b[2];
219
220 return result;
221 }
222 /******************************************************************/
add3(const double a[],const double b[],double result[])223 double* add3( const double a[], const double b[], double result[] )
224 {
225 result[0] = a[0] + b[0];
226 result[1] = a[1] + b[1];
227 result[2] = a[2] + b[2];
228
229 return result;
230 }
231 /******************************************************************/
mult3(const double a[],double b,double result[])232 double* mult3( const double a[], double b, double result[] )
233 {
234 result[0] = a[0] * b;
235 result[1] = a[1] * b;
236 result[2] = a[2] * b;
237
238 return result;
239 }
240 /*************************************************************/
copy3(const double a[],double result[])241 double* copy3( const double a[], double result[] )
242 {
243 result[0] = a[0];
244 result[1] = a[1];
245 result[2] = a[2];
246
247 return result;
248 }
249 /*************************************************************/
change_sign3(const double a[],double result[])250 double* change_sign3( const double a[], double result[] )
251 {
252 result[0] = -a[0];
253 result[1] = -a[1];
254 result[2] = -a[2];
255
256 return result;
257 }
258 /*************************************************************/
dot_prod3(const double a[],const double b[])259 double dot_prod3( const double a[], const double b[] )
260 {
261 return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
262 }
263 /*************************************************************/
dot_prodchar3(const S_CHAR a[],const S_CHAR b[])264 int dot_prodchar3( const S_CHAR a[], const S_CHAR b[] )
265 {
266 int prod = ((int)a[0]*(int)b[0] + (int)a[1]*(int)b[1] + (int)a[2]*(int)b[2])/100;
267 if ( prod > 100 )
268 prod = 100;
269 else
270 if ( prod < -100 )
271 prod = -100;
272 return prod;
273 }
274 /*************************************************************/
cross_prod3(const double a[],const double b[],double result[])275 double* cross_prod3( const double a[], const double b[], double result[] )
276 {
277 double tmp[3];
278
279 tmp[0] = (a[1]*b[2]-a[2]*b[1]);
280 tmp[1] = -(a[0]*b[2]-a[2]*b[0]);
281 tmp[2] = (a[0]*b[1]-a[1]*b[0]);
282
283 result[0] = tmp[0];
284 result[1] = tmp[1];
285 result[2] = tmp[2];
286
287 return result;
288 }
289 /*************************************************************/
triple_prod(double a[],double b[],double c[],double * sine_value)290 double triple_prod( double a[], double b[], double c[], double *sine_value )
291 {
292 double ab[3], dot_prod_ab_c, abs_c, abs_ab;
293 cross_prod3( a, b, ab );
294 /* ab[0] = (a[1]*b[2]-a[2]*b[1]); */
295 /* ab[1] = -(a[0]*b[2]-a[2]*b[0]); */
296 /* ab[2] = (a[0]*b[1]-a[1]*b[0]); */
297 dot_prod_ab_c = dot_prod3( ab, c );
298 /* dot_prod_ab_c = ab[0]*c[0] + ab[1]*c[1] + ab[2]*c[2]; */
299 if ( sine_value ) {
300 abs_c = len3( c );
301 /* abs_c = sqrt( c[0]*c[0] + c[1]*c[1] + c[2]*c[2] ); */
302 abs_ab = len3( ab );
303 /* abs_ab = sqrt( ab[0]*ab[0] + ab[1]*ab[1] + ab[2]*ab[2] ); */
304
305 if ( abs_c > 1.e-7 /* otherwise c has zero length */ && abs_ab > 1.e-7 /* otherwise a is parallel to b*/ ) {
306 *sine_value = MPY_SINE * dot_prod_ab_c / ( abs_c * abs_ab);
307 /* *sine_value = dot_prod_ab_c / ( abs_c * abs_ab); */
308 } else {
309 *sine_value = 0.0;
310 }
311 }
312 return dot_prod_ab_c;
313 }
314 /*************************************************************/
CompDble(const void * a1,const void * a2)315 int CompDble( const void *a1, const void *a2 )
316 {
317 double diff = pDoubleForSort[*(const int*)a1] - pDoubleForSort[*(const int*)a2];
318 if ( diff > 0.0 )
319 return 1;
320 if ( diff < 0.0 )
321 return -1;
322 return 0;
323 }
324 /*************************************************************/
325 #define T2D_OKAY 1
326 #define T2D_WARN 2
327 #define T2D_UNDF 4
Get2DTetrahedralAmbiguity(double at_coord[][3],int bAddExplicitNeighbor,int bFix2DstereoBorderCase)328 int Get2DTetrahedralAmbiguity( double at_coord[][3], int bAddExplicitNeighbor, int bFix2DstereoBorderCase )
329 {
330 /* const double one_pi = 2.0*atan2(1.0 , 0.0 ); */
331 const double one_pi = 3.14159265358979323846; /* M_PI */
332 const double two_pi = 2.0*one_pi;
333 const double dAngleAndPiMaxDiff = 2.0*atan2(1.0, sqrt(7.0)); /* min sine between 2 InPlane bonds */
334 int nBondType[MAX_NUM_STEREO_ATOM_NEIGH], nBondOrder[MAX_NUM_STEREO_ATOM_NEIGH];
335 double dBondDirection[MAX_NUM_STEREO_ATOM_NEIGH];
336 volatile double dAngle, dAlpha, dLimit, dBisector;
337 /* 2010-02-10 added 'volatile': workaround ensuring proper behavior for gcc 32-bit */
338 /* cml-enabled compiles at >=O1 for SID484922 and alike (both lin&win had problems) */
339 int nNumNeigh = MAX_NUM_STEREO_ATOM_NEIGH - (bAddExplicitNeighbor != 0);
340 int i, num_Up, num_Dn, bPrev_Up, cur_len_Up, cur_first_Up, len_Up, first_Up;
341 int ret=0;
342
343 for ( i = 0, num_Up = num_Dn = 0; i < nNumNeigh; i ++ )
344 {
345 dAngle = atan2( at_coord[i][1], at_coord[i][0] ); /* range from -pi to +pi */
346 if ( dAngle < 0.0 ) {
347 dAngle += two_pi;
348 }
349 dBondDirection[i] = dAngle;
350 nBondType[i] = (at_coord[i][2] > 0.0)? 1 : (at_coord[i][2] < 0.0)? -1 : 0; /* z-coord sign */
351 if ( nBondType[i] > 0 ) {
352 num_Up ++;
353 } else
354 if ( nBondType[i] < 0 ) {
355 num_Dn ++;
356 }
357 nBondOrder[i] = i;
358 }
359 if ( num_Up < num_Dn ) {
360 for ( i = 0; i < nNumNeigh; i ++ ) {
361 nBondType[i] = -nBondType[i];
362 }
363 inchi_swap( (char*)&num_Dn, (char*)&num_Up, sizeof(num_Dn) );
364 }
365 if ( !num_Up ) {
366 return T2D_UNDF;
367 }
368
369 /* sort according to the bond orientations */
370 pDoubleForSort = dBondDirection;
371 insertions_sort( nBondOrder, (unsigned) nNumNeigh, sizeof(nBondOrder[0]), CompDble );
372
373 /* find the longest contiguous sequence of Up bonds */
374 if ( num_Up == nNumNeigh ) {
375 /* all bonds are Up */
376 len_Up = cur_len_Up = nNumNeigh; /* added cur_len_Up initialization 1/8/2002 */
377 first_Up = 0;
378 } else {
379 /* at least one bond is not Up */
380 cur_len_Up = len_Up = bPrev_Up = 0;
381 /* prev. cycle header version ---
382 for ( i = 0; 1; i ++ ) {
383 if ( i >= nNumNeigh && !bPrev_Up ) {
384 break;
385 }
386 ----------} */
387 /* look at all bonds and continue (circle therough the beginning) as long as the current bond is Up */
388 for ( i = 0; i < nNumNeigh || bPrev_Up; i ++ ) {
389 if ( nBondType[nBondOrder[i % nNumNeigh]] > 0 ) {
390 if ( bPrev_Up ) {
391 cur_len_Up ++; /* uncrement number of Up bonds in current contiguous sequence of them */
392 } else {
393 bPrev_Up = 1; /* start new contiguous sequence of Up bonds */
394 cur_len_Up = 1;
395 cur_first_Up = i % nNumNeigh;
396 }
397 } else
398 if ( bPrev_Up ) { /* end of contiguous sequence of Up bonds */
399 if ( cur_len_Up > len_Up ) {
400 first_Up = cur_first_Up; /* store the sequence because it is longer than the ptrvious one */
401 len_Up = cur_len_Up;
402 }
403 bPrev_Up = 0;
404 }
405 }
406 }
407 #if ( FIX_2D_STEREO_BORDER_CASE == 1 )
408 /* check if the bonds with ordering numbers first_Up+len_Up and first_Up+len_Up+1 */
409 /* have identical angles. In this case switch their order to enlarge the Up sequence */
410 #define ZERO_ANGLE 0.000001
411 if ( nNumNeigh - len_Up >= 2 ) {
412 int next1, next2;
413 for ( i = 1; i < nNumNeigh - len_Up; i ++ ) {
414 next2 = (first_Up+len_Up + i) % nNumNeigh; /* the 2nd after Up sequence */
415 if ( nBondType[nBondOrder[next2]] > 0 ) {
416 next1 = (first_Up+len_Up) % nNumNeigh; /* the 1st after Up sequence */
417 dAngle = dBondDirection[nBondOrder[next1]] - dBondDirection[nBondOrder[next2]];
418 if ( fabs(dAngle) < ZERO_ANGLE ) {
419 inchi_swap( (char*)&nBondOrder[next1], (char*)&nBondOrder[next2], sizeof(nBondOrder[0]) );
420 len_Up ++;
421 break;
422 }
423 }
424 }
425 }
426 /* check whether the not-Up bond (located before the found first-Up) has */
427 /* same angle as the Up bond that precedes this not-Up bond */
428 if ( nNumNeigh - len_Up >= 2 ) {
429 int next1, next2;
430 for ( i = 1; i < nNumNeigh - len_Up; i ++ ) {
431 next2 = (first_Up+nNumNeigh - i - 1 ) % nNumNeigh; /* the 2nd before Up sequence */
432 if ( nBondType[nBondOrder[next2]] > 0 ) {
433 next1 = (first_Up+nNumNeigh-1) % nNumNeigh; /* the 1st before Up sequence */
434 dAngle = dBondDirection[nBondOrder[next1]] - dBondDirection[nBondOrder[next2]];
435 if ( fabs(dAngle) < ZERO_ANGLE ) {
436 inchi_swap( (char*)&nBondOrder[next1], (char*)&nBondOrder[next2], sizeof(nBondOrder[0]) );
437 first_Up = next1;
438 len_Up ++;
439 break;
440 }
441 }
442 }
443 }
444 #else
445 if ( bFix2DstereoBorderCase ) {
446 /* check if the bonds with ordering numbers first_Up+len_Up and first_Up+len_Up+1 */
447 /* have identical angles. In this case switch their order to enlarge the Up sequence */
448 #define ZERO_ANGLE 0.000001
449 if ( nNumNeigh - len_Up >= 2 ) {
450 int next1, next2;
451 for ( i = 1; i < nNumNeigh - len_Up; i ++ ) {
452 next2 = (first_Up+len_Up + i) % nNumNeigh; /* the 2nd after Up sequence */
453 if ( nBondType[nBondOrder[next2]] > 0 ) {
454 next1 = (first_Up+len_Up) % nNumNeigh; /* the 1st after Up sequence */
455 dAngle = dBondDirection[nBondOrder[next1]] - dBondDirection[nBondOrder[next2]];
456 if ( fabs(dAngle) < ZERO_ANGLE ) {
457 inchi_swap( (char*)&nBondOrder[next1], (char*)&nBondOrder[next2], sizeof(nBondOrder[0]) );
458 len_Up ++;
459 break;
460 }
461 }
462 }
463 }
464 /* check whether the not-Up bond (located before the found first-Up) has */
465 /* same angle as the Up bond that precedes this not-Up bond */
466 if ( nNumNeigh - len_Up >= 2 ) {
467 int next1, next2;
468 for ( i = 1; i < nNumNeigh - len_Up; i ++ ) {
469 next2 = (first_Up+nNumNeigh - i - 1 ) % nNumNeigh; /* the 2nd before Up sequence */
470 if ( nBondType[nBondOrder[next2]] > 0 ) {
471 next1 = (first_Up+nNumNeigh-1) % nNumNeigh; /* the 1st before Up sequence */
472 dAngle = dBondDirection[nBondOrder[next1]] - dBondDirection[nBondOrder[next2]];
473 if ( fabs(dAngle) < ZERO_ANGLE ) {
474 inchi_swap( (char*)&nBondOrder[next1], (char*)&nBondOrder[next2], sizeof(nBondOrder[0]) );
475 first_Up = next1;
476 len_Up ++;
477 break;
478 }
479 }
480 }
481 }
482 }
483 #endif
484 /* Turn all the bonds around the center so that */
485 /* the 1st Up bond has zero radian direction */
486 dAlpha = dBondDirection[nBondOrder[first_Up]];
487 for ( i = 0; i < nNumNeigh; i ++ ) {
488 if ( i == nBondOrder[first_Up] ) {
489 dBondDirection[i] = 0.0;
490 } else {
491 dAngle = dBondDirection[i] - dAlpha;
492 if ( dAngle < 0.0 ) {
493 dAngle += two_pi;
494 }
495 dBondDirection[i] = dAngle;
496 }
497 }
498
499 /********************************************************
500 * Process particular cases
501 ********************************************************/
502
503
504 if ( nNumNeigh == 3 ) /************************ 3 bonds ************************/
505 {
506
507 switch( num_Up )
508 {
509
510
511
512 case 0: /* 0 Up */
513 return T2D_UNDF;
514
515
516
517
518 case 1: /* 1 Up */
519 if ( num_Dn )
520 {
521 #ifdef _DEBUG
522 if ( num_Dn != 1 ) /* debug only */
523 return -1;
524 #endif
525 ret = (T2D_UNDF | T2D_WARN);
526 }
527 else
528 {
529 dAngle = dBondDirection[nBondOrder[(first_Up + 2) % nNumNeigh]] -
530 dBondDirection[nBondOrder[(first_Up + 1) % nNumNeigh]];
531
532 if ( dAngle < 0.0 )
533 dAngle += two_pi;
534 if ( dAngle - one_pi < -MIN_ANGLE || dAngle - one_pi > MIN_ANGLE )
535 {
536 ret = T2D_OKAY;
537 }
538 else
539 {
540 ret = (T2D_UNDF | T2D_WARN);
541 }
542 }
543 break;
544
545
546
547
548 case 2: /* 2 Up */
549 if ( num_Dn )
550 {
551 dAlpha = dBondDirection[nBondOrder[(first_Up + 1) % nNumNeigh]] -
552 dBondDirection[nBondOrder[(first_Up ) % nNumNeigh]];
553
554 if ( dAlpha < 0.0 )
555 dAlpha += two_pi;
556
557 if ( dAlpha > one_pi - MIN_ANGLE )
558 {
559 ret = T2D_OKAY;
560 }
561 else if ( dAlpha < two_pi / 3.0 - MIN_ANGLE )
562 {
563 ret = (T2D_UNDF | T2D_WARN);
564 }
565 else
566 {
567 /* angle between 2 Up bonds is between 120 and 180 degrees */
568 /* direction of the (Alpha angle bisector) + 180 degrees */
569 dBisector = dBondDirection[nBondOrder[(first_Up ) % nNumNeigh]];
570 dBisector+= dBondDirection[nBondOrder[(first_Up + 1 ) % nNumNeigh]];
571 dBisector/= 2.0;
572 dBisector-= one_pi;
573 if ( dBisector < 0.0 )
574 {
575 dBisector += two_pi;
576 }
577 if ( dAlpha < two_pi / 3.0 + MIN_ANGLE )
578 {
579 /* dAlpha is inside ( 2pi/3 - eps, 2pi/3 + eps ) interval */
580 dLimit = MIN_ANGLE * 3.0 / 2.0;
581 }
582 else
583 {
584 dLimit = dAlpha * 3.0 / 2.0 - one_pi;
585 }
586
587 dAngle = dBondDirection[nBondOrder[(first_Up + 2 ) % nNumNeigh]];
588
589 if ( dBisector - dAngle < -dLimit ||
590 dBisector - dAngle > dLimit )
591 {
592 ret = (T2D_UNDF | T2D_WARN);
593 }
594 else
595 {
596 ret = T2D_OKAY;
597 }
598 }
599 } /* if ( num_Dn ) */
600 else
601 {
602 ret = T2D_OKAY;
603 }
604 break;
605
606
607
608 case 3: /* 3 Up */
609 ret = T2D_OKAY;
610 break;
611
612
613 default:/* other Up */
614 return -1;
615
616 } /* eof switch( num_Up ) at nNumNeigh == 3 */
617
618 }
619
620
621 else if ( nNumNeigh == 4) /******************************* 4 bonds ********************/
622 {
623 switch( num_Up )
624 {
625
626 case 0: /* 0 Up */
627 return T2D_UNDF;
628
629
630 case 1: /* 1 Up */
631 if ( num_Dn )
632 {
633 if ( nBondType[nBondOrder[(first_Up + 2) % nNumNeigh]] < 0 )
634 {
635 /*
636 * Up, In Plane, Dn, In Plane. Undefined if angle between
637 * two In Plane bonds is wuthin pi +/- 2*arcsine(1/sqrt(8)) interval
638 * That is, 138.5 to 221.4 degrees; for certainty the interval is
639 * increased by 5.7 degrees at each end to
640 * 134.8 to 227.1 degrees
641 */
642 dAngle = dBondDirection[nBondOrder[(first_Up + 3) % nNumNeigh]] -
643 dBondDirection[nBondOrder[(first_Up + 1) % nNumNeigh]];
644 if ( dAngle < 0.0 ) {
645 dAngle += two_pi;
646 }
647 if ( fabs( dAngle - one_pi ) < dAngleAndPiMaxDiff + MIN_ANGLE ) {
648 ret = (T2D_UNDF | T2D_WARN);
649 }
650 else
651 {
652 ret = T2D_OKAY;
653 }
654 }
655 else
656 {
657 ret = T2D_OKAY;
658 }
659 #ifdef _DEBUG
660 if ( num_Dn != 1 ) /* debug only */
661 return -1;
662 #endif
663 }
664 else
665 {
666 ret = T2D_OKAY;
667 dAngle = dBondDirection[nBondOrder[(first_Up + 3) % nNumNeigh]] -
668 dBondDirection[nBondOrder[(first_Up + 1) % nNumNeigh]];
669 if ( dAngle < 0.0 )
670 {
671 dAngle += two_pi;
672 }
673 if ( dAngle < one_pi - MIN_ANGLE )
674 {
675 ret |= T2D_WARN;
676 }
677 }
678 break;
679
680
681 case 2: /* 2 Up */
682 #if ( FIX_2D_STEREO_BORDER_CASE == 1 )
683 if ( len_Up == 1 )
684 {
685 ret = T2D_OKAY;
686 }
687 else
688 {
689 dAngle = dBondDirection[nBondOrder[(first_Up + 3) % nNumNeigh]] -
690 dBondDirection[nBondOrder[(first_Up + 0) % nNumNeigh]];
691 dAngle = fabs(two_pi - dAngle);
692 dAlpha = dBondDirection[nBondOrder[(first_Up + 2) % nNumNeigh]] -
693 dBondDirection[nBondOrder[(first_Up + 1) % nNumNeigh]];
694 dAlpha = fabs(dAlpha);
695 if ( dAngle < 2.0 * ZERO_ANGLE && dAlpha > MIN_ANGLE ||
696 dAlpha < 2.0 * ZERO_ANGLE && dAngle > MIN_ANGLE )
697 {
698 ret = (T2D_OKAY | T2D_WARN);
699 }
700 else
701 {
702 ret = (T2D_UNDF | T2D_WARN);
703 }
704 }
705 #else
706 if ( bFix2DstereoBorderCase )
707 {
708 /* bug fix */
709 if ( len_Up == 1 )
710 {
711 ret = T2D_OKAY;
712 }
713 else
714 {
715 dAngle = dBondDirection[nBondOrder[(first_Up + 3) % nNumNeigh]] -
716 dBondDirection[nBondOrder[(first_Up + 0) % nNumNeigh]];
717 dAngle = fabs(two_pi - dAngle);
718 dAlpha = dBondDirection[nBondOrder[(first_Up + 2) % nNumNeigh]] -
719 dBondDirection[nBondOrder[(first_Up + 1) % nNumNeigh]];
720 dAlpha = fabs(dAlpha);
721 if ( (dAngle < 2.0 * ZERO_ANGLE && dAlpha > MIN_ANGLE) ||
722 (dAlpha < 2.0 * ZERO_ANGLE && dAngle > MIN_ANGLE) )
723 {
724 ret = (T2D_OKAY | T2D_WARN);
725 }
726 else
727 {
728 ret = (T2D_UNDF | T2D_WARN);
729 }
730 }
731 }
732 else
733 {
734 /* original InChI v. 1 bug */
735 if ( cur_len_Up == 1 )
736 {
737 ret = T2D_OKAY;
738 }
739 else
740 {
741 ret = (T2D_UNDF | T2D_WARN);
742 }
743 }
744 #endif
745 break;
746
747
748 case 3: /* 3 Up */
749 ret = T2D_OKAY;
750 dAngle = dBondDirection[nBondOrder[(first_Up + 2) % nNumNeigh]] -
751 dBondDirection[nBondOrder[(first_Up + 0) % nNumNeigh]];
752 if ( dAngle < 0.0 )
753 {
754 dAngle += two_pi;
755 }
756 if ( dAngle < one_pi - MIN_ANGLE )
757 {
758 ret |= T2D_WARN;
759 }
760 break;
761
762 case 4: /* 4 Up */
763 ret = (T2D_UNDF | T2D_WARN);
764 break;
765
766 default:/* other Up */
767 return -1; /* program error */
768
769 } /* eof switch( num_Up ) at nNumNeigh == 4 */
770
771 if ( ret == T2D_OKAY )
772 {
773 /* check whether all bonds are inside a less than 180 degrees sector */
774 for ( i = 0; i < nNumNeigh; i ++ )
775 {
776 dAngle = dBondDirection[nBondOrder[(i + nNumNeigh - 1) % nNumNeigh]] -
777 dBondDirection[nBondOrder[ i % nNumNeigh]];
778 if ( dAngle < 0.0 )
779 {
780 dAngle += two_pi;
781 }
782 if ( dAngle < one_pi - MIN_ANGLE )
783 {
784 ret |= T2D_WARN;
785 break;
786 }
787 }
788 }
789
790 } /* eof nNumNeigh == 4 */
791
792 else /*************************** number of bonds != 3 or 4 ******************/
793 {
794
795 return -1; /* error */
796 }
797
798
799 return ret;
800
801 }
802
803 /*************************************************************/
triple_prod_and_min_abs_sine2(double at_coord[][3],double central_at_coord[],int bAddedExplicitNeighbor,double * min_sine,int * bAmbiguous)804 double triple_prod_and_min_abs_sine2(double at_coord[][3], double central_at_coord[], int bAddedExplicitNeighbor, double *min_sine, int *bAmbiguous)
805 {
806 double min_sine_value=9999.0, sine_value, min_edge_len, max_edge_len, min_edge_len_NoExplNeigh, max_edge_len_NoExplNeigh;
807 double s0, s1, s2, s3, e01, e02, e03, e12, e13, e23, tmp[3], e[3][3];
808 double prod, ret, central_prod[4];
809 int bLongEdges;
810
811 if ( !min_sine ) {
812 return triple_prod( at_coord[0], at_coord[1], at_coord[2], NULL );
813 }
814
815 ret = triple_prod( at_coord[0], at_coord[1], at_coord[2], &sine_value );
816 sine_value = MPY_SINE * fabs( sine_value );
817
818 diff3( at_coord[1], at_coord[0], e[2] );
819 diff3( at_coord[0], at_coord[2], e[1] );
820 diff3( at_coord[2], at_coord[1], e[0] );
821
822 /* lengths of the 6 edges of the tetrahedra */
823 e03 = len3( at_coord[0] ); /* 1 */
824 e13 = len3( at_coord[1] );
825 e23 = len3( at_coord[2] ); /* includes added neighbor if bAddedExplicitNeighbor*/
826 e02 = len3( e[1] ); /* includes added neighbor if bAddedExplicitNeighbor*/
827 e12 = len3( e[0] ); /* includes added neighbor if bAddedExplicitNeighbor*/
828 e01 = len3( e[2] );
829
830 /* min & max edge length */
831 max_edge_len =
832 min_edge_len = e03;
833
834 if ( min_edge_len > e13 )
835 min_edge_len = e13;
836 if ( min_edge_len > e01 )
837 min_edge_len = e01;
838 min_edge_len_NoExplNeigh = min_edge_len;
839
840 if ( min_edge_len > e23 )
841 min_edge_len = e23;
842 if ( min_edge_len > e02 )
843 min_edge_len = e02;
844 if ( min_edge_len > e12 )
845 min_edge_len = e12;
846
847 if ( max_edge_len < e13 )
848 max_edge_len = e13;
849 if ( max_edge_len < e01 )
850 max_edge_len = e01;
851 max_edge_len_NoExplNeigh = max_edge_len;
852
853 if ( max_edge_len < e23 )
854 max_edge_len = e23;
855 if ( max_edge_len < e02 )
856 max_edge_len = e02;
857 if ( max_edge_len < e12 )
858 max_edge_len = e12;
859
860 if ( !bAddedExplicitNeighbor ) {
861 min_edge_len_NoExplNeigh = min_edge_len;
862 max_edge_len_NoExplNeigh = max_edge_len;
863 }
864
865 bLongEdges = bAddedExplicitNeighbor?
866 ( max_edge_len_NoExplNeigh < MAX_EDGE_RATIO * min_edge_len_NoExplNeigh ) :
867 ( max_edge_len < MAX_EDGE_RATIO * min_edge_len );
868
869 if ( sine_value > MIN_SINE && ( min_sine || bAmbiguous ) ) {
870 if ( min_sine ) {
871 prod = fabs( ret );
872 /* tetrahedra height = volume(prod) / area of a plane(cross_prod) */
873 /* (instead of a tetrahedra calculate parallelogram/parallelepiped area/volume) */
874
875 /* 4 heights from each of the 4 vertices to the opposite plane */
876 s0 = prod / len3( cross_prod3( at_coord[1], at_coord[2], tmp ) );
877 s1 = prod / len3( cross_prod3( at_coord[0], at_coord[2], tmp ) );
878 s2 = prod / len3( cross_prod3( at_coord[0], at_coord[1], tmp ) );
879 s3 = prod / len3( cross_prod3( e[0], e[1], tmp ) );
880 /* abs. value of a sine of an angle between each tetrahedra edge and plane */
881 /* sine = height / edge length */
882 if ( (sine_value = s0/e01) < min_sine_value )
883 min_sine_value = sine_value;
884 if ( (sine_value = s0/e02) < min_sine_value )
885 min_sine_value = sine_value;
886 if ( (sine_value = s0/e03) < min_sine_value )
887 min_sine_value = sine_value;
888
889 if ( (sine_value = s1/e01) < min_sine_value )
890 min_sine_value = sine_value;
891 if ( (sine_value = s1/e12) < min_sine_value )
892 min_sine_value = sine_value;
893 if ( (sine_value = s1/e13) < min_sine_value )
894 min_sine_value = sine_value;
895
896 if ( (sine_value = s2/e02) < min_sine_value )
897 min_sine_value = sine_value;
898 if ( (sine_value = s2/e12) < min_sine_value )
899 min_sine_value = sine_value;
900 if ( (sine_value = s2/e23) < min_sine_value )
901 min_sine_value = sine_value;
902
903 if ( (sine_value = s3/e03) < min_sine_value )
904 min_sine_value = sine_value;
905 if ( (sine_value = s3/e13) < min_sine_value )
906 min_sine_value = sine_value;
907 if ( (sine_value = s3/e23) < min_sine_value )
908 min_sine_value = sine_value;
909 /* actually use triple sine */
910 *min_sine = sine_value = MPY_SINE * min_sine_value;
911 }
912
913 if ( bAmbiguous && sine_value >= MIN_SINE ) {
914 /* check whether the central atom is outside the tetrahedra (0,0,0), at_coord[0,1,2] */
915 /* compare the tetrahedra volume and the volume of a tetrahedra having central_at_coord[] vertex */
916 int i;
917 diff3( central_at_coord, at_coord[0], tmp );
918 central_prod[0] = triple_prod( at_coord[0], at_coord[1], central_at_coord, NULL );
919 central_prod[1] = triple_prod( at_coord[1], at_coord[2], central_at_coord, NULL );
920 central_prod[2] = triple_prod( at_coord[2], at_coord[0], central_at_coord, NULL );
921 central_prod[3] = triple_prod( e[2], e[1], tmp, NULL );
922 for ( i = 0; i <= 3; i ++ ) {
923 if ( central_prod[i] / ret < -MIN_SINE_OUTSIDE ) {
924 *bAmbiguous |= AMBIGUOUS_STEREO;
925 break;
926 }
927 }
928 }
929 #if ( STEREO_CENTER_BONDS_NORM == 1 )
930
931 if ( bLongEdges && !bAddedExplicitNeighbor && max_edge_len >= MIN_LEN_STRAIGHT ) {
932 /* possible planar tetragon */
933 if ( sine_value < MIN_SINE_SQUARE ) {
934 *min_sine = MIN_SINE / 2.0; /* force parity to be undefined */
935 if ( bAmbiguous && !*bAmbiguous ) {
936 *bAmbiguous |= AMBIGUOUS_STEREO;
937 }
938 }
939 }
940
941 if ( bLongEdges && sine_value < MIN_SINE_SQUARE && sine_value < MIN_SINE_EDGE * min_edge_len_NoExplNeigh ) {
942 *min_sine = MIN_SINE / 2.0; /* force parity to be undefined */
943 if ( bAmbiguous && !*bAmbiguous ) {
944 *bAmbiguous |= AMBIGUOUS_STEREO;
945 }
946 }
947 #endif
948
949 } else
950 if ( min_sine ) {
951 *min_sine = sine_value;
952 }
953
954 return ret;
955 }
956 /*************************************************************/
triple_prod_and_min_abs_sine(double at_coord[][3],double * min_sine)957 double triple_prod_and_min_abs_sine(double at_coord[][3], double *min_sine)
958 {
959 double min_sine_value=9999.0, sine_value;
960 double prod=0.0;
961
962 if ( !min_sine ) {
963 return triple_prod( at_coord[0], at_coord[1], at_coord[2], NULL );
964 }
965
966 prod = triple_prod( at_coord[0], at_coord[1], at_coord[2], &sine_value );
967 sine_value = fabs( sine_value );
968 min_sine_value = inchi_min( min_sine_value, sine_value );
969
970 prod = triple_prod( at_coord[1], at_coord[2], at_coord[0], &sine_value );
971 sine_value = fabs( sine_value );
972 min_sine_value = inchi_min( min_sine_value, sine_value );
973
974 prod = triple_prod( at_coord[2], at_coord[0], at_coord[1], &sine_value );
975 sine_value = fabs( sine_value );
976 min_sine_value = inchi_min( min_sine_value, sine_value );
977
978 *min_sine = min_sine_value;
979
980 return prod;
981 }
982 /*************************************************************/
983 /* Find if point (0,0,0)a and 3 atoms are in one plane */
are_3_vect_in_one_plane(double at_coord[][3],double min_sine)984 int are_3_vect_in_one_plane( double at_coord[][3], double min_sine)
985 {
986 double actual_min_sine;
987 double prod;
988 prod = triple_prod_and_min_abs_sine( at_coord, &actual_min_sine);
989 return actual_min_sine <= min_sine;
990 }
991 /*************************************************************/
992 /* Find if 4 atoms are in one plane */
are_4at_in_one_plane(double at_coord[][3],double min_sine)993 int are_4at_in_one_plane( double at_coord[][3], double min_sine)
994 {
995 double actual_min_sine, min_actual_min_sine;
996 double coord[3][3], prod;
997 int i, k, j;
998 for ( k = 0; k < 4; k ++ ) { /* cycle added 4004-08-15 */
999 for ( i = j = 0; i < 4; i ++ ) {
1000 if ( i != k ) {
1001 diff3( at_coord[i], at_coord[k], coord[j] );
1002 j ++;
1003 }
1004 }
1005 prod = triple_prod_and_min_abs_sine( coord, &actual_min_sine);
1006 if ( !k || actual_min_sine < min_actual_min_sine ) {
1007 min_actual_min_sine = actual_min_sine;
1008 }
1009 }
1010 return min_actual_min_sine <= min_sine;
1011 }
1012 /*************************************************************/
triple_prod_char(inp_ATOM * at,int at_1,int i_next_at_1,S_CHAR * z_dir1,int at_2,int i_next_at_2,S_CHAR * z_dir2)1013 int triple_prod_char( inp_ATOM *at, int at_1, int i_next_at_1, S_CHAR *z_dir1,
1014 int at_2, int i_next_at_2, S_CHAR *z_dir2 )
1015 {
1016 inp_ATOM *at1, *at2;
1017 double pnt[3][3], len;
1018 int i;
1019 int ret = 0;
1020
1021 at1 = at + at_1;
1022 at2 = at + at[at_1].neighbor[i_next_at_1];
1023
1024 pnt[0][0] = at2->x - at1->x;
1025 pnt[0][1] = at2->y - at1->y;
1026 pnt[0][2] = at2->z - at1->z;
1027
1028 at2 = at + at_2;
1029 at1 = at + at[at_2].neighbor[i_next_at_2];
1030
1031 pnt[1][0] = at2->x - at1->x;
1032 pnt[1][1] = at2->y - at1->y;
1033 pnt[1][2] = at2->z - at1->z;
1034 /*
1035 * resultant pnt vector directions:
1036 *
1037 * pnt[0] pnt[1]
1038 *
1039 * [at_1]---->[...] [...]---->[at_2]
1040 *
1041 *
1042 * add3 below: (pnt[0] + pnt[1]) -> pnt[1]
1043 */
1044 add3( pnt[0], pnt[1], pnt[1] );
1045
1046
1047
1048 for ( i = 0; i < 3; i ++ ) {
1049 pnt[0][i] = (double)z_dir1[i];
1050 pnt[2][i] = (double)z_dir2[i];
1051 }
1052 for ( i = 0; i < 3; i ++ ) {
1053 len = len3( pnt[i] );
1054 if ( len < MIN_BOND_LEN ) {
1055 if ( i == 1 && (at[at_1].bUsed0DParity || at[at_2].bUsed0DParity) ) {
1056 pnt[i][0] = 0.0;
1057 pnt[i][1] = 1.0;
1058 pnt[i][2] = 0.0;
1059 len = 1.0; /* standard at_1-->at_2 vector coordinates in case of 0D allene */
1060 } else {
1061 goto exit_function; /* too short bond */
1062 }
1063 }
1064 mult3( pnt[i], 1.0/len, pnt[i] );
1065 }
1066 len = 100.0*triple_prod(pnt[0], pnt[1], pnt[2], NULL );
1067 /*
1068 * ^ pnt[0]
1069 * | The orientation on this diagram
1070 * | produces len = -100
1071 * [at_1]------>[at_2]
1072 * pnt[1] /
1073 * /
1074 * / pnt[2] (up from the plane)
1075 * v
1076 *
1077 * Note: len is invariant upon at_1 <--> at_2 transposition because
1078 * triple product changes sign upon pnt[0]<-->pnt[2] transposition and
1079 * triple product changes sign upon pnt[1]--> -pnt[1] change of direction:
1080 *
1081 * triple_prod(pnt[0], pnt[1], pnt[2], NULL ) =
1082 * triple_prod(pnt[2], -pnt[1], pnt[0], NULL )
1083 *
1084 */
1085
1086 ret = len >= 0.0? (int)(floor(len+0.5)) : -(int)(floor(0.5-len));
1087
1088 exit_function:
1089
1090 return ret;
1091 }
1092
1093
1094 /****************************************************************/
1095
1096 #if ( NEW_STEREOCENTER_CHECK == 1 ) /* { */
1097
1098 /********************************************************************************************/
bInpAtomHasRequirdNeigh(inp_ATOM * at,int cur_at,int RequirdNeighType,int NumDbleBonds)1099 int bInpAtomHasRequirdNeigh ( inp_ATOM *at, int cur_at, int RequirdNeighType, int NumDbleBonds )
1100 {
1101 /* RequirdNeighType:
1102 reqired neighbor types (bitmap):
1103 0 => any neighbors
1104 1 => no terminal hydrogen atom neighbors
1105 2 => no terminal -X and -XH together (don't care about -X, -XH bond type, charge, radical)
1106 (X = tautomeric endpoint atom)
1107 NumDbleBonds:
1108 if non-zero then allow double, alternating and tautomeric bonds
1109 */
1110 int i, j, ni, nj, bond_type, num_1s, num_mult, num_other;
1111
1112 if ( at[cur_at].endpoint ) { /* tautomeric endpoint cannot be a stereo center */
1113 return 0;
1114 }
1115
1116 if ( (1 & RequirdNeighType) && at[cur_at].num_H ) {
1117 return 0;
1118 }
1119
1120 if ( 2 & RequirdNeighType ) {
1121 for ( i = 0; i < at[cur_at].valence; i ++ ) {
1122 ni = (int)at[cur_at].neighbor[i];
1123 if ( at[ni].valence != 1 ||
1124 !get_endpoint_valence( at[ni].el_number ) ) {
1125 continue;
1126 }
1127 for ( j = i+1; j < at[cur_at].valence; j ++ ) {
1128 nj = (int)at[cur_at].neighbor[j];
1129 if ( at[nj].valence != 1 ||
1130 at[ni].el_number != at[nj].el_number ||
1131 !get_endpoint_valence( at[nj].el_number ) ) {
1132 continue;
1133 }
1134 /*
1135 * if (at[ni].num_H != at[nj].num_H) then the atoms (neighbors of at[cur_at]
1136 * are tautomeric endpoints and are indistinguishable => cur_at is not stereogenic
1137 * if (at[ni].num_H == at[nj].num_H) then the neighbors are indistinguishable
1138 * and cur_at will be found non-sterogenic later
1139 * get_endpoint_valence() check will not allow the neighbors to be carbons
1140 * Therefore the following "if" is not needed; we may just return 0.
1141 */
1142 if ( at[ni].num_H != at[nj].num_H && strcmp(at[ni].elname, "C" ) ) {
1143 return 0; /* found -X and -XH neighbors */
1144 }
1145 }
1146 }
1147 }
1148
1149 num_1s = num_mult = num_other = 0;
1150
1151 for ( i = 0; i < at[cur_at].valence; i ++ ) {
1152 bond_type = (at[cur_at].bond_type[i] & ~BOND_MARK_ALL);
1153 switch( bond_type ) {
1154 case BOND_SINGLE:
1155 num_1s ++;
1156 break;
1157 case BOND_DOUBLE:
1158 case BOND_ALTERN:
1159 case BOND_TAUTOM:
1160 case BOND_ALT12NS:
1161 num_mult ++;
1162 break;
1163 default:
1164 num_other ++;
1165 break;
1166 }
1167 }
1168
1169 if ( num_other ) {
1170 return 0;
1171 }
1172
1173 if ( (NumDbleBonds && NumDbleBonds > num_mult) ||
1174 (!NumDbleBonds && at[cur_at].valence != num_1s) ) {
1175 return 0;
1176 }
1177 return 1;
1178 }
1179 /********************************************************************************************/
bCanInpAtomBeAStereoCenter(inp_ATOM * at,int cur_at,int bPointedEdgeStereo)1180 int bCanInpAtomBeAStereoCenter( inp_ATOM *at, int cur_at, int bPointedEdgeStereo )
1181 {
1182
1183 /*************************************************************************************
1184 * current version
1185 *************************************************************************************
1186 * Use #define to split the stereocenter description table into parts
1187 * to make it easier to read
1188 *
1189 * --------- 4 single bonds stereocenters -------
1190 * 0 1 2 3 4 5
1191 *
1192 * | | | | | |
1193 * -C- -Si- -Ge- -Sn- >As[+] >B[-]
1194 * | | | | | |
1195 */
1196 #define SZELEM1 "C\000","Si", "Ge", "Sn", "As", "B\000",
1197 #define CCHARGE1 0, 0, 0, 0, 1, -1,
1198 #define CNUMBONDSANDH1 4, 4, 4, 4, 4, 4,
1199 #define CCHEMVALENCEH1 4, 4, 4, 4, 4, 4,
1200 #define CHAS3MEMBRING1 0, 0, 0, 0, 0, 0,
1201 #define CREQUIRDNEIGH1 0, 0, 0, 0, 3, 0,
1202 /*
1203 * --------------- S, Se stereocenters ----------
1204 * 6 7 8 9 10 11 12 13
1205 *
1206 * | | || | | ||
1207 * -S= =S= -S[+] >S[+] -Se= =Se= -Se[+] >Se[+]
1208 * | | | | | | | |
1209 */
1210 #define SZELEM2 "S\000","S\000","S\000","S\000","Se", "Se", "Se", "Se",
1211 #define CCHARGE2 0, 0, 1, 1, 0, 0, 1, 1,
1212 #define CNUMBONDSANDH2 3, 4, 3, 4, 3, 4, 3, 4,
1213 #define CCHEMVALENCEH2 4, 6, 3, 5, 4, 6, 3, 5,
1214 #define CHAS3MEMBRING2 0, 0, 0, 0, 0, 0, 0, 0,
1215 #define CREQUIRDNEIGH2 3, 3, 3, 3, 3, 3, 3, 3,
1216 /*
1217 * ------------------ N, P stereocenters -----------------
1218 * 14 15 16 17 18 19 20
1219 *
1220 * Phosphine Arsine
1221 * X---Y
1222 * | | \ / | | \ / \ /
1223 * =N- >N[+] N >P[+] =P- P As
1224 * | | | | | | |
1225 */
1226 #define SZELEM3 "N\000","N\000","N\000","P\000","P\000","P\000", "As",
1227 #define CCHARGE3 0, 1, 0, 1, 0, 0, 0,
1228 #define CNUMBONDSANDH3 4, 4, 3, 4, 4, 3, 3,
1229 #define CCHEMVALENCEH3 5, 4, 3, 4, 5, 3, 3,
1230 #define CHAS3MEMBRING3 0, 0, 1, 0, 0, 0, 0,
1231 #define CREQUIRDNEIGH3 3, 3, 1, 3, 3, 2, 2,
1232
1233 #define PHOSPHINE_STEREO 19 /* the number must match Phosphine number in the comments, see above */
1234 #define ARSINE_STEREO 20 /* the number must match Arsine number in the comments, see above */
1235
1236 static char szElem[][3]={ SZELEM1 SZELEM2 SZELEM3 };
1237 static S_CHAR cCharge[]={ CCHARGE1 CCHARGE2 CCHARGE3 };
1238 static S_CHAR cNumBondsAndH[]={ CNUMBONDSANDH1 CNUMBONDSANDH2 CNUMBONDSANDH3 };
1239 static S_CHAR cChemValenceH[]={ CCHEMVALENCEH1 CCHEMVALENCEH2 CCHEMVALENCEH3 };
1240 static S_CHAR cHas3MembRing[]={ CHAS3MEMBRING1 CHAS3MEMBRING2 CHAS3MEMBRING3 };
1241 static S_CHAR cRequirdNeigh[]={ CREQUIRDNEIGH1 CREQUIRDNEIGH2 CREQUIRDNEIGH3 };
1242
1243 static int n = sizeof(szElem)/sizeof(szElem[0]);
1244 /* reqired neighbor types (bitmap):
1245 0 => check bonds only
1246 1 => no terminal hydrogen atom neighbors
1247 2 => no terminal -X and -XH together (don't care the bond type, charge, radical)
1248 (X = tautomeric endpoint atom)
1249 Note: whenever cChemValenceH[] > cNumBondsAndH[]
1250 the tautomeric and/or alternating bonds
1251 are permitted
1252
1253 */
1254 int i, ret = 0;
1255 for ( i = 0; i < n; i++ ) {
1256 if ( !strcmp( at[cur_at].elname, szElem[i]) &&
1257 at[cur_at].charge == cCharge[i] &&
1258 (!at[cur_at].radical || at[cur_at].radical == 1) &&
1259 at[cur_at].valence +at[cur_at].num_H == cNumBondsAndH[i] &&
1260 at[cur_at].chem_bonds_valence+at[cur_at].num_H == cChemValenceH[i] &&
1261 (cHas3MembRing[i]? is_atom_in_3memb_ring( at, cur_at ) : 1) &&
1262 bInpAtomHasRequirdNeigh ( at, cur_at, cRequirdNeigh[i], cChemValenceH[i]-cNumBondsAndH[i]) ) {
1263 ret = cNumBondsAndH[i];
1264 break;
1265 }
1266 }
1267
1268 if ( i == PHOSPHINE_STEREO && !(bPointedEdgeStereo & PES_BIT_PHOSPHINE_STEREO) )
1269 ret = 0;
1270 if ( i == ARSINE_STEREO && !(bPointedEdgeStereo & PES_BIT_ARSINE_STEREO) )
1271 ret = 0;
1272 return ret;
1273 }
1274
1275 #else /* } NEW_STEREOCENTER_CHECK { */
1276
1277 /********************************************************************************************/
bCanAtomBeAStereoCenter(char * elname,S_CHAR charge,S_CHAR radical)1278 int bCanAtomBeAStereoCenter( char *elname, S_CHAR charge, S_CHAR radical )
1279 {
1280 static const char szElem[][3] = { "C\000", "Si", "Ge", "N\000", "P\000", "As", "B\000" };
1281 static const S_CHAR cCharge[] = { 0, 0, 0, 1, 1, 1, -1 };
1282 int i, ret = 0;
1283 for ( i = 0; i < sizeof(szElem)/sizeof(szElem[0]); i++ ) {
1284 if ( !strcmp( elname, szElem[i] ) && (charge == cCharge[i]) ) {
1285 ret = (!radical || radical == RADICAL_SINGLET);
1286 break;
1287 }
1288 }
1289 return ret;
1290 }
1291 #endif /* } NEW_STEREOCENTER_CHECK */
1292
1293 /****************************************************************/
1294 /* used for atoms adjacent to stereogenic bonds only */
bAtomHasValence3(char * elname,S_CHAR charge,S_CHAR radical)1295 int bAtomHasValence3( char *elname, S_CHAR charge, S_CHAR radical )
1296 {
1297 static const char szElem[][3] = { "N\000" };
1298 static const S_CHAR cCharge[] = { 0, };
1299 int i, ret = 0;
1300 for ( i = 0; i < (int)(sizeof(szElem)/sizeof(szElem[0])); i++ ) {
1301 if ( !strcmp( elname, szElem[i] ) && (charge == cCharge[i]) ) {
1302 ret = ( !radical || radical == RADICAL_SINGLET );
1303 break;
1304 }
1305 }
1306 return ret;
1307 }
1308
1309 /****************************************************************/
1310 /* used for atoms adjacent to stereogenic bonds only */
bCanAtomHaveAStereoBond(char * elname,S_CHAR charge,S_CHAR radical)1311 int bCanAtomHaveAStereoBond( char *elname, S_CHAR charge, S_CHAR radical )
1312 {
1313 static const char szElem[][3] = { "C\000", "Si", "Ge", "N\000", "N\000" };
1314 static const S_CHAR cCharge[] = { 0, 0, 0, 0, 1, };
1315 static const int n = sizeof(szElem)/sizeof(szElem[0]);
1316 int i, ret = 0;
1317 for ( i = 0; i < n; i++ ) {
1318 if ( !strcmp( elname, szElem[i] ) && (charge == cCharge[i]) ) {
1319 ret = (!radical || radical == RADICAL_SINGLET);
1320 break;
1321 }
1322 }
1323 return ret;
1324 }
1325 /****************************************************************/
1326 /* used for atoms adjacent to stereogenic bonds only */
bCanAtomBeMiddleAllene(char * elname,S_CHAR charge,S_CHAR radical)1327 int bCanAtomBeMiddleAllene( char *elname, S_CHAR charge, S_CHAR radical )
1328 {
1329 static const char szElem[][3] = { "C\000", "Si", "Ge", };
1330 static const S_CHAR cCharge[] = { 0, 0, 0, };
1331 static const int n = sizeof(szElem)/sizeof(szElem[0]);
1332 int i, ret = 0;
1333 for ( i = 0; i < n; i++ ) {
1334 if ( !strcmp( elname, szElem[i] ) && (charge == cCharge[i]) ) {
1335 ret = (!radical || radical == RADICAL_SINGLET);
1336 break;
1337 }
1338 }
1339 return ret;
1340 }
1341 /*****************************************************************/
bIsSuitableHeteroInpAtom(inp_ATOM * at)1342 int bIsSuitableHeteroInpAtom( inp_ATOM *at )
1343 {
1344 int val, num_H;
1345 if ( 0 == at->charge &&
1346 (!at->radical || RADICAL_SINGLET == at->radical) &&
1347 0 < (val=get_endpoint_valence( at->el_number ) )) {
1348 num_H = at->num_H;
1349 if ( val == at->chem_bonds_valence + num_H ) {
1350 switch( val ) {
1351 case 2: /* O */
1352 if ( !num_H && 1 == at->valence )
1353 return 0; /* =O */
1354 break; /* not found */
1355 case 3: /* N */
1356 if ( (1 == at->valence && 1 == num_H) ||
1357 (2 == at->valence && 0 == num_H) )
1358 return 1; /* =N- or =NH */
1359 break; /* not found */
1360 }
1361 }
1362 }
1363 return -1;
1364 }
1365 /****************************************************************/
bIsOxide(inp_ATOM * at,int cur_at)1366 int bIsOxide( inp_ATOM *at, int cur_at )
1367 {
1368 int i, bond_type;
1369 inp_ATOM *a = at + cur_at, *an;
1370 for ( i = 0; i < a->valence; i ++ ) {
1371 bond_type = (a->bond_type[i] &= ~BOND_MARK_ALL);
1372 if ( bond_type == BOND_DOUBLE ) {
1373 an = at + (int)a->neighbor[i];
1374 if ( 1 == an->valence &&
1375 !an->charge && !an->num_H && !an->radical &&
1376 2 == get_endpoint_valence( an->el_number ) ) {
1377 return 1;
1378 }
1379 } else
1380 if ( bond_type == BOND_TAUTOM || bond_type == BOND_ALT12NS ) {
1381 an = at + (int)a->neighbor[i];
1382 if ( 1 == an->valence &&
1383 2 == get_endpoint_valence( an->el_number ) ) {
1384 return 1;
1385 }
1386 }
1387 }
1388 return 0;
1389 }
1390 /****************************************************************/
1391 /* used for atoms adjacent to stereogenic bonds only */
bCanAtomBeTerminalAllene(char * elname,S_CHAR charge,S_CHAR radical)1392 int bCanAtomBeTerminalAllene( char *elname, S_CHAR charge, S_CHAR radical )
1393 {
1394 static const char szElem[][3] = { "C\000", "Si", "Ge", };
1395 static const S_CHAR cCharge[] = { 0, 0, 0, };
1396 static const int n = sizeof(szElem)/sizeof(szElem[0]);
1397 int i, ret = 0;
1398 for ( i = 0; i < n; i++ ) {
1399 if ( !strcmp( elname, szElem[i] ) && (charge == cCharge[i]) ) {
1400 ret = (!radical || radical == RADICAL_SINGLET);
1401 break;
1402 }
1403 }
1404 return ret;
1405 }
1406 /************************************************************************/
GetHalfStereobond0DParity(inp_ATOM * at,int cur_at,AT_NUMB nSbNeighOrigAtNumb[],int nNumExplictAttachments,int bond_parity,int nFlag)1407 int GetHalfStereobond0DParity( inp_ATOM *at, int cur_at, AT_NUMB nSbNeighOrigAtNumb[],
1408 int nNumExplictAttachments, int bond_parity, int nFlag )
1409 {
1410 int m, last_parity, cur_parity;
1411 int i, icur2nxt, icur2neigh, cur_order_parity, nxt_at;
1412 AT_NUMB nNextSbAtOrigNumb;
1413 /* find atom parities for all valid streobonds incident to at[cur_at] */
1414 for ( m = 0, last_parity = 0; m < MAX_NUM_STEREO_BONDS && at[cur_at].sb_parity[m]; m ++ ) {
1415 icur2nxt = icur2neigh = -1; /* ordering number of neighbors in nSbNeighOrigAtNumb[] */
1416 cur_parity = 0; /* parity for mth stereobond incident to the cur_at */
1417 if ( 0 <= at[cur_at].sb_ord[m] && at[cur_at].sb_ord[m] < at[cur_at].valence &&
1418 0 <= (nxt_at = at[cur_at].neighbor[(int)at[cur_at].sb_ord[m]]) &&
1419 at[nxt_at].valence <= MAX_NUM_STEREO_BONDS && /* make sure it is a valid stereobond */
1420 (nNextSbAtOrigNumb = at[nxt_at].orig_at_number) ) {
1421 /* since at[cur_at].sn_ord[m] = -1 for explicit H use at[cur_at].sn_orig_at_num[m] */
1422 for ( i = 0; i < nNumExplictAttachments; i ++ ) {
1423 if ( at[cur_at].sn_orig_at_num[m] == nSbNeighOrigAtNumb[i] ) {
1424 icur2neigh = i; /* neighbor */
1425 } else
1426 if ( nNextSbAtOrigNumb == nSbNeighOrigAtNumb[i] ) {
1427 icur2nxt = i; /* atom connected by a stereobond */
1428 }
1429 }
1430 if ( icur2neigh >= 0 && icur2nxt >= 0 ) {
1431 if ( ATOM_PARITY_WELL_DEF(at[cur_at].sb_parity[m]) ) {
1432 /* parity of at[cur_atom] neighbor permutation to reach this order: { next_atom, neigh_atom, ...} */
1433 cur_order_parity = (icur2nxt + icur2neigh + (icur2nxt > icur2neigh) - 1) % 2;
1434 cur_parity = 2 - (cur_order_parity + at[cur_at].sb_parity[m]) % 2;
1435 } else {
1436 /* unknowm/undef parities do not depend on the neighbor order */
1437 cur_parity = at[cur_at].sb_parity[m];
1438 }
1439 }
1440 } else {
1441 continue;
1442 }
1443 /* use a well-known parity if available; if not then use preferably the unknown */
1444 if ( !last_parity ) {
1445 last_parity = cur_parity;
1446 } else
1447 if ( last_parity != cur_parity && cur_parity ) {
1448 if ( ATOM_PARITY_WELL_DEF(last_parity) ) {
1449 if ( ATOM_PARITY_WELL_DEF(cur_parity) ) {
1450 last_parity = 0; /* error: all well-defined parities should be same */
1451 break;
1452 }
1453 } else
1454 if ( ATOM_PARITY_WELL_DEF(cur_parity) ) {
1455 /* replace unknown/undefined parity with well-known */
1456 last_parity = cur_parity;
1457 } else {
1458 /* select min unknown/undefined parity (out of AB_PARITY_UNKN and AB_PARITY_UNDF) */
1459 last_parity = inchi_min(cur_parity, last_parity);
1460 }
1461 }
1462 }
1463 if ( last_parity ) {
1464 bond_parity = last_parity;
1465 at[cur_at].bUsed0DParity |= nFlag; /* set flag: used stereobond 0D parity */
1466 }
1467 return bond_parity;
1468 }
1469 /*******************************************************************************************/
FixSb0DParities(inp_ATOM * at,int chain_length,int at_1,int i_next_at_1,S_CHAR z_dir1[],int at_2,int i_next_at_2,S_CHAR z_dir2[],int * pparity1,int * pparity2)1470 int FixSb0DParities( inp_ATOM *at, /* inp_ATOM *at_removed_H, int num_removed_H,*/ int chain_length,
1471 int at_1, int i_next_at_1, S_CHAR z_dir1[],
1472 int at_2, int i_next_at_2, S_CHAR z_dir2[],
1473 int *pparity1, int *pparity2 )
1474 {
1475 int k, parity1, parity2, abs_parity1, abs_parity2;
1476 int j1, j2, parity_sign;
1477 /*
1478 AT_NUMB nSbNeighOrigAtNumb1[MAX_NUM_STEREO_BOND_NEIGH], nSbNeighOrigAtNumb2[MAX_NUM_STEREO_BOND_NEIGH];
1479 int nNumExplictAttachments1, nNumExplictAttachments2;
1480 */
1481 parity1 = parity2 = AB_PARITY_NONE;
1482 j1 = j2 = -1;
1483 parity_sign = ( *pparity1 < 0 || *pparity2 < 0 )? -1 : 1;
1484
1485 abs_parity1 = abs(*pparity1);
1486 abs_parity2 = abs(*pparity2);
1487
1488 for ( k = 0; k < MAX_NUM_STEREO_BONDS && at[at_1].sb_parity[k]; k ++ ) {
1489 if ( at[at_1].sb_ord[k] == i_next_at_1 ) {
1490 parity1 = at[at_1].sb_parity[k];
1491 j1 = k;
1492 }
1493 }
1494 for ( k = 0; k < MAX_NUM_STEREO_BONDS && at[at_2].sb_parity[k]; k ++ ) {
1495 if ( at[at_2].sb_ord[k] == i_next_at_2 ) {
1496 parity2 = at[at_2].sb_parity[k];
1497 j2 = k;
1498 }
1499 }
1500 switch( (j1 >= 0) + 2*(j2 >= 0) ) {
1501 case 0:
1502 /* the bond has no 0D parity */
1503 *pparity1 = *pparity2 = parity_sign * AB_PARITY_UNDF;
1504 return 0;
1505 case 1:
1506 case 2:
1507 /* 0D parity data error */
1508 *pparity1 = *pparity2 = AB_PARITY_NONE;
1509 return -1;
1510 case 3:
1511 /* the bond has 0D parity */
1512 switch ( !(ATOM_PARITY_WELL_DEF( abs_parity1 ) && ATOM_PARITY_WELL_DEF( parity1 )) +
1513 2 * !(ATOM_PARITY_WELL_DEF( abs_parity2 ) && ATOM_PARITY_WELL_DEF( parity2 )) ) {
1514 case 0:
1515 /* both parities are well-defined; continue */
1516 break;
1517 case 1:
1518 /* 0D parity not well-defined for at_1 */
1519 *pparity1 = parity_sign * (ATOM_PARITY_WELL_DEF( parity1 )? abs_parity1 :
1520 ATOM_PARITY_WELL_DEF( abs_parity1 )? parity1 :
1521 inchi_min(abs_parity1, parity1));
1522 *pparity2 = parity_sign * abs_parity2;
1523 return -1;
1524 case 2:
1525 /* 0D parity not well-defined for at_2 */
1526 *pparity1 = parity_sign * abs_parity1;
1527 *pparity2 = parity_sign * (ATOM_PARITY_WELL_DEF( parity2 )? abs_parity2 :
1528 ATOM_PARITY_WELL_DEF( abs_parity2 )? parity2 :
1529 inchi_min(abs_parity2, parity2));
1530 return -1;
1531 case 3:
1532 abs_parity1 = (ATOM_PARITY_WELL_DEF( parity1 )? abs_parity1 :
1533 ATOM_PARITY_WELL_DEF( abs_parity1 )? parity1 :
1534 inchi_min(abs_parity1, parity1));
1535 abs_parity2 = (ATOM_PARITY_WELL_DEF( parity2 )? abs_parity2 :
1536 ATOM_PARITY_WELL_DEF( abs_parity2 )? parity2 :
1537 inchi_min(abs_parity2, parity2));
1538 *pparity1 = *pparity2 = parity_sign * inchi_min(abs_parity1, abs_parity2);
1539 /*return (parity1 == parity2)? 0 : -1;*/
1540 return -1;
1541 }
1542 break;
1543 }
1544 /* we are here if both end-atoms of the bond have well-defined 0D parities */
1545 /*
1546 nNumExplictAttachments1 = GetSbNeighOrigAtNumb( at, at_1, at_removed_H, num_removed_H, nSbNeighOrigAtNumb1 );
1547 nNumExplictAttachments2 = GetSbNeighOrigAtNumb( at, at_2, at_removed_H, num_removed_H, nSbNeighOrigAtNumb2 );
1548 parity1 = GetHalfStereobond0DParity( at, at_1, nSbNeighOrigAtNumb1, nNumExplictAttachments1, *pparity1, 0 );
1549 parity2 = GetHalfStereobond0DParity( at, at_2, nSbNeighOrigAtNumb2, nNumExplictAttachments2, *pparity2, 0 );
1550 */
1551 *pparity1 = parity_sign * abs_parity1;
1552 *pparity2 = parity_sign * abs_parity2;
1553
1554 if ( chain_length % 2 ) {
1555 /* allene; chain_length = (number of double bonds) - 1 */
1556 /*
1557 int zer1 = ( !z_dir1[0] && !z_dir1[1] && !z_dir1[2] );
1558 int zer2 = ( !z_dir2[0] && !z_dir2[1] && !z_dir2[2] );
1559 */
1560 int bWrong_z_dir1 = (0 != (at[at_1].bUsed0DParity & FlagSB_0D));
1561 int bWrong_z_dir2 = (0 != (at[at_2].bUsed0DParity & FlagSB_0D));
1562
1563 if ( bWrong_z_dir1 && bWrong_z_dir2 ) {
1564 goto set_default;
1565 } else
1566 if ( bWrong_z_dir1 || bWrong_z_dir2 ) {
1567 double r12[3], zi1[3], zi2[3], abs_r12, abs_zi2;
1568 int at_i1, at_i2, j;
1569 S_CHAR z_dir[3];
1570 r12[0] = at[at_2].x - at[at_1].x;
1571 r12[1] = at[at_2].y - at[at_1].y;
1572 r12[2] = at[at_2].z - at[at_1].z;
1573 abs_r12 = len3( r12 );
1574 if ( abs_r12 < MIN_BOND_LEN ) {
1575 goto set_default;
1576 }
1577 /* make r12[] point to the atom with 'good' z_dir[] */
1578 if ( bWrong_z_dir1 ) {
1579 at_i1 = at_2; /* has good z_dir2[] */
1580 at_i2 = at_1; /* has bad z_dir1[] */
1581 zi1[0] = z_dir2[0];
1582 zi1[1] = z_dir2[1];
1583 zi1[2] = z_dir2[2];
1584 mult3( r12, 1.0/abs_r12, r12 ); /* make length = 1 */
1585 } else {
1586 at_i1 = at_1; /* has good z_dir1[] */
1587 at_i2 = at_2; /* has bad z_dir2[] */
1588 zi1[0] = z_dir1[0];
1589 zi1[1] = z_dir1[1];
1590 zi1[2] = z_dir1[2];
1591 mult3( r12, -1.0/abs_r12, r12 ); /* make length = 1 */
1592 }
1593 cross_prod3( r12, zi1, zi2 );
1594 abs_zi2 = len3( zi2 );
1595 mult3( zi2, 100.0/abs_zi2, zi2 ); /* make length = 100 */
1596 for ( j = 0; j < 3; j ++ ) {
1597 z_dir[j] = (S_CHAR) (zi2[j]>= 0.0? floor(0.5 + zi2[j]) :
1598 -floor(0.5 - zi2[j])); /* abs(z_dir) = 100 */
1599 }
1600 if ( bWrong_z_dir1 ) {
1601 memcpy( z_dir1, z_dir, sizeof(z_dir) );
1602 } else {
1603 memcpy( z_dir2, z_dir, sizeof(z_dir) );
1604 }
1605 }
1606 return 0;
1607
1608 set_default:
1609 /* z_dir1[] = x-direction; z_dir2[] = z-direction; r12[] = y-direction */
1610 z_dir1[0] = 100;
1611 z_dir1[1] = z_dir1[2] = 0;
1612 z_dir2[0] = z_dir2[1] = 0;
1613 z_dir2[2] = 100;
1614 }
1615 return 0;
1616 }
1617 /**********************************************************/
1618 /* without this InChI fails on reconstructed CID=450438 */
1619 /* (isotopic, Unknown SB adjacent to SB with known parity) */
1620 /**********************************************************/
FixUnkn0DStereoBonds(inp_ATOM * at,int num_at)1621 int FixUnkn0DStereoBonds(inp_ATOM *at, int num_at)
1622 {
1623 int i, m, num=0;
1624
1625 /* add usual Unknown stereobond descriptors to each Unknown bond */
1626 for( i = 0; i < num_at; i ++ ) {
1627 for ( m = 0; m < MAX_NUM_STEREO_BONDS && at[i].sb_parity[m]; m ++ ) {
1628 if ( AB_PARITY_UNKN == at[i].sb_parity[m] ) {
1629 at[i].bond_stereo[ (int)at[i].sb_ord[m] ] = STEREO_DBLE_EITHER;
1630 num ++;
1631 }
1632 }
1633 }
1634 #ifdef NEVER
1635 if ( num ) {
1636 int j;
1637 /* how to remove Unknown stereo bond parities */
1638 for( i = 0; i < num_at; i ++ ) {
1639 for ( m = 0; m < MAX_NUM_STEREO_BONDS && at[i].sb_parity[m]; m ++ ) {
1640 if ( AB_PARITY_UNKN == at[i].sb_parity[m] ) {
1641 for ( j = m+1; j < MAX_NUM_STEREO_BONDS; j ++ ) {
1642 at[i].sb_parity[j-1] = at[i].sb_parity[j];
1643 at[i].sb_ord[j-1] = at[i].sb_ord[j];
1644 at[i].sn_ord[j-1] = at[i].sn_ord[j];
1645 at[i].sn_orig_at_num[j-1] = at[i].sn_orig_at_num[j];
1646 }
1647 at[i].sb_parity[j-1] = 0;
1648 at[i].sb_ord[j-1] = 0;
1649 at[i].sn_ord[j-1] = 0;
1650 at[i].sn_orig_at_num[j-1] = 0;
1651 }
1652 }
1653 }
1654 }
1655 #endif
1656 return num;
1657 }
1658 /*======================================================================================================
1659
1660 half_stereo_bond_parity() General Description:
1661
1662 A) find projections of 3 bonds on a reasonable plane defined
1663 by a vector z_dir perpendicular to the plane
1664 B) calculate parity
1665
1666 half_stereo_bond_parity() Detailed Description:
1667
1668 1) Find at_coord[] = vectors from the central atoms to its neighbors
1669 2) If only 2 neighbors are present, then create a reasonable 3rd neighbor
1670 (an implicit H or a fictitious atom in case of =NX) coordinates
1671 3) Normalize at_coord[] to unit length
1672 4) Find unit vector pnt[2] perpendicular to the plane containing
1673 at_coord[] arrow ends.
1674 Even though it is not necessary, make z-coordinate of pnt[2] positive.
1675 ** pnt[2] has the new z-axis direction **
1676 5) Let pnt[0] = perpendicular to pnt[2] component of at_coord[0];
1677 Normalize pnt[0] to unit length.
1678 ** pnt[0] has the new x-axis direction **
1679 6) Let pnt[1] = pnt[2] x pnt[0] (cross-product);
1680 ** pnt[1] has the new y-axis direction **
1681 7) Find at_coord[] in the new xyz-basis and normalize their xy-projections
1682 to a unit length
1683 8) In the new xy-plane find (counterclockwise) angles:
1684 tmp1 = (from at_coord[0] to at_coord[1])
1685 tmp2 = (from at_coord[0] to at_coord[2])
1686 9) Calculate the parity: if tmp1 < tmp2 then 1 (odd) else 2 (even)
1687 (even: looking from the arrow end of the new z-axis, 0, 1, and 2 neighbors
1688 are in clockwise order)
1689 10) Calculate z_dir = 100*pnt[2].
1690
1691 Note1. If z_dir vectors of atoms located at the opposite ends of a double bond have approximately
1692 opposite directions (that is, their dot-product is negative) then the parity of the
1693 stereogenic bond calculated from half-bond-parities should be inverted
1694
1695 Note2. In case of a tetrahedral cumulene a triple product (z_dir1, (1->2), z_dir2) is used instead
1696 of the dot-product. (1->2) is a vector from the atom#1 to the atom #2. This triple product
1697 is invariant with respect to the atom numbering because it does not change upon (1,2)
1698 permutation.
1699
1700 Stereo ambiguity in case of 2 neighbors:
1701 ----------------------------------------
1702 Undefined: single-double bond angle > pi - arcsin(0.03) = 178.28164199834454285275613218975 degrees
1703 Ambiguous: single-double bond angle > 175 degrees = pi - 0.087156 Rad
1704
1705 Return values
1706 (cases: I=only in case of isotopic H atoms the neighbors are different,
1707 N=in case of non-isotopic H atoms the neighbors are different)
1708
1709 -4 = AB_PARITY_UNDF => atom is adjacent to a stereogenic bond, but the geometry is undefined, I
1710 -3 = AB_PARITY_UNKN => atom is adjacent to a stereogenic bond, but the geometry is not known to the iuser, I
1711 -2 =-AB_PARITY_EVEN => parity of an atom adjacent to a stereogenic bond, I
1712 -1 =-AB_PARITY_ODD => parity of an atom adjacent to a stereogenic bond, I
1713 0 = AB_PARITY_NONE => the atom is not adjacent to a stereogenic bond
1714 1 = AB_PARITY_ODD => parity of an atom adjacent to a stereogenic bond, N&I
1715 2 = AB_PARITY_EVEN => parity of an atom adjacent to a stereogenic bond, N&I
1716 3 = AB_PARITY_UNKN => atom is adjacent to a stereogenic bond, but the geometry is not known to the iuser, N&I
1717 4 = AB_PARITY_UNDF => atom is adjacent to a stereogenic bond, but the geometry is undefined, N&I
1718 5 = AB_PARITY_IISO => atom constitutionally equivalent to this atom may be adjacent to a stereogenic bond, I
1719
1720
1721 =====================================================================================================*/
1722
half_stereo_bond_parity(inp_ATOM * at,int cur_at,inp_ATOM * at_removed_H,int num_removed_H,S_CHAR * z_dir,int bPointedEdgeStereo,int vABParityUnknown)1723 int half_stereo_bond_parity( inp_ATOM *at, int cur_at, inp_ATOM *at_removed_H,
1724 int num_removed_H, S_CHAR *z_dir,
1725 int bPointedEdgeStereo, int vABParityUnknown )
1726 {
1727 double at_coord[MAX_NUM_STEREO_BOND_NEIGH][3], c, s, tmp[3], tmp1, tmp2, min_tmp, max_tmp, z;
1728 double temp[3], pnt[3][3];
1729 int j, k, p0, p1, p2, next, bValence3=0, num_z, nType, num_either_single, num_either_double;
1730 int nNumExplictAttachments;
1731 int bond_parity = AB_PARITY_UNDF;
1732 int num_H=0, num_iH, num_eH=0, num_nH=0 /* = num_iso_H[0] */;
1733 int num_iso_H[NUM_H_ISOTOPES+1];
1734 int index_H[5]; /* cannot have more than 4 elements: 1 H, 1 1H, 1 D, 1 T atom(s) */
1735 /* const double one_pi = 2.0*atan2(1.0 , 0.0 ); */
1736 const double one_pi = 3.14159265358979323846; /* M_PI */
1737 const double two_pi = 2.0*one_pi;
1738 int bIgnoreIsotopicH = (0 != (at[cur_at].cFlags & AT_FLAG_ISO_H_POINT));
1739 AT_NUMB nSbNeighOrigAtNumb[MAX_NUM_STEREO_BOND_NEIGH];
1740
1741
1742 if ( z_dir && !z_dir[0] && !z_dir[1] && !z_dir[2] ) {
1743 z_dir[2]=100;
1744 }
1745
1746 num_H = at[cur_at].num_H;
1747 if ( num_H > NUM_H_ISOTOPES )
1748 return 0; /* at least 2 H atoms are isotopically identical */
1749
1750 if ( MAX_NUM_STEREO_BOND_NEIGH < at[cur_at].valence + num_H ||
1751 MIN_NUM_STEREO_BOND_NEIGH > at[cur_at].valence + num_H )
1752 return 0;
1753
1754 if ( !bCanAtomHaveAStereoBond( at[cur_at].elname, at[cur_at].charge, at[cur_at].radical ) )
1755 return 0;
1756 if ( !bIgnoreIsotopicH ) {
1757 for ( j = 0, num_nH = num_H; j < NUM_H_ISOTOPES; j ++ ) {
1758 if ( (k = (int)at[cur_at].num_iso_H[j]) > 1 ) {
1759 return AB_PARITY_IISO; /* two or more identical isotopic H atoms */
1760 }
1761 num_nH -= k;
1762 }
1763 }
1764 /* at this point num_nH = number of non-isotopic H atoms */
1765 if ( num_nH > 1 )
1766 return AB_PARITY_IISO; /* two or more identical non-isotopic H atoms */
1767 if ( num_nH < 0 )
1768 return CT_ISO_H_ERR; /* program error */ /* <BRKPT> */
1769
1770 /********************************************************************
1771 * Note. At this point all (implicit and explicit) isotopic
1772 * terminal H neighbors are either different or not present.
1773 ********************************************************************/
1774
1775 /* locate explicit hydrogen atoms */
1776 /* (at_removed_H are sorted in ascending isotopic H mass order, non-isotopic first) */
1777 memset( num_iso_H, 0, sizeof(num_iso_H) );
1778 if ( at_removed_H && num_removed_H > 0 ) {
1779 for ( j = 0; j < num_removed_H; j ++ ) {
1780 if ( at_removed_H[j].neighbor[0] == cur_at ) {
1781 k = bIgnoreIsotopicH? 0 : at_removed_H[j].iso_atw_diff;
1782 if ( 0 <= k && k <= NUM_H_ISOTOPES ) {
1783 if ( ++num_iso_H[k] > 1 ) /* num_iso_H[0] = number of non-isotopic H atoms */
1784 return CT_ISO_H_ERR; /* program error in counting hydrogens */ /* <BRKPT> */
1785 index_H[num_eH++] = j;
1786 } else {
1787 return CT_ISO_H_ERR; /* program error */ /* <BRKPT> */
1788 }
1789 }
1790 }
1791 num_iH = num_H - num_eH; /* number of implicit non-isotopic and isotopic H atoms */
1792 if ( num_iH > 1 ) {
1793 /* more than one implicit H: cannot reconstruct the geometry */
1794 bond_parity = -AB_PARITY_UNDF;
1795 goto exit_function;
1796 }
1797 } else {
1798 num_iH = num_H;
1799 }
1800 /* at this point num_iH = number of implicit non-isotopic and isotopic H atoms */
1801 if ( at[cur_at].valence + num_eH < MIN_NUM_STEREO_BOND_NEIGH ) {
1802 /* =NH or =CHD when no explicit H is present */
1803 return num_H == 1? AB_PARITY_UNDF : -AB_PARITY_UNDF;
1804 }
1805
1806 bValence3 = bAtomHasValence3( at[cur_at].elname, at[cur_at].charge, at[cur_at].radical );
1807 /*
1808 * Can one explicit hydrogen be added to make asymmetric configuration?
1809 * For now we can add 1 H atom in case of an appropriate geometry if:
1810 * (a) one non-isotopic H (even if explicit isotopic H atoms are present), or
1811 * (b) one isotopic or non-isotopic H if NO explicit isotopic or non-isotopic H atom is present
1812 * This makes sense only in case chem. valence = 4. In case of chem. valence = 3, do not check.
1813 */
1814 if ( at[cur_at].valence + num_eH == MIN_NUM_STEREO_BOND_NEIGH && !bValence3 &&
1815 !(/*(a)*/ (1 == num_nH && !num_iso_H[0]) ||
1816 /*(b)*/ (1 == num_H && !num_eH))
1817 ) {
1818 goto exit_function;
1819 /* return num_H == 1? AB_PARITY_UNDF : -AB_PARITY_UNDF; */
1820 }
1821
1822 /* store neighbors coordinates */
1823 num_z = num_either_single = num_either_double = 0;
1824 for ( k = nNumExplictAttachments = 0; k < 2; k ++ ) {
1825 switch( k ) {
1826 case 0:
1827 for ( j = 0; j < num_eH; j ++, nNumExplictAttachments ++ ) {
1828 next = index_H[j];
1829 at_coord[nNumExplictAttachments][0] = at_removed_H[next].x - at[cur_at].x;
1830 at_coord[nNumExplictAttachments][1] = at_removed_H[next].y - at[cur_at].y;
1831 nSbNeighOrigAtNumb[nNumExplictAttachments] = at_removed_H[next].orig_at_number;
1832 /* use the fact that (at_removed_H - at) = (number of atoms except removed explicit H) */
1833 z = -get_z_coord( at, (at_removed_H-at)+next, 0 /*neighbor #*/, &nType, -(bPointedEdgeStereo & PES_BIT_POINT_EDGE_STEREO) );
1834 switch ( nType ) {
1835 case ZTYPE_EITHER:
1836 num_either_single ++; /* bond in "Either" direction. */
1837 break;
1838 case ZTYPE_UP:
1839 case ZTYPE_DOWN:
1840 nType = -nType; /* at_removed_H[] contains bonds TO the center, not from */
1841 z = len2( at_coord[nNumExplictAttachments] );
1842 /*
1843 z = sqrt( at_coord[nNumExplictAttachments][0]*at_coord[nNumExplictAttachments][0]
1844 + at_coord[nNumExplictAttachments][1]*at_coord[nNumExplictAttachments][1] );
1845 */
1846 if ( nType == ZTYPE_DOWN )
1847 z = -z;
1848 /* no break; here */
1849 case ZTYPE_3D:
1850 num_z ++;
1851 }
1852 at_coord[nNumExplictAttachments][2] = z;
1853 }
1854 break;
1855 case 1:
1856 for ( j = 0; j < at[cur_at].valence; j ++, nNumExplictAttachments ++ ) {
1857 next = at[cur_at].neighbor[j];
1858 at_coord[nNumExplictAttachments][0] = at[next].x - at[cur_at].x;
1859 at_coord[nNumExplictAttachments][1] = at[next].y - at[cur_at].y;
1860 nSbNeighOrigAtNumb[nNumExplictAttachments] = at[next].orig_at_number;
1861
1862 z = get_z_coord( at, cur_at, j /*neighbor #*/, &nType, (bPointedEdgeStereo & PES_BIT_POINT_EDGE_STEREO) );
1863 switch ( nType ) {
1864 case ZTYPE_EITHER:
1865 num_either_single ++; /* bond in "Either" direction. */
1866 break;
1867 case ZTYPE_UP:
1868 case ZTYPE_DOWN:
1869 z = len2( at_coord[nNumExplictAttachments] );
1870 /*
1871 z = sqrt( at_coord[nNumExplictAttachments][0]*at_coord[nNumExplictAttachments][0]
1872 + at_coord[nNumExplictAttachments][1]*at_coord[nNumExplictAttachments][1] );
1873 */
1874 if ( nType == ZTYPE_DOWN )
1875 z = -z;
1876 /* no break; here */
1877 case ZTYPE_3D:
1878 num_z ++;
1879 }
1880 at_coord[nNumExplictAttachments][2] = z;
1881 }
1882 break;
1883 }
1884 }
1885
1886 if ( num_either_single ) {
1887 bond_parity = vABParityUnknown /*AB_PARITY_UNKN*/; /* single bond is 'unknown' */
1888 goto exit_function;
1889 }
1890
1891 /* nNumExplictAttachments is a total number of attachments, including removed explicit terminal hydrogens */
1892 if ( nNumExplictAttachments == 2 ) {
1893 /* create coordinates of the implicit hydrogen (or a fictitious atom in case of ==N-X ), */
1894 /* coord[2][], attached to the cur_at. */
1895 for ( j = 0; j < 3; j ++ ) {
1896 at_coord[2][j] = - ( at_coord[0][j] + at_coord[1][j] );
1897 }
1898 nSbNeighOrigAtNumb[nNumExplictAttachments] = 0; /* implicit H or lone pair */
1899 }
1900 for ( j = 0; j < 3; j ++ ) {
1901 tmp[j] = len3( at_coord[j] );
1902 }
1903 min_tmp = inchi_min( tmp[0], inchi_min(tmp[1], tmp[2]) );
1904 max_tmp = inchi_max( tmp[0], inchi_max(tmp[1], tmp[2]) );
1905 if ( min_tmp < MIN_BOND_LEN || min_tmp < MIN_SINE*max_tmp ) {
1906 /* all bonds or some of bonds are too short */
1907 if ( at[cur_at].sb_parity[0] ) {
1908 /* use bond psrity; the reconciliation in ReconcileAllCmlBondParities()
1909 * has made all ways to calculate parity produce same result
1910 */
1911 bond_parity = GetHalfStereobond0DParity( at, cur_at, nSbNeighOrigAtNumb,
1912 nNumExplictAttachments, bond_parity, FlagSB_0D );
1913 }
1914
1915 goto exit_function;
1916 }
1917 /* normalize lengths to 1 */
1918 for ( j = 0; j < 3; j ++ ) {
1919 mult3( at_coord[j], 1.0/tmp[j], at_coord[j] );
1920 }
1921
1922 /* find projections of at_coord vector differences on the plane containing their arrowhead ends */
1923 for ( j = 0; j < 3; j ++ ) {
1924 /* pnt[0..2] = {0-1, 1-2, 2-0} */
1925 tmp[j] = len3(diff3( at_coord[j], at_coord[(j+1)%3], pnt[j] ));
1926 if ( tmp[j] < MIN_SINE ) {
1927 goto exit_function; /* angle #i-cur_at-#j is too small */
1928 }
1929 mult3( pnt[j], 1.0/tmp[j], pnt[j] ); /* 2003-10-06 */
1930 }
1931 /* find pnt[p2], a vector perpendicular to the plane, and its length tmp[p2] */
1932 /* replace previous pnt[p2], tmp[p2] with new values; the old values do not have any additional */
1933 /* information because pnt[p0]+pnt[p1]+pnt[p2]=0 */
1934 /* 10-6-2003: a cross-product of one pair pnt[j], pnt[(j+1)%3] can be very small. Find the larges one */
1935 tmp1 = len3( cross_prod3( pnt[0], pnt[1], temp ) );
1936 for (j = 1, k = 0; j < 3; j ++ ) {
1937 tmp2 = len3( cross_prod3( pnt[j], pnt[(j+1)%3], temp ) );
1938 if ( tmp2 > tmp1 ) {
1939 tmp1 = tmp2;
1940 k = j;
1941 }
1942 }
1943 /* previously p0=0, p1=1, p2=2 */
1944 p0 = k;
1945 p1 = (k+1)%3;
1946 p2 = (k+2)%3;
1947 tmp[p2] = len3( cross_prod3( pnt[p0], pnt[p1], pnt[p2] ) );
1948 if ( tmp[p2] < MIN_SINE*tmp[p0]*tmp[p1] ) {
1949 goto exit_function; /* pnt[p0] is almost colinear to pnt[p1] */
1950 }
1951 /* new basis: pnt[p0], pnt[p1], pnt[p2]; set z-coord sign and make abs(pnt[p2]) = 1 */
1952 mult3( pnt[p2], (pnt[p2][2]>0.0? 1.0:-1.0)/tmp[p2], pnt[p2] ); /* unit vector in the new z-axis direction */
1953
1954 min_tmp = dot_prod3( at_coord[0], pnt[p2] ); /* non-planarity measure (sine): hight of at_coord[] pyramid */
1955 mult3( pnt[p2], min_tmp, pnt[p0] ); /* vector height of the pyramid, ideally 0 */
1956 /* find new pnt[p0] = projection of at_coord[p0] on plane orthogonal to pnt[p2] */
1957 tmp[p0] = len3(diff3( at_coord[0], pnt[p0], pnt[p0] ));
1958 mult3( pnt[p0], 1.0/tmp[p0], pnt[p0] ); /* new x axis basis vector */
1959 cross_prod3( pnt[p2], pnt[p0], pnt[p1] ); /* new y axis basis vector */
1960 /* find at_coord in the new basis of {pnt[p0], pnt[p1], pnt[p2]} */
1961 for ( j = 0; j < 3; j ++ ) {
1962 copy3( at_coord[j], temp );
1963 for ( k = 0; k < 3; k ++ ) {
1964 at_coord[j][k] = dot_prod3( temp, pnt[(k+p0)%3] );
1965 }
1966 /* new xy plane projection length */
1967 tmp[j] = sqrt(at_coord[j][0]*at_coord[j][0] + at_coord[j][1]*at_coord[j][1]);
1968 /* make new xy plane projection length = 1 */
1969 mult3( at_coord[j], 1.0/tmp[j], at_coord[j] );
1970 }
1971
1972 s = fabs( at_coord[1][0]*at_coord[2][1] - at_coord[1][1]*at_coord[2][0] ); /* 1-2 sine */
1973 c = at_coord[1][0]*at_coord[2][0] + at_coord[1][1]*at_coord[2][1]; /* 1-2 cosine */
1974 if ( s < MIN_SINE && c > 0.5 ) {
1975 goto exit_function; /* bonds to neigh. 1 and 2 have almost same direction; relative angles are undefined */
1976 }
1977 c = at_coord[0][0]; /* cosine of the angle between new Ox axis and a bond to the neighbor 0. Should be 1 */
1978 s = at_coord[0][1]; /* sine. Should be 0 */
1979 /* turn vectors so that vector #1 (at_coord[0]) becomes {1, 0} */
1980 for ( j = 0; j < MAX_NUM_STEREO_BOND_NEIGH; j ++ ) {
1981 tmp1 = c*at_coord[j][0] + s*at_coord[j][1];
1982 tmp2 = -s*at_coord[j][0] + c*at_coord[j][1];
1983 at_coord[j][0] = tmp1;
1984 at_coord[j][1] = tmp2;
1985 }
1986 /* counterclockwise angles from the direction to neigh 0 to to directions to neighbors 1 and 2: */
1987 tmp1 = atan2( at_coord[1][1], at_coord[1][0] ); /* range -pi and +pi */
1988 tmp2 = atan2( at_coord[2][1], at_coord[2][0] );
1989 if ( tmp1 < 0.0 )
1990 tmp1 += two_pi; /* range 0 to 2*pi */
1991 if ( tmp2 < 0.0 )
1992 tmp2 += two_pi;
1993 /*-----------------------------------
1994 Example
1995 1 \ case tmp1 < tmp2
1996 \ parity is odd
1997 \ (counterclockwise)
1998 A------- 0
1999 /
2000 /
2001 2 /
2002
2003 ------------------------------------*/
2004 bond_parity = 2 - ( tmp1 < tmp2 );
2005 for ( j = 0; j < 3; j ++ ) {
2006 z_dir[j] = (S_CHAR) (pnt[p2][j]>= 0.0? floor(0.5 + 100.0 * pnt[p2][j]) :
2007 -floor(0.5 - 100.0 * pnt[p2][j])); /* abs(z_dir) = 100 */
2008 }
2009 /* check for ambiguity */
2010 if ( nNumExplictAttachments > 2 ) {
2011 min_tmp = inchi_min( tmp1, tmp2 );
2012 max_tmp = inchi_max( tmp1, tmp2 );
2013 if ( min_tmp > one_pi-MIN_SINE || max_tmp < one_pi+MIN_SINE || max_tmp-min_tmp > one_pi - MIN_SINE ) {
2014 at[cur_at].bAmbiguousStereo |= AMBIGUOUS_STEREO;
2015 } else /* 3D ambiguity 8-28-2002 */
2016 if ( fabs(at_coord[0][2]) > MAX_SINE ) { /* all fabs(at_coord[j][2] (j=0..2) must be equal */
2017 at[cur_at].bAmbiguousStereo |= AMBIGUOUS_STEREO;
2018 }
2019 } else
2020 if ( nNumExplictAttachments == 2 ) { /* 10-6-2003: added */
2021 min_tmp = fabs(tmp1 - one_pi);
2022 if ( min_tmp < MIN_SINE ) {
2023 bond_parity = AB_PARITY_UNDF; /* consider as undefined 10-6-2003 */
2024 } else
2025 if ( min_tmp < MIN_ANGLE_DBOND ) {
2026 at[cur_at].bAmbiguousStereo |= AMBIGUOUS_STEREO;
2027 }
2028 }
2029
2030
2031 /* for 3 neighbors moving implicit H to the index=0 from index=2 position */
2032 /* can be done in 2 transpositions and does not change atom's parity */
2033 exit_function:
2034 if ( num_H > 1 && bond_parity > 0 && !(bond_parity & AB_PARITY_0D) /*&& PARITY_WELL_DEF(bond_parity)*/ ) {
2035 /*
2036 * stereo only if isotopes are counted. Do not inverse
2037 * Examples: sign for this:
2038 * H D
2039 * / / H
2040 * ==C or ==CH /
2041 * \ ==N (bValence3=1)
2042 * D
2043 * two explicit one explicit H isotope (D),
2044 * isotopic H atoms one implicit H
2045 */
2046 bond_parity = -bond_parity; /* refers to isotopically substituted structure only */
2047 }
2048 return bond_parity;
2049 }
2050
2051 /*************************************************************/
save_a_stereo_bond(int z_prod,int result_action,int at1,int ord1,AT_NUMB * stereo_bond_neighbor1,S_CHAR * stereo_bond_ord1,S_CHAR * stereo_bond_z_prod1,S_CHAR * stereo_bond_parity1,int at2,int ord2,AT_NUMB * stereo_bond_neighbor2,S_CHAR * stereo_bond_ord2,S_CHAR * stereo_bond_z_prod2,S_CHAR * stereo_bond_parity2)2052 int save_a_stereo_bond( int z_prod, int result_action,
2053 int at1, int ord1, AT_NUMB *stereo_bond_neighbor1, S_CHAR *stereo_bond_ord1, S_CHAR *stereo_bond_z_prod1, S_CHAR *stereo_bond_parity1,
2054 int at2, int ord2, AT_NUMB *stereo_bond_neighbor2, S_CHAR *stereo_bond_ord2, S_CHAR *stereo_bond_z_prod2, S_CHAR *stereo_bond_parity2 )
2055 {
2056 int i1, i2;
2057 for ( i1 = 0; i1 < MAX_NUM_STEREO_BONDS && stereo_bond_neighbor1[i1]; i1 ++ )
2058 ;
2059 for ( i2 = 0; i2 < MAX_NUM_STEREO_BONDS && stereo_bond_neighbor2[i2]; i2 ++ )
2060 ;
2061 if ( i1 == MAX_NUM_STEREO_BONDS || i2 == MAX_NUM_STEREO_BONDS )
2062 return 0;
2063
2064 stereo_bond_parity1[i1] =
2065 stereo_bond_parity2[i2] = result_action;
2066
2067 stereo_bond_neighbor1[i1] = (AT_NUMB) (at2+1);
2068 stereo_bond_ord1[i1] = (S_CHAR)ord1;
2069 stereo_bond_neighbor2[i2] = (AT_NUMB) (at1+1);
2070 stereo_bond_ord2[i2] = (S_CHAR)ord2;
2071 stereo_bond_z_prod1[i1] =
2072 stereo_bond_z_prod2[i2] = (S_CHAR)z_prod;
2073 return 1;
2074 }
2075 /***************************************************************/
get_allowed_stereo_bond_type(int bond_type)2076 int get_allowed_stereo_bond_type( int bond_type )
2077 {
2078 #if (ALLOW_TAUT_ATTACHMENTS_TO_STEREO_BONDS == 0 )
2079 if ( (bond_type & ~BOND_MARK_ALL) == BOND_TAUTOM )
2080 return 0; /* no tautomer bonds allowed */
2081 else
2082 #endif
2083 #if ( EXCL_ALL_AROM_BOND_PARITY == 1 ) /* { */
2084 /* a stereo bond cannot belong to an aromatic atom */
2085 if ( (bond_type &= ~BOND_MARK_ALL) == BOND_ALTERN )
2086 {
2087 return 0;
2088 }
2089 #else /* } { */
2090 #if ( ADD_6MEMB_AROM_BOND_PARITY == 1 )
2091 /* accept any aromatic bond as a stereo bond */
2092 if ( (bond_type &= ~BOND_MARK_ALL) == BOND_ALTERN )
2093 #else
2094 /* accept only aromatic bonds in non-6-member rings */
2095 if ( (bond_type &= ~BOND_MARK_ALL) == BOND_ALTERN ) )
2096 #endif
2097 {
2098 return BOND_ALTERN;
2099 }
2100 #endif /* } */
2101 else
2102 /* at this point BOND_MARK_ALL bits have been removed from bond_type */
2103 if ( bond_type == BOND_DOUBLE || bond_type == BOND_SINGLE ) {
2104 return bond_type;
2105 }
2106 #if (ALLOW_TAUT_ATTACHMENTS_TO_STEREO_BONDS == 1 )
2107 else
2108 if ( bond_type == BOND_TAUTOM ) {
2109 return BOND_TAUTOM;
2110 }
2111 #endif
2112
2113 return 0; /* wrong bond type */
2114 }
2115
2116 /*************************************************************/
can_be_a_stereo_bond_with_isotopic_H(inp_ATOM * at,int cur_at,INCHI_MODE nMode)2117 int can_be_a_stereo_bond_with_isotopic_H( inp_ATOM *at, int cur_at, INCHI_MODE nMode )
2118 {
2119 int i, j, next_at, num_stereo_bonds, bFound;
2120 int bond_type, num_2s, num_alt;
2121 int num_2s_next, num_alt_next, num_wrong_bonds_1, num_wrong_bonds_2;
2122 #if ( N_V_STEREOBONDS == 1 )
2123 int n2sh, num_2s_hetero[2], num_2s_hetero_next[2], next_next_at, type_N, type_N_next;
2124 #endif
2125 if ( MAX_NUM_STEREO_BOND_NEIGH < at[cur_at].valence+at[cur_at].num_H ||
2126 MIN_NUM_STEREO_BOND_NEIGH > at[cur_at].valence+at[cur_at].num_H )
2127 return 0;
2128 if ( !bCanAtomHaveAStereoBond( at[cur_at].elname, at[cur_at].charge, at[cur_at].radical ) )
2129 return 0;
2130 /* count bonds and find the second atom on the stereo bond */
2131 num_2s = num_alt = num_wrong_bonds_1 = 0;
2132 #if ( N_V_STEREOBONDS == 1 )
2133 num_2s_hetero[0] = num_2s_hetero[1] = type_N = 0;
2134 if ( 0 == at[cur_at].num_H && 0 == at[cur_at].charge && 0 == at[cur_at].radical &&
2135 3 == get_endpoint_valence( at[cur_at].el_number ) ) {
2136 if ( 2 == at[cur_at].valence && 3 == at[cur_at].chem_bonds_valence ) {
2137 type_N = 1;
2138 } else
2139 if ( 3 == at[cur_at].valence && 5 == at[cur_at].chem_bonds_valence ) {
2140 type_N = 2; /* unfortunately includes >N# */
2141 }
2142 }
2143 #endif
2144 for ( i = 0, num_stereo_bonds = 0; i < at[cur_at].valence; i ++ ) {
2145 bFound = 0;
2146 next_at = at[cur_at].neighbor[i];
2147 bond_type = get_allowed_stereo_bond_type( (int)at[cur_at].bond_type[i] );
2148 if ( bond_type == BOND_ALTERN ) {
2149 num_alt ++;
2150 if ( cur_at > next_at && !(nMode & CMODE_NO_ALT_SBONDS) )
2151 bFound = 1;
2152 } else
2153 if ( bond_type == BOND_DOUBLE ) {
2154 num_2s ++;
2155 #if ( N_V_STEREOBONDS == 1 )
2156 if ( 0 <= (n2sh = bIsSuitableHeteroInpAtom( at + next_at )) ) {
2157 num_2s_hetero[n2sh] ++; /* n2sh=0 -> =N- or =NH; n2sh=1 -> =O */
2158 }
2159 #endif
2160 if ( cur_at > next_at )
2161 bFound = 1;
2162 } else
2163 if ( bond_type != BOND_SINGLE && bond_type != BOND_TAUTOM ) {
2164 num_wrong_bonds_1 ++;
2165 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2166 if ( num_wrong_bonds_1 > 1 || (num_wrong_bonds_1 && 2 >= at[cur_at].valence) ) {
2167 return 0; /* wrong bond type */
2168 } else {
2169 continue;
2170 }
2171 #else
2172 return 0; /* wrong bond type */
2173 #endif
2174 }
2175
2176 if ( bFound ) {
2177 /* check "next_at" atom on the opposite side of the bond */
2178 if ( MAX_NUM_STEREO_BOND_NEIGH < at[next_at].valence+at[next_at].num_H ||
2179 MIN_NUM_STEREO_BOND_NEIGH > at[next_at].valence+at[next_at].num_H )
2180 continue;
2181 if ( !bCanAtomHaveAStereoBond( at[next_at].elname, at[next_at].charge, at[next_at].radical ) )
2182 continue;
2183 /* next atom neighbors */
2184 num_2s_next = num_alt_next = num_wrong_bonds_2 = 0;
2185 #if ( N_V_STEREOBONDS == 1 )
2186 num_2s_hetero_next[0] = num_2s_hetero_next[1] = type_N_next = 0;
2187 if ( 0 == at[next_at].num_H && 0 == at[next_at].charge && 0 == at[next_at].radical &&
2188 3 == get_endpoint_valence( at[next_at].el_number ) ) {
2189 if ( 2 == at[next_at].valence && 3 == at[next_at].chem_bonds_valence ) {
2190 type_N_next = 1; /* -N= */
2191 } else
2192 if ( 3 == at[next_at].valence && 5 == at[next_at].chem_bonds_valence ) {
2193 type_N_next = 2; /* unfortunately includes >N# */
2194 }
2195 }
2196 #endif
2197 for ( j = 0; j < at[next_at].valence; j ++ ) {
2198 bond_type = get_allowed_stereo_bond_type( (int)at[next_at].bond_type[j] );
2199 if ( bond_type == BOND_ALTERN )
2200 num_alt_next ++;
2201 else
2202 if ( bond_type == BOND_DOUBLE ) {
2203 num_2s_next ++;
2204 #if ( N_V_STEREOBONDS == 1 )
2205 next_next_at = at[next_at].neighbor[j];
2206 if ( 0 <= (n2sh = bIsSuitableHeteroInpAtom( at + next_next_at )) ) {
2207 num_2s_hetero_next[n2sh] ++; /* n2sh=0 -> =N- or =NH; n2sh=1 -> =O */
2208 }
2209 #endif
2210 } else
2211 if ( bond_type != BOND_SINGLE && bond_type != BOND_TAUTOM ) {
2212 num_wrong_bonds_2 ++;
2213 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2214 if ( num_wrong_bonds_1 > 1 || (num_wrong_bonds_1 && 2 >= at[cur_at].valence) ) {
2215 break; /* wrong bond type */
2216 } else {
2217 continue;
2218 }
2219 #else
2220 break; /* wrong bond type */
2221 #endif
2222 }
2223 }
2224 /* figure out whether the at[cur_at]--at[next_at] bond may not be stereogenic */
2225
2226 #if ( N_V_STEREOBONDS == 1 )
2227 if ( 3 == (type_N | type_N_next) &&
2228 ( (2 == type_N && !bIsOxide( at, cur_at )) ||
2229 (2 == type_N_next && !bIsOxide( at, next_at )) ) ) {
2230 bFound = 0;
2231 } else
2232 #endif
2233 if ( j < at[next_at].valence || /* at[next_at] has a wrong bond type*/
2234 (num_alt_next>0) + (num_2s_next>0) != 1 /* only one type of stereogenic bond permitted */
2235 ) {
2236 bFound = 0;
2237 } else
2238 if ( 2 < num_2s_next ) {
2239 bFound = 0;
2240 } else
2241 if ( 2 == num_2s_next ) {
2242 if ( 2 == at[next_at].valence ) {
2243 ; /* only one double bond permitted except cumulenes */
2244 #if ( N_V_STEREOBONDS == 1 )
2245 } else
2246 if ( 1 == (num_2s_hetero_next[0] | num_2s_hetero_next[1]) &&
2247 3 == at[next_at].valence + at[next_at].num_H &&
2248 5 == at[next_at].chem_bonds_valence + at[next_at].num_H &&
2249 3 == get_endpoint_valence( at[next_at].el_number ) &&
2250 (!type_N || bIsOxide( at, next_at )) ) {
2251 ; /*
2252 * found:
2253 *
2254 * \ / \ / \ /
2255 * \ / \ / \ /
2256 * N==C or N==C or N==N
2257 * // \ // \ // \
2258 * O ^ \ N ^ \ O ^ \
2259 * | | |
2260 * | | |
2261 * at[next_at] at[next_at] at[next_at]
2262 */
2263 #endif
2264 } else {
2265 bFound = 0;
2266 }
2267 }
2268
2269 }
2270 if ( bFound ) {
2271 num_stereo_bonds++;
2272 }
2273 }
2274
2275 if ( (num_alt>0) + (num_2s>0) != 1 || !num_stereo_bonds )
2276 return 0;
2277 if ( num_2s > 1 ) {
2278 #if ( N_V_STEREOBONDS == 1 )
2279 if ( 2 == num_2s &&
2280 1 == (num_2s_hetero[0] | num_2s_hetero[1]) &&
2281 3 == at[cur_at].valence + at[cur_at].num_H &&
2282 5 == at[cur_at].chem_bonds_valence + at[cur_at].num_H &&
2283 3 == get_endpoint_valence( at[cur_at].el_number ) ) {
2284 ;
2285 } else {
2286 return 0;
2287 }
2288 #else
2289 return 0;
2290 #endif
2291 }
2292
2293 return num_stereo_bonds;
2294 }
2295 /*************************************************************/
half_stereo_bond_action(int nParity,int bUnknown,int bIsotopic,int vABParityUnknown)2296 int half_stereo_bond_action( int nParity, int bUnknown, int bIsotopic, int vABParityUnknown )
2297 {
2298 #define AB_NEGATIVE 0x10
2299 #define AB_UNKNOWN 0x20
2300 int nAction;
2301
2302 if ( nParity == AB_PARITY_NONE )
2303 return AB_PARITY_NONE;
2304
2305 /* Unknown (type 1) in the parity value may come from the 'Either' single bond only */
2306 /* Treat it as a known single bond geometry and unknown (Either) double bond */
2307 if ( nParity == vABParityUnknown /*AB_PARITY_UNKN*/ )
2308 nParity = AB_PARITY_ODD | AB_UNKNOWN;
2309 if ( nParity == -vABParityUnknown /*AB_PARITY_UNKN*/ )
2310 nParity = AB_PARITY_ODD | AB_UNKNOWN | AB_NEGATIVE;
2311
2312 /* make positive, replace AB_PARITY_EVEN with AB_PARITY_ODD */
2313 if ( nParity < 0 )
2314 nParity = ((nParity == -AB_PARITY_EVEN)? AB_PARITY_ODD : (-nParity)) | AB_NEGATIVE;
2315 else
2316 if (nParity == AB_PARITY_EVEN)
2317 nParity = AB_PARITY_ODD;
2318
2319 /* Unknown (type 2): was detected in the double bond attribute */
2320 /* (this 'unknown' came from 'Either' double bond) */
2321 /* Treat both unknowns in the same way */
2322 if ( bUnknown )
2323 nParity |= AB_UNKNOWN;
2324
2325 if ( bIsotopic ) {
2326 switch ( nParity ) {
2327 case AB_PARITY_ODD:
2328 case AB_PARITY_ODD | AB_NEGATIVE:
2329 nAction = AB_PARITY_CALC;
2330 break;
2331 case AB_PARITY_ODD | AB_UNKNOWN:
2332 case AB_PARITY_UNDF | AB_UNKNOWN:
2333 case AB_PARITY_ODD | AB_UNKNOWN | AB_NEGATIVE:
2334 case AB_PARITY_UNDF | AB_UNKNOWN | AB_NEGATIVE:
2335 nAction = vABParityUnknown /*AB_PARITY_UNKN*/;
2336 break;
2337 case AB_PARITY_IISO:
2338 case AB_PARITY_IISO | AB_UNKNOWN:
2339 nAction = AB_PARITY_NONE;
2340 break;
2341 case AB_PARITY_UNDF:
2342 case AB_PARITY_UNDF | AB_NEGATIVE:
2343 nAction = AB_PARITY_UNDF;
2344 break;
2345 default:
2346 nAction = -1; /* program error */
2347 }
2348 } else {
2349 /* Non-isotopic */
2350 switch ( nParity ) {
2351 case AB_PARITY_ODD:
2352 nAction = AB_PARITY_CALC;
2353 break;
2354 case AB_PARITY_ODD | AB_UNKNOWN:
2355 case AB_PARITY_UNDF | AB_UNKNOWN:
2356 nAction = vABParityUnknown /*AB_PARITY_UNKN*/;
2357 break;
2358 /* case AB_PARITY_ODD | AB_UNKNOWN | AB_NEGATIVE: */
2359 case AB_PARITY_UNDF:
2360 nAction = AB_PARITY_UNDF;
2361 break;
2362 case AB_PARITY_ODD | AB_UNKNOWN | AB_NEGATIVE:
2363 case AB_PARITY_ODD | AB_NEGATIVE:
2364 case AB_PARITY_IISO:
2365 case AB_PARITY_IISO | AB_UNKNOWN:
2366 case AB_PARITY_UNDF | AB_NEGATIVE:
2367 case AB_PARITY_UNDF | AB_UNKNOWN | AB_NEGATIVE:
2368 nAction = AB_PARITY_NONE;
2369 break;
2370 default:
2371 nAction = -1; /* program error */
2372 }
2373 }
2374 return nAction;
2375 #undef AB_NEGATIVE
2376 #undef AB_UNKNOWN
2377 }
2378 /*************************************************************/
set_stereo_bonds_parity(sp_ATOM * out_at,inp_ATOM * at,int at_1,inp_ATOM * at_removed_H,int num_removed_H,INCHI_MODE nMode,QUEUE * q,AT_RANK * nAtomLevel,S_CHAR * cSource,AT_RANK min_sb_ring_size,int bPointedEdgeStereo,int vABParityUnknown)2379 int set_stereo_bonds_parity( sp_ATOM *out_at, inp_ATOM *at, int at_1, inp_ATOM *at_removed_H, int num_removed_H,
2380 INCHI_MODE nMode, QUEUE *q, AT_RANK *nAtomLevel, S_CHAR *cSource,
2381 AT_RANK min_sb_ring_size, int bPointedEdgeStereo, int vABParityUnknown )
2382 {
2383 int j, k, next_at_1, i_next_at_1, i_next_at_2, at_2, next_at_2, num_stereo_bonds, bFound, bAllene;
2384 int bond_type, num_2s_1, num_alt_1;
2385 int num_2s_2, num_alt_2;
2386 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2387 int num_wrong_bonds_1, num_wrong_bonds_2;
2388 #endif
2389 #if ( N_V_STEREOBONDS == 1 )
2390 int n2sh, num_2s_hetero[2], num_2s_hetero_next[2], next_next_at, type_N, type_N_next;
2391 #endif
2392 int num_stored_stereo_bonds, num_stored_isotopic_stereo_bonds;
2393 int chain_length, num_chains, cur_chain_length;
2394 int all_at_2[MAX_NUM_STEREO_BONDS];
2395 int all_pos_1[MAX_NUM_STEREO_BONDS], all_pos_2[MAX_NUM_STEREO_BONDS];
2396 S_CHAR all_unkn[MAX_NUM_STEREO_BONDS];
2397 int /*at_1_parity, at_2_parity,*/ nUnknown, stop=0;
2398
2399 /* at_1_parity = AB_PARITY_NONE; */ /* do not know */
2400 /* check valence */
2401 if ( MAX_NUM_STEREO_BOND_NEIGH < at[at_1].valence+at[at_1].num_H ||
2402 MIN_NUM_STEREO_BOND_NEIGH > at[at_1].valence+at[at_1].num_H )
2403 return 0;
2404 if ( !bCanAtomHaveAStereoBond( at[at_1].elname, at[at_1].charge, at[at_1].radical ) )
2405 return 0;
2406 if ( at[at_1].c_point )
2407 return 0; /* rejects atoms that can lose or gain a (positive) charge. 01-24-2003 */
2408
2409 /* middle cumulene atoms, for example, =C=, should be ignored here */
2410 /* only atoms at the ends of cumulene chains are considered. */
2411 if ( !at[at_1].num_H && 2 == at[at_1].valence &&
2412 BOND_DOUBLE == get_allowed_stereo_bond_type( (int)at[at_1].bond_type[0] ) &&
2413 BOND_DOUBLE == get_allowed_stereo_bond_type( (int)at[at_1].bond_type[1] ) ) {
2414 return 0;
2415 }
2416
2417 /* count bonds and find the second atom on the stereo bond */
2418 num_2s_1 = num_alt_1 = 0;
2419 chain_length = 0;
2420 num_chains = 0;
2421 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2422 num_wrong_bonds_1 = 0;
2423 #endif
2424 #if ( N_V_STEREOBONDS == 1 )
2425 num_2s_hetero[0] = num_2s_hetero[1] = type_N = 0;
2426 if ( 0 == at[at_1].num_H && 0 == at[at_1].charge && 0 == at[at_1].radical &&
2427 3 == get_endpoint_valence( at[at_1].el_number ) ) {
2428 if ( 2 == at[at_1].valence && 3 == at[at_1].chem_bonds_valence ) {
2429 type_N = 1;
2430 } else
2431 if ( 3 == at[at_1].valence && 5 == at[at_1].chem_bonds_valence ) {
2432 type_N = 2; /* unfortunately includes >N# */
2433 }
2434 }
2435 #endif
2436 for ( i_next_at_1 = 0, num_stereo_bonds = 0; i_next_at_1 < at[at_1].valence; i_next_at_1 ++ ) {
2437 nUnknown = (at[at_1].bond_stereo[i_next_at_1] == STEREO_DBLE_EITHER);
2438 bond_type = get_allowed_stereo_bond_type( (int)at[at_1].bond_type[i_next_at_1] );
2439 at_2 = -1; /* not found */
2440 if ( bond_type == BOND_ALTERN ||
2441 bond_type == BOND_DOUBLE ) {
2442 next_at_1 = at_2 = at[at_1].neighbor[i_next_at_1];
2443 next_at_2 = at_1;
2444 }
2445 switch ( bond_type ) {
2446 case BOND_ALTERN:
2447 num_alt_1 ++;
2448 #if ( FIND_RING_SYSTEMS == 1 )
2449 if ( at[at_1].nRingSystem != at[at_2].nRingSystem )
2450 continue; /* reject alt. bond connecting different ring systems */
2451 #endif
2452 if ( (nMode & CMODE_NO_ALT_SBONDS) ||
2453 !bCanAtomHaveAStereoBond( at[at_2].elname, at[at_2].charge, at[at_2].radical ) ) {
2454 continue; /* reject non-stereogenic bond to neighbor ord. #i_next_at_1 */
2455 }
2456 break;
2457 case BOND_DOUBLE:
2458 /* check for cumulene/allene */
2459 num_2s_1++;
2460 cur_chain_length = 0;
2461 if ( bCanAtomBeTerminalAllene( at[at_1].elname, at[at_1].charge, at[at_1].radical ) ) {
2462 /*
2463 * Example of cumulene
2464 * chain length = 2: >X=C=C=Y<
2465 * | | | |
2466 * 1st cumulene atom= at_1 | | at_2 =last cumlene chain atom
2467 * next to at_1= next_at_1 next_at_2 =previous to at_2
2468 *
2469 * chain length odd: stereocenter on the middle atom ( 1=> allene )
2470 * chain length even: "long stereogenic bond"
2471 */
2472 while ((bAllene =
2473 !at[at_2].num_H && at[at_2].valence == 2 &&
2474 BOND_DOUBLE == get_allowed_stereo_bond_type( (int)at[at_2].bond_type[0] ) &&
2475 BOND_DOUBLE == get_allowed_stereo_bond_type( (int)at[at_2].bond_type[1] )) &&
2476 bCanAtomBeMiddleAllene( at[at_2].elname, at[at_2].charge, at[at_2].radical ) ) {
2477 k = ((int)at[at_2].neighbor[0]==next_at_2); /* opposite neighbor position */
2478 next_at_2 = at_2;
2479 nUnknown += (at[at_2].bond_stereo[k] == STEREO_DBLE_EITHER);
2480 at_2 = (int)at[at_2].neighbor[k];
2481 cur_chain_length ++; /* count =C= atoms */
2482 }
2483 if ( cur_chain_length ) {
2484 num_chains ++;
2485 if ( bAllene /* at the end of the chain atom Y is =Y=, not =Y< or =Y- */ ||
2486 !bCanAtomBeTerminalAllene( at[at_2].elname, at[at_2].charge, at[at_2].radical ) ) {
2487 cur_chain_length = 0;
2488 continue; /* ignore: does not fit cumulene description; go to check next at_1 neighbor */
2489 }
2490 chain_length = cur_chain_length; /* accept a stereogenic cumulele */
2491 }
2492 }
2493 #if ( N_V_STEREOBONDS == 1 )
2494 if ( !cur_chain_length &&
2495 0 <= (n2sh = bIsSuitableHeteroInpAtom( at + at_2 )) ) {
2496 num_2s_hetero[n2sh] ++; /* n2sh=0 -> =N- or =NH; n2sh=1 -> =O */
2497 }
2498 #endif
2499 if ( !cur_chain_length &&
2500 !bCanAtomHaveAStereoBond( at[at_2].elname, at[at_2].charge, at[at_2].radical ) ) {
2501 continue; /* reject non-stereogenic bond to neighbor #i_next_at_1 */
2502 }
2503
2504 break;
2505
2506 case BOND_SINGLE:
2507 case BOND_TAUTOM:
2508 continue; /* reject non-stereogenic bond to neighbor #i_next_at_1 */
2509 default:
2510 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2511 num_wrong_bonds_1 ++;
2512 continue;
2513 #else
2514 return 0; /* wrong bond type; */
2515 #endif
2516 }
2517
2518 /* check atom at the opposite end of possibly stereogenic bond */
2519
2520 bFound = (at_2 >= 0 && at_1 > at_2 ); /* i_next_at_1 = at_1 stereogenic bond neighbor attachment number */
2521
2522 if ( bFound ) {
2523 /* check "at_2" atom on the opposite side of the bond or cumulene chain */
2524 if ( MAX_NUM_STEREO_BOND_NEIGH < at[at_2].valence+at[at_2].num_H ||
2525 MIN_NUM_STEREO_BOND_NEIGH > at[at_2].valence+at[at_2].num_H )
2526 continue;
2527
2528 /* check at_2 neighbors and bonds */
2529 num_2s_2 = num_alt_2 = 0;
2530 #if ( N_V_STEREOBONDS == 1 )
2531 num_2s_hetero_next[0] = num_2s_hetero_next[1] = type_N_next = 0;
2532 if ( 0 == at[at_2].num_H && 0 == at[at_2].charge && 0 == at[at_2].radical &&
2533 3 == get_endpoint_valence( at[at_2].el_number ) ) {
2534 if ( 2 == at[at_2].valence && 3 == at[at_2].chem_bonds_valence ) {
2535 type_N_next = 1; /* -N= */
2536 } else
2537 if ( 3 == at[at_2].valence && 5 == at[at_2].chem_bonds_valence ) {
2538 type_N_next = 2; /* unfortunately includes >N# */
2539 }
2540 }
2541 #endif
2542 i_next_at_2 = -1; /* unassigned mark */
2543 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2544 num_wrong_bonds_2 = 0;
2545 #endif
2546 for ( j = 0; j < at[at_2].valence; j ++ ) {
2547 bond_type = get_allowed_stereo_bond_type( (int)at[at_2].bond_type[j] );
2548 if ( !bond_type ) {
2549 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2550 num_wrong_bonds_2 ++;
2551 continue; /* this bond type is not allowed to be adjacent to a stereo bond */
2552 #else
2553 break;
2554 #endif
2555 }
2556 if ( bond_type == BOND_DOUBLE ) {
2557 num_2s_2 ++;
2558 #if ( N_V_STEREOBONDS == 1 )
2559 next_next_at = at[at_2].neighbor[j];
2560 if ( 0 <= (n2sh = bIsSuitableHeteroInpAtom( at + next_next_at )) ) {
2561 num_2s_hetero_next[n2sh] ++; /* n2sh=0 -> =N- or =NH; n2sh=1 -> =O */
2562 }
2563 #endif
2564 } else {
2565 num_alt_2 += ( bond_type == BOND_ALTERN );
2566 }
2567 if ( (int)at[at_2].neighbor[j] == next_at_2 )
2568 i_next_at_2 = j; /* assigned */
2569 }
2570 if (
2571 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2572 num_wrong_bonds_2 > 1 || (num_wrong_bonds_2 && 2 >= at[at_2].valence) ||
2573 #else
2574 j < at[at_2].valence /* "next" has a wrong bond type*/ ||
2575 #endif
2576 (num_alt_2>0) + (num_2s_2>0) != 1 || /* all double XOR all alt bonds only */
2577 /* num_2s_2 > 1 ||*/ /* only one double bond permitted */
2578 i_next_at_2 < 0 /* atom next to the opposite atom not found */ ) {
2579 bFound = 0;
2580 } else
2581 if ( at[at_2].c_point ) {
2582 bFound = 0; /* rejects atoms that can lose or gain a (positive) charge. 01-24-2003 */
2583 } else
2584 if ( num_2s_2 > 2 ) {
2585 bFound = 0;
2586 } else
2587 #if ( N_V_STEREOBONDS == 1 )
2588 if ( 3 == (type_N | type_N_next) &&
2589 ( (2 == type_N && !bIsOxide( at, at_1 )) ||
2590 (2 == type_N_next && !bIsOxide( at, at_2 )) ) ) {
2591 bFound = 0;
2592 } else
2593 #endif
2594 if ( 2 == num_2s_2 ) {
2595 #if ( N_V_STEREOBONDS == 1 )
2596 if ( !chain_length &&
2597 1 == (num_2s_hetero_next[0] | num_2s_hetero_next[1]) &&
2598 3 == at[at_2].valence + at[at_2].num_H &&
2599 5 == at[at_2].chem_bonds_valence + at[at_2].num_H &&
2600 3 == get_endpoint_valence( at[at_2].el_number ) &&
2601 (!type_N || bIsOxide( at, at_2 )) ) {
2602 /*
2603 * found:
2604 *
2605 * \ / \ / \ /
2606 * \ / \ / \ /
2607 * N==C or N==C or N==N
2608 * // \ // \ // \
2609 * O ^ \ N ^ \ O ^ \
2610 * | | |
2611 * | | |
2612 * at[at_2] at[at_2] at[at_2]
2613 */
2614 ;
2615 } else {
2616 bFound = 0;
2617 }
2618 #else
2619 bFound = 0;
2620 #endif
2621 }
2622
2623
2624 if ( chain_length && num_alt_2 )
2625 return 0; /* allow no alt bonds in cumulenes */
2626 }
2627 if ( bFound ) {
2628 all_pos_1[num_stereo_bonds] = i_next_at_1; /* neighbor to at_1 position */
2629 all_pos_2[num_stereo_bonds] = i_next_at_2; /* neighbor to at_2 position */
2630 all_at_2[num_stereo_bonds] = at_2; /* at_2 */
2631 all_unkn[num_stereo_bonds] = nUnknown; /* stereogenic bond has Unknown configuration */
2632 /*
2633 if ( (at[at_1].bUsed0DParity & 2) || (at[at_2].bUsed0DParity & 2) ) {
2634 for ( k = 0; k < MAX_NUM_STEREO_BONDS && at[at_1].sb_parity[k]; k ++ ) {
2635 if ( at[at_1].sb_neigh[k] == i_next_at_1 ) {
2636 if ( at[at_1].sb_parity[k] == AB_PARITY_UNKN && !nUnknown ) {
2637 all_unkn[num_stereo_bonds] = 1;
2638 }
2639 break;
2640 }
2641 }
2642 }
2643 */
2644 num_stereo_bonds ++;
2645 }
2646 }
2647 if ( num_chains > 1 ) {
2648 return 0; /* cannot be more than 1 cumulene chain. */
2649 }
2650 #if ( ONE_BAD_SB_NEIGHBOR == 1 )
2651 if ( num_wrong_bonds_1 > 1 || (num_wrong_bonds_1 && 2 >= at[at_1].valence) ) {
2652 return 0; /* wrong bond type */
2653 }
2654 #endif
2655 /* accept only short chains for now */
2656 /* chain_length=1: >C=C=C< tetrahedral center, allene */
2657 /* chain_length=2: >C=C=C=C< stereogenic bond, cumulene */
2658 if ( chain_length && (num_stereo_bonds != 1 || num_alt_1 || chain_length > MAX_CUMULENE_LEN) ) {
2659 return 0;
2660 }
2661
2662 /* we need 1 double bond/chain XOR up to 3 arom. bonds */
2663 /* to have a stereogenic bond */
2664 if ( (num_alt_1>0) + (num_2s_1>0) != 1 || !num_stereo_bonds /*|| num_2s_1 > 1*/ )
2665 return 0;
2666
2667 if ( num_2s_1 > 1 ) {
2668 #if ( N_V_STEREOBONDS == 1 )
2669 if ( 2 == num_2s_1 &&
2670 2 == type_N &&
2671 1 == (num_2s_hetero[0] | num_2s_hetero[1]) &&
2672 3 == at[at_1].valence + at[at_1].num_H &&
2673 5 == at[at_1].chem_bonds_valence + at[at_1].num_H &&
2674 3 == get_endpoint_valence( at[at_1].el_number ) ) {
2675 ;
2676 } else {
2677 return 0;
2678 }
2679 #else
2680 return 0;
2681 #endif
2682 }
2683
2684 /* ================== calculate parities ====================== */
2685
2686
2687 /* find possibly stereo bonds and save them */
2688 num_stored_isotopic_stereo_bonds = 0;
2689 num_stored_stereo_bonds = 0;
2690 for ( k = 0; k < num_stereo_bonds; k ++ ) {
2691
2692 int cur_parity, next_parity, abs_cur_parity, abs_next_parity, dot_prod_z;
2693 S_CHAR z_dir1[3], z_dir2[3]; /* 3D vectors for half stereo bond parity direction */
2694 int chain_len_bits = MAKE_BITS_CUMULENE_LEN(chain_length);
2695 int cur_parity_defined, next_parity_defined;
2696 int cur_action, next_action, result_action;
2697
2698 at_2 = all_at_2[k];
2699 i_next_at_1 = all_pos_1[k];
2700
2701 #if ( MIN_SB_RING_SIZE > 0 )
2702 if( at[at_1].nRingSystem == at[at_2].nRingSystem ) {
2703 /* check min. ring size only if both double bond/cumulene */
2704 /* ending atoms belong to the same ring system */
2705 j = is_bond_in_Nmax_memb_ring( at, at_1, i_next_at_1, q, nAtomLevel, cSource, min_sb_ring_size );
2706 if ( j > 0 ) {
2707 continue;
2708 } else
2709 if ( j < 0 ) {
2710 return CT_STEREOBOND_ERROR;
2711 }
2712 }
2713 #endif
2714
2715 i_next_at_2 = all_pos_2[k];
2716 nUnknown = all_unkn[k];
2717 memset(z_dir1, 0, sizeof(z_dir1));
2718 memset(z_dir2, 0, sizeof(z_dir2));
2719
2720 /********************************************************************************
2721 * find atom parities (negative means parity due to H-isotopes only)
2722 * and half stereo bond parity directions z_dir1, z_dir2.
2723 *
2724 * Bond can have unknown or undefined parity or no parity because of:
2725 * 1. Geometry (poorly defined, cannot calculate, for example linear =C-F
2726 * or =CHD with no geometry) -- Undefined parity
2727 * H
2728 * 2. Identical H atoms (no parity in principle, for example =C< )
2729 * -- No parity H
2730 *
2731 * 3. The user said double bond stereo is unknown
2732 * or at least one of single bonds is in unknown direction
2733 * -- Unknown parity
2734 *
2735 * These 3 cases (see above) are referred below as 1, 2, 3.
2736 * Each of the cases may be present or not (2 possibilities)
2737 * Total number of combination is 2*2*2=8
2738 *
2739 * Since a case when all 3 are not present is a well-defined parity,
2740 * we do not consider this case here. Then 2*2*2-1=7 cases are left.
2741 *
2742 * If several cases are present, list them below separated by "+".
2743 * For example, 1+2 means (1) undefined geometry and (2) no parity
2744 * is possible because of identical H atoms.
2745 *
2746 * N) Decision table, Non-isotopic, 2*2*2-1=7 cases:
2747 * =================================================
2748 * none : 2+any: 1+2(e.g.=CH2); 1+2+3; 2; 2+3 AB_PARITY_NONE=0
2749 * undefined: 1 AB_PARITY_UNDF
2750 * unknown : 1+3; 3 AB_PARITY_UNKN
2751 *
2752 * I) Decision table, Isotopic, 2*2*2-1=7 cases:
2753 * =============================================
2754 * none : none
2755 * undefined: 1; 1+2; 1+2+3; 2; 2+3
2756 * unknown : 1+3; 3
2757 *
2758 * Note: When defining identical atoms H atoms in case 2,
2759 * Isotopic and Non-isotopic cases are different:
2760 * N: do NOT take into account the isotopic composition of H atoms
2761 * I: DO take into account the isotopic composition of H atoms
2762 * (it is assumed that H isotopes are always different)
2763 *
2764 * half_stereo_bond_parity() returns:
2765 * ==================================
2766 * Note: half_stereo_bond_parity() is unaware of case 3.
2767 *
2768 * can't be a half of a stereo bond AB_PARITY_NONE
2769 * 1, isotopic & non-isotopic: AB_PARITY_UNDF
2770 * 1, isotopic only -AB_PARITY_UNDF
2771 * 2, no parity: identical H isotopes AB_PARITY_IISO
2772 * 3, 'Either' single bond(s) AB_PARITY_UNKN ???
2773 * 3, 'Either' single bond(s), iso H -AB_PARITY_UNKN ???
2774 * defined parity AB_PARITY_ODD, AB_PARITY_EVEN
2775 * defined parity for isotopic only: -AB_PARITY_ODD, -AB_PARITY_EVEN
2776 *
2777 * Resultant value for the stereo bond parity
2778 * ---+-------------------+-------+--------+----------------+
2779 * 3? | half_stereo_bond_ | N or I| case 1,| bond parity |
2780 * | parity()= | | 2 or 3 | |
2781 * ---+-------------------+-------+--------+----------------+
2782 * ( AB_PARITY_ODD/EVEN) => N&I: - => AB_PARITY_CALC (=6, calc.later)
2783 * 3+( AB_PARITY_ODD/EVEN) => N&I: 3 => AB_PARITY_UNKN (=3)
2784 * (-AB_PARITY_ODD/EVEN) => N: 2 => AB_PARITY_NONE (=0)
2785 * (-AB_PARITY_ODD/EVEN) => I: - => AB_PARITY_CALC
2786 * 3+(-AB_PARITY_ODD/EVEN) => N: 2+3 => AB_PARITY_UNDF (=4)
2787 * 3+(-AB_PARITY_ODD/EVEN) => I: 3 => AB_PARITY_UNKN
2788 * ( AB_PARITY_IISO ) => N: 1+2, 2 => AB_PARITY_NONE (=0)
2789 * ( AB_PARITY_IISO ) => I: 1+2, 2 => AB_PARITY_UNDF
2790 * 3+( AB_PARITY_IISO ) => N: 1+2+3,2+3=> AB_PARITY_NONE
2791 * 3+( AB_PARITY_IISO ) => I: 1+2+3,2+3=> AB_PARITY_UNDF
2792 * ( AB_PARITY_UNDF ) => N&I: 1 => AB_PARITY_UNDF
2793 * 3+( AB_PARITY_UNDF ) => N&I: 1+3 => AB_PARITY_UNKN
2794 * (-AB_PARITY_UNDF ) => N: 1+2 => AB_PARITY_NONE
2795 * (-AB_PARITY_UNDF ) => I: 1 => AB_PARITY_UNDF
2796 * 3+(-AB_PARITY_UNDF ) => N: 1+2+3 => AB_PARITY_NONE
2797 * 3+(-AB_PARITY_UNDF ) => I: 1+3 => AB_PARITY_UNKN
2798 * ---+-------------------+-------+--------+----------------+
2799
2800 * If bond parity is undefined because abs(dot_prod_z) < MIN_DOT_PROD
2801 * then replace: AB_PARITY_CALC
2802 * with: AB_PARITY_UNDF
2803 * Joining two half_bond_parity() results:
2804 *
2805 *
2806 * atom1 \ atom2 | AB_PARITY_NONE AB_PARITY_UNKN AB_PARITY_UNDF AB_PARITY_CALC
2807 * ----------------+---------------------------------------------------------------
2808 *0=AB_PARITY_NONE | AB_PARITY_NONE AB_PARITY_NONE AB_PARITY_NONE AB_PARITY_NONE
2809 *3=AB_PARITY_UNKN | AB_PARITY_UNKN AB_PARITY_UNKN AB_PARITY_UNKN
2810 *4=AB_PARITY_UNDF | AB_PARITY_UNDF AB_PARITY_UNDF
2811 *6=AB_PARITY_CALC | AB_PARITY_CALC
2812 *
2813 * that is, take min out of the two
2814 *********************************************************************************/
2815
2816 cur_parity = half_stereo_bond_parity( at, at_1, at_removed_H, num_removed_H,
2817 z_dir1, bPointedEdgeStereo, vABParityUnknown );
2818 next_parity = half_stereo_bond_parity( at, at_2, at_removed_H, num_removed_H,
2819 z_dir2, bPointedEdgeStereo, vABParityUnknown );
2820
2821 if ( RETURNED_ERROR(cur_parity) || RETURNED_ERROR(next_parity) ) {
2822 return CT_CALC_STEREO_ERR;
2823 }
2824 if ( (at[at_1].bUsed0DParity & FlagSB_0D) || (at[at_2].bUsed0DParity & FlagSB_0D) ) {
2825 FixSb0DParities( at, /* at_removed_H, num_removed_H,*/ chain_length,
2826 at_1, i_next_at_1, z_dir1,
2827 at_2, i_next_at_2, z_dir2, &cur_parity, &next_parity );
2828 }
2829
2830 if ( cur_parity == AB_PARITY_NONE || abs(cur_parity) == AB_PARITY_IISO ) {
2831 continue;
2832 }
2833 if ( next_parity == AB_PARITY_NONE || abs(next_parity) == AB_PARITY_IISO ) {
2834 continue;
2835 }
2836
2837 cur_action = half_stereo_bond_action( cur_parity, nUnknown, 0, vABParityUnknown ); /* -1 => program error */
2838 next_action = half_stereo_bond_action( next_parity, nUnknown, 0, vABParityUnknown );
2839 result_action = inchi_min(cur_action, next_action);
2840
2841 if ( result_action == -1 ) {
2842 stop = 1; /* program error <BRKPT> */
2843 }
2844
2845 abs_cur_parity = abs( cur_parity );
2846 abs_next_parity = abs( next_parity );
2847 cur_parity_defined = ATOM_PARITY_WELL_DEF(abs_cur_parity);
2848 next_parity_defined = ATOM_PARITY_WELL_DEF(abs_next_parity);
2849
2850
2851 if ( cur_parity_defined && next_parity_defined ) {
2852 /* find how the whole bond parity depend on geometry */
2853 /* if dot_prod_z < 0 then bond_parity := 3-bond_parity */
2854 /* can be done only for a well-defined geometry */
2855 /*
2856 dot_prod_z = (chain_len_bits & BIT_CUMULENE_CHI)?
2857 triple_prod_char( at, at_1, i_next_at_1, z_dir1, at_2, i_next_at_2, z_dir2 ) :
2858 dot_prodchar3(z_dir1, z_dir2);
2859 */
2860 dot_prod_z = (chain_len_bits && BOND_CHAIN_LEN(chain_len_bits)%2)?
2861 triple_prod_char( at, at_1, i_next_at_1, z_dir1, at_2, i_next_at_2, z_dir2 ) :
2862 dot_prodchar3(z_dir1, z_dir2);
2863
2864 if ( abs(dot_prod_z) < MIN_DOT_PROD ) {
2865 /* The geometry is not well-defined. Eliminate AB_PARITY_CALC */
2866 result_action = inchi_min( result_action, AB_PARITY_UNDF );
2867 }
2868 } else {
2869 dot_prod_z = 0;
2870 }
2871
2872 if ( result_action != AB_PARITY_NONE && result_action != -1 ) {
2873 /* stereo, no isotopes (only positive) */
2874 if ( cur_parity > 0 && next_parity > 0 ) {
2875 if ( save_a_stereo_bond( dot_prod_z, result_action | chain_len_bits,
2876 at_1, i_next_at_1, out_at[at_1].stereo_bond_neighbor,
2877 out_at[at_1].stereo_bond_ord, out_at[at_1].stereo_bond_z_prod,
2878 out_at[at_1].stereo_bond_parity,
2879 at_2, i_next_at_2, out_at[at_2].stereo_bond_neighbor,
2880 out_at[at_2].stereo_bond_ord, out_at[at_2].stereo_bond_z_prod,
2881 out_at[at_2].stereo_bond_parity) ) {
2882 if ( !out_at[at_1].parity ||
2883 (cur_parity_defined && !ATOM_PARITY_WELL_DEF(abs(out_at[at_1].parity))) ) {
2884 out_at[at_1].parity = cur_parity;
2885 memcpy( out_at[at_1].z_dir, z_dir1, sizeof(out_at[0].z_dir) );
2886 }
2887 if ( !out_at[at_2].parity ||
2888 (next_parity_defined && !ATOM_PARITY_WELL_DEF(abs(out_at[at_2].parity))) ) {
2889 out_at[at_2].parity = next_parity;
2890 memcpy( out_at[at_2].z_dir, z_dir2, sizeof(out_at[0].z_dir) );
2891 }
2892 out_at[at_1].bAmbiguousStereo |= at[at_1].bAmbiguousStereo;
2893 out_at[at_2].bAmbiguousStereo |= at[at_2].bAmbiguousStereo;
2894 num_stored_stereo_bonds ++;
2895 }
2896 }
2897 }
2898
2899 /* stereo + isotopic (all non-zero) */
2900 cur_action = half_stereo_bond_action( cur_parity, nUnknown, 1, vABParityUnknown ); /* -1 => program error */
2901 next_action = half_stereo_bond_action( next_parity, nUnknown, 1, vABParityUnknown );
2902 result_action = inchi_min(cur_action, next_action);
2903 cur_parity = abs_cur_parity;
2904 next_parity = abs_next_parity;
2905 if ( result_action != AB_PARITY_NONE && result_action != -1 ) {
2906 /* stero, isotopic */
2907 if ( cur_parity > 0 && next_parity > 0 ) {
2908 if( save_a_stereo_bond( dot_prod_z, result_action | chain_len_bits,
2909 at_1, i_next_at_1, out_at[at_1].stereo_bond_neighbor2,
2910 out_at[at_1].stereo_bond_ord2, out_at[at_1].stereo_bond_z_prod2,
2911 out_at[at_1].stereo_bond_parity2,
2912 at_2, i_next_at_2, out_at[at_2].stereo_bond_neighbor2,
2913 out_at[at_2].stereo_bond_ord2, out_at[at_2].stereo_bond_z_prod2,
2914 out_at[at_2].stereo_bond_parity2) ) {
2915 if ( !out_at[at_1].parity2 ||
2916 (cur_parity_defined && !ATOM_PARITY_WELL_DEF(abs(out_at[at_1].parity2))) ) {
2917 out_at[at_1].parity2 = cur_parity /*| chain_len_bits*/;
2918 if ( !out_at[at_1].parity ) {
2919 memcpy( out_at[at_1].z_dir, z_dir1, sizeof(out_at[0].z_dir) );
2920 }
2921 }
2922 if ( !out_at[at_2].parity2 || /* next line changed from abs(out_at[at_2].parity) 2006-03-05 */
2923 (next_parity_defined && !ATOM_PARITY_WELL_DEF(abs(out_at[at_2].parity2))) ) {
2924 out_at[at_2].parity2 = next_parity /*| chain_len_bits*/;
2925 if ( !out_at[at_2].parity ) {
2926 memcpy( out_at[at_2].z_dir, z_dir2, sizeof(out_at[0].z_dir) );
2927 }
2928 }
2929 out_at[at_1].bAmbiguousStereo |= at[at_1].bAmbiguousStereo;
2930 out_at[at_2].bAmbiguousStereo |= at[at_2].bAmbiguousStereo;
2931 num_stored_isotopic_stereo_bonds ++;
2932 }
2933 }
2934 } else
2935 if ( result_action == -1 ) {
2936 stop = 1; /* program error? <BRKPT> */
2937 }
2938
2939 }
2940 if ( stop ) {
2941 return CT_CALC_STEREO_ERR;
2942 }
2943 return /*num_stored_stereo_bonds+*/ num_stored_isotopic_stereo_bonds;
2944 }
2945
2946
2947 /*********************************************************************/
2948 /* if isotopic H, D, T added, can the atom be a stereo center? */
2949 #if ( NEW_STEREOCENTER_CHECK == 1 )
2950 /* int bCanInpAtomBeAStereoCenter( inp_ATOM *at, int cur_at ) */
can_be_a_stereo_atom_with_isotopic_H(inp_ATOM * at,int cur_at,int bPointedEdgeStereo)2951 int can_be_a_stereo_atom_with_isotopic_H( inp_ATOM *at, int cur_at, int bPointedEdgeStereo )
2952 {
2953 int nNumNeigh;
2954 if ( (nNumNeigh = bCanInpAtomBeAStereoCenter( at, cur_at, bPointedEdgeStereo )) &&
2955 at[cur_at].valence + at[cur_at].num_H == nNumNeigh &&
2956 at[cur_at].num_H <= NUM_H_ISOTOPES
2957 ) {
2958 return 1;
2959 }
2960 return 0;
2961 }
2962 #else
can_be_a_stereo_atom_with_isotopic_H(inp_ATOM * at,int cur_at)2963 int can_be_a_stereo_atom_with_isotopic_H( inp_ATOM *at, int cur_at )
2964 {
2965 int j, ret = 0;
2966 if ( bCanAtomBeAStereoCenter( at[cur_at].elname, at[cur_at].charge, at[cur_at].radical ) &&
2967 at[cur_at].valence + at[cur_at].num_H == MAX_NUM_STEREO_ATOM_NEIGH &&
2968 at[cur_at].num_H < MAX_NUM_STEREO_ATOM_NEIGH
2969 ) {
2970
2971 for ( j = 0, ret=1; ret && j < at[cur_at].valence; j ++ ) {
2972 if ( (at[cur_at].bond_type[j] & ~BOND_MARK_ALL) != BOND_SINGLE ) {
2973 ret = 0;
2974 }
2975 }
2976 }
2977 return ret;
2978 }
2979 #endif
2980 /***************************************************************/
GetStereocenter0DParity(inp_ATOM * at,int cur_at,int j1,AT_NUMB nSbNeighOrigAtNumb[],int nFlag)2981 int GetStereocenter0DParity( inp_ATOM *at, int cur_at, int j1, AT_NUMB nSbNeighOrigAtNumb[], int nFlag )
2982 {
2983 int parity = AB_PARITY_NONE;
2984 if ( at[cur_at].p_parity && (j1 == MAX_NUM_STEREO_ATOM_NEIGH-1 || j1 == MAX_NUM_STEREO_ATOM_NEIGH) ) {
2985 int i, num_trans_inp, num_trans_neigh;
2986 AT_NUMB nInpNeighOrigAtNumb[MAX_NUM_STEREO_ATOM_NEIGH];
2987 for ( i = 0; i < MAX_NUM_STEREO_ATOM_NEIGH; i ++ ) {
2988 nInpNeighOrigAtNumb[i] = at[cur_at].p_orig_at_num[i];
2989 if ( nInpNeighOrigAtNumb[i] == at[cur_at].orig_at_number ) {
2990 nInpNeighOrigAtNumb[i] = 0; /* lone pair or explicit H */
2991 }
2992 }
2993 num_trans_inp = insertions_sort( nInpNeighOrigAtNumb, MAX_NUM_STEREO_ATOM_NEIGH, sizeof(nInpNeighOrigAtNumb[0]), comp_AT_NUMB );
2994 num_trans_neigh = insertions_sort( nSbNeighOrigAtNumb, j1, sizeof(nSbNeighOrigAtNumb[0]), comp_AT_NUMB );
2995 if ( j1 == MAX_NUM_STEREO_ATOM_NEIGH-1 ) {
2996 ; /*num_trans_neigh += j1;*/ /* the lone pair or implicit H is implicitly at the top of the list */
2997 }
2998 if ( !memcmp( nInpNeighOrigAtNumb + MAX_NUM_STEREO_ATOM_NEIGH-j1, nSbNeighOrigAtNumb, j1*sizeof(AT_NUMB) ) ) {
2999 if ( ATOM_PARITY_WELL_DEF(at[cur_at].p_parity) ) {
3000 parity = 2 - (num_trans_inp + num_trans_neigh + at[cur_at].p_parity) % 2;
3001 } else {
3002 parity = at[cur_at].p_parity;
3003 }
3004 at[cur_at].bUsed0DParity |= nFlag; /* 0D parity used for streocenter parity */
3005 }
3006 }
3007 return parity;
3008 }
3009
3010 /***************************************************************
3011 * Get stereo atom parity for the current order of attachments
3012 * The result in at[cur_at].parity is valid for previously removed
3013 * explicit hydrogen atoms, including isotopic ones, that are located in at_removed_H[]
3014 * The return value is a calculated parity.
3015 */
3016 #define ADD_EXPLICIT_HYDROGEN_NEIGH 1
3017 #define ADD_EXPLICIT_LONE_PAIR_NEIGH 2
set_stereo_atom_parity(sp_ATOM * out_at,inp_ATOM * at,int cur_at,inp_ATOM * at_removed_H,int num_removed_H,int bPointedEdgeStereo,int vABParityUnknown)3018 int set_stereo_atom_parity( sp_ATOM *out_at, inp_ATOM *at, int cur_at, inp_ATOM *at_removed_H,
3019 int num_removed_H, int bPointedEdgeStereo, int vABParityUnknown)
3020 {
3021 int j, k, next_at, num_z, j1, nType, num_explicit_H, tot_num_iso_H, nMustHaveNumNeigh;
3022 int num_explicit_iso_H[NUM_H_ISOTOPES+1]; /* numbers of removed hydrogen atoms */
3023 int index_H[MAX_NUM_STEREO_ATOM_NEIGH]; /* cannot have more than 4 elements: 1 H, 1 D, 1 T atom(s) */
3024 double z, sum_xyz[3], min_sine, triple_product;
3025 double at_coord[MAX_NUM_STEREO_ATOM_NEIGH][3];
3026 double bond_len_xy[4], rmax=0.0, rmin=0.0;
3027 double at_coord_center[3];
3028 int parity, bAmbiguous = 0, bAddExplicitNeighbor = 0, b2D = 0, n2DTetrahedralAmbiguity = 0;
3029 int bIgnoreIsotopicH = (0 != (at[cur_at].cFlags & AT_FLAG_ISO_H_POINT));
3030 AT_NUMB nSbNeighOrigAtNumb[MAX_NUM_STEREO_ATOM_NEIGH];
3031
3032 out_at[cur_at].parity =
3033 out_at[cur_at].parity2 =
3034 out_at[cur_at].stereo_atom_parity =
3035 out_at[cur_at].stereo_atom_parity2 = AB_PARITY_NONE;
3036 parity = AB_PARITY_NONE;
3037
3038 memset(num_explicit_iso_H, 0, sizeof(num_explicit_iso_H));
3039 num_explicit_H = 0;
3040
3041 #if ( NEW_STEREOCENTER_CHECK == 1 )
3042 if ( !(nMustHaveNumNeigh = bCanInpAtomBeAStereoCenter( at, cur_at, bPointedEdgeStereo ) ) ||
3043 at[cur_at].num_H > NUM_H_ISOTOPES
3044 ) {
3045 goto exit_function;
3046 }
3047 #else
3048 nMustHaveNumNeigh = MAX_NUM_STEREO_ATOM_NEIGH;
3049 if ( !bCanAtomBeAStereoCenter( at[cur_at].elname, at[cur_at].charge, at[cur_at].radical ) ||
3050 at[cur_at].valence + at[cur_at].num_H != nMustHaveNumNeigh ||
3051 at[cur_at].num_H > NUM_H_ISOTOPES
3052 ) {
3053 goto exit_function;
3054 }
3055 for ( j = 0; j < at[cur_at].valence; j ++ ) {
3056 if ( (at[cur_at].bond_type[j] & ~BOND_MARK_ALL) != BOND_SINGLE ) {
3057 goto exit_function;
3058 }
3059 }
3060 #endif
3061
3062 /* numbers of isotopic H atoms */
3063 for ( j = 0, tot_num_iso_H = 0; j < NUM_H_ISOTOPES; j ++ ) {
3064 if ( at[cur_at].num_iso_H[j] > 1 ) {
3065 goto exit_function; /* two or more identical hydrogen isotopic neighbors */
3066 }
3067 tot_num_iso_H += at[cur_at].num_iso_H[j];
3068 }
3069 if ( bIgnoreIsotopicH ) {
3070 tot_num_iso_H = 0; /* isotopic H considered subject to exchange => ignore isotopic */
3071 }
3072 /* number of non-isotopic H atoms */
3073 if ( at[cur_at].num_H - tot_num_iso_H > 1 ) {
3074 goto exit_function; /* two or more identical hydrogen non-isotopic neighbors */
3075 }
3076
3077 /* count removed explicit terminal hydrogens attached to at[cur_at]. */
3078 /* the result is num_explicit_H. */
3079 /* Removed hydrogens are sorted in increasing isotopic shift order */
3080 if ( at_removed_H && num_removed_H > 0 ) {
3081 for ( j = 0; j < num_removed_H; j ++ ) {
3082 if ( at_removed_H[j].neighbor[0] == cur_at ) {
3083 k = at_removed_H[j].iso_atw_diff;
3084 /* iso_atw_diff values: H=>0, 1H=>1, D=2H=>2, T=3H=>3 */
3085 if ( k < 0 || k > NUM_H_ISOTOPES || bIgnoreIsotopicH )
3086 k = 0; /* treat wrong H isotopes as non-isotopic H */
3087 num_explicit_iso_H[k] ++;
3088 index_H[num_explicit_H++] = j;
3089 }
3090 }
3091 }
3092
3093 /* coordinates initialization */
3094 num_z = 0;
3095 sum_xyz[0] = sum_xyz[1] = sum_xyz[2] = 0.0;
3096
3097 at_coord_center[0] =
3098 at_coord_center[1] =
3099 at_coord_center[2] = 0.0;
3100
3101 /* fill out stereo center neighbors coordinates */
3102 /* and obtain the parity from the geometry */
3103
3104 for ( k = 0, j1 = 0; k < 2; k ++ ) {
3105 switch( k ) {
3106
3107 case 0:
3108 /* add coordinates of removed hydrogens */
3109 for ( j = 0; j < num_explicit_H; j ++, j1 ++ ) {
3110 next_at = index_H[j];
3111 /* use bond description located at removed_H atom */
3112 /* minus sign at get_z_coord: at_removed_H[] contains bonds TO at[cur_at], not FROM it. */
3113 /* Note: &at[(at_removed_H-at)+ next_at] == &at_removed_H[next_at] */
3114 z = -get_z_coord( at, (at_removed_H-at)+ next_at, 0 /*neighbor #*/, &nType, -(bPointedEdgeStereo & PES_BIT_POINT_EDGE_STEREO) );
3115 switch ( nType ) {
3116 case ZTYPE_EITHER:
3117 parity = vABParityUnknown /*AB_PARITY_UNKN*/ ; /* no parity: bond in "Either" direction. */
3118 goto exit_function;
3119 case ZTYPE_UP:
3120 case ZTYPE_DOWN:
3121 nType = -nType; /* at_removed_H[] contains bonds TO the center, not from */
3122 b2D ++;
3123 /* no break; here */
3124 case ZTYPE_3D:
3125 num_z ++;
3126 }
3127
3128 nSbNeighOrigAtNumb[j1] = at_removed_H[next_at].orig_at_number;
3129 at_coord[j1][0] = at_removed_H[next_at].x-at[cur_at].x;
3130 at_coord[j1][1] = at_removed_H[next_at].y-at[cur_at].y;
3131 bond_len_xy[j1] = len2(at_coord[j1]);
3132 /* bond_len_xy[j1] = sqrt(at_coord[j1][0]*at_coord[j1][0]+at_coord[j1][1]*at_coord[j1][1]); */
3133 at_coord[j1][2] = (nType==ZTYPE_3D? z :
3134 nType==ZTYPE_UP? bond_len_xy[j1] :
3135 nType==ZTYPE_DOWN? -bond_len_xy[j1] : 0.0 );
3136 }
3137 break;
3138 case 1:
3139 /* add all coordinates of other neighboring atoms */
3140 for ( j = 0; j < at[cur_at].valence; j ++, j1 ++ ) {
3141 next_at = at[cur_at].neighbor[j];
3142 z = get_z_coord( at, cur_at, j, &nType, (bPointedEdgeStereo & PES_BIT_POINT_EDGE_STEREO) );
3143 switch ( nType ) {
3144 case ZTYPE_EITHER:
3145 parity = vABParityUnknown /*AB_PARITY_UNKN*/; /* unknown parity: bond in "Either" direction. */
3146 goto exit_function;
3147 case ZTYPE_UP:
3148 case ZTYPE_DOWN:
3149 b2D ++;
3150 case ZTYPE_3D:
3151 num_z ++;
3152 }
3153
3154 nSbNeighOrigAtNumb[j1] = at[next_at].orig_at_number;
3155 at_coord[j1][0] = at[next_at].x-at[cur_at].x;
3156 at_coord[j1][1] = at[next_at].y-at[cur_at].y;
3157 bond_len_xy[j1] = len2(at_coord[j1]);
3158 /* bond_len_xy[j1] = sqrt(at_coord[j1][0]*at_coord[j1][0]+at_coord[j1][1]*at_coord[j1][1]); */
3159 at_coord[j1][2] = (nType==ZTYPE_3D? z :
3160 nType==ZTYPE_UP? bond_len_xy[j1] :
3161 nType==ZTYPE_DOWN? -bond_len_xy[j1] : 0.0 );
3162 }
3163 break;
3164 }
3165 }
3166 /* j1 is the number of explicit neighbors (that is, all neighbors except implicit H) */
3167
3168 b2D = (b2D == num_z && num_z); /* 1 => two-dimensional */
3169
3170 if ( MAX_NUM_STEREO_ATOM_NEIGH != at[cur_at].valence+num_explicit_H &&
3171 MAX_NUM_STEREO_ATOM_NEIGH-1 != at[cur_at].valence+num_explicit_H ) {
3172 /* not enough geometry data to find the central atom parity */
3173 if ( nMustHaveNumNeigh == at[cur_at].valence+at[cur_at].num_H &&
3174 at[cur_at].num_H > 1 ) {
3175 /* only isotopic parity is possible; no non-isotopic parity */
3176 if ( parity == vABParityUnknown /*AB_PARITY_UNKN*/ ) {
3177 parity = -vABParityUnknown /*AB_PARITY_UNKN*/; /* the user marked the center as "unknown" */
3178 } else {
3179 parity = -AB_PARITY_UNDF; /* not enough geometry; only isotopic parity is possible */
3180 }
3181 } else {
3182 parity = AB_PARITY_NONE; /* not a stereocenter at all */
3183 }
3184 goto exit_function;
3185 }
3186 /* make all vector lengths equal to 1; exit if too short. 9-10-2002 */
3187 for ( j = 0; j < j1; j ++ ) {
3188 z = len3( at_coord[j] );
3189 if ( z < MIN_BOND_LEN ) {
3190 /* bond length is too small: use 0D parities */
3191 if ( AB_PARITY_NONE == (parity = GetStereocenter0DParity( at, cur_at, j1, nSbNeighOrigAtNumb, FlagSC_0D )) ) {
3192 parity = AB_PARITY_UNDF;
3193 }
3194 goto exit_function;
3195 }
3196 #if ( STEREO_CENTER_BONDS_NORM == 1 )
3197 else {
3198 mult3( at_coord[j], 1.0/z, at_coord[j] );
3199 }
3200 #endif
3201 rmax = j? inchi_max( rmax, z) : z;
3202 rmin = j? inchi_min( rmin, z) : z;
3203 }
3204 if ( rmin / rmax < MIN_SINE ) {
3205 /* bond ratio is too small: use 0D parities */
3206 if ( AB_PARITY_NONE == (parity = GetStereocenter0DParity( at, cur_at, j1, nSbNeighOrigAtNumb, FlagSC_0D )) ) {
3207 parity = AB_PARITY_UNDF;
3208 }
3209 goto exit_function;
3210 }
3211 for ( j = 0; j < j1; j ++ ) {
3212 add3( sum_xyz, at_coord[j], sum_xyz );
3213 }
3214
3215
3216
3217 /* here j1 is a number of neighbors including explicit terminal isotopic H */
3218 /* num_explicit_iso_H[0] = number of explicit non-isotopic hydrogen atom neighbors */
3219 j = j1;
3220 /* Add Explicit Neighbor */
3221 if ( j1 == MAX_NUM_STEREO_ATOM_NEIGH-1 ) {
3222 /* add an explicit neighbor if possible */
3223 if ( nMustHaveNumNeigh == MAX_NUM_STEREO_ATOM_NEIGH-1 ) {
3224 bAddExplicitNeighbor = ADD_EXPLICIT_LONE_PAIR_NEIGH;
3225 } else
3226 if ( nMustHaveNumNeigh == MAX_NUM_STEREO_ATOM_NEIGH ) {
3227 /* check whether an explicit non-isotopic hydrogen can be added */
3228 /* to an atom that is a stereogenic atom */
3229 if ( 1 == at[cur_at].num_H - num_explicit_H && /* the atom has only one one implicit hydrogen */
3230 1 == at[cur_at].num_H - tot_num_iso_H ) { /* this hydrogen is non-isotopic */
3231 bAddExplicitNeighbor = ADD_EXPLICIT_HYDROGEN_NEIGH;
3232 }
3233 }
3234 }
3235
3236 if ( bAddExplicitNeighbor ) {
3237 /***********************************************************
3238 * May happen only if (j1 == MAX_NUM_STEREO_ATOM_NEIGH-1)
3239 * 3 neighbors only, no H-neighbors. Create and add coordinates of an implicit H
3240 * or a fake 4th neighbor, that is, a lone pair
3241 */
3242 if ( parity == vABParityUnknown /*AB_PARITY_UNKN*/ ) {
3243 goto exit_function; /* the user insists the parity is unknown and the isotopic */
3244 /* composition of the neighbors does not contradict */
3245 } else
3246 if ( num_z == 0 || are_3_vect_in_one_plane(at_coord, MIN_SINE) ) {
3247 /* "hydrogen down" rule is needed to resolve an ambiguity */
3248 if ( num_z > 0 ) {
3249 bAmbiguous |= AMBIGUOUS_STEREO;
3250 }
3251 #if ( APPLY_IMPLICIT_H_DOWN_RULE == 1 ) /* { */
3252 /* Although H should be at the top of the list, add it to the bottom. */
3253 /* This will be taken care of later by inverting parity 1<->2 */
3254 at_coord[j][0] = 0.0;
3255 at_coord[j][1] = 0.0;
3256 #if ( STEREO_CENTER_BONDS_NORM == 1 )
3257 at_coord[j][2] = -1.0;
3258 #else
3259 at_coord[j][2] = -(bond_len_xy[0]+bond_len_xy[1]+bond_len_xy[2])/3.0;
3260 #endif
3261 #else /* } APPLY_IMPLICIT_H_DOWN_RULE { */
3262 #if (ALWAYS_SET_STEREO_PARITY == 1)
3263 parity = AB_PARITY_EVEN; /* suppose atoms are pre-sorted (testing) */
3264 #else
3265 /* all 3 bonds are in one plain: try to get 0D parities */
3266 if ( AB_PARITY_NONE == (parity = GetStereocenter0DParity( at, cur_at, j1, nSbNeighOrigAtNumb, FlagSC_0D )) ) {
3267 parity = AB_PARITY_UNDF;
3268 }
3269 /*parity = AB_PARITY_UNDF;*/ /* no parity can be calculated found */
3270 #endif
3271 goto exit_function;
3272 #endif /* } APPLY_IMPLICIT_H_DOWN_RULE */
3273 } else {
3274 /* we have enough information to find implicit hydrogen coordinates */
3275 /*
3276 at_coord[j][0] = -sum_x;
3277 at_coord[j][1] = -sum_y;
3278 at_coord[j][2] = -sum_z;
3279 */
3280 copy3( sum_xyz, at_coord[j] );
3281 change_sign3( at_coord[j], at_coord[j] );
3282 z = len3( at_coord[j] );
3283 #if ( FIX_STEREO_SCALING_BUG == 1 )
3284 if ( z > 1.0 ) {
3285 rmax *= z;
3286 } else {
3287 rmin *= z;
3288 }
3289 #else
3290 /* Comparing the original bond lengths to lenghts derived from normalized to 1 */
3291 /* This bug leads to pronouncing legitimate stereogenic atoms */
3292 /* connected by 3 bonds "undefined" if in a nicely drawn 2D structure */
3293 /* bond lengths are about 20 or greater. Reported by Reinhard Dunkel 2005-08-05 */
3294 if ( bPointedEdgeStereo & PES_BIT_FIX_SP3_BUG ) {
3295 /* coordinate scaling bug fixed here */
3296 if ( z > 1.0 ) {
3297 rmax *= z;
3298 } else {
3299 rmin *= z;
3300 }
3301 } else {
3302 /* original InChI v.1 bug */
3303 rmax = inchi_max( rmax, z );
3304 rmin = inchi_min( rmin, z );
3305 }
3306 #endif
3307 if ( z < MIN_BOND_LEN || rmin/rmax < MIN_SINE ) {
3308 /* the new 4th bond is too short: try to get 0D parities */
3309 if ( AB_PARITY_NONE == (parity = GetStereocenter0DParity( at, cur_at, j1, nSbNeighOrigAtNumb, FlagSC_0D )) ) {
3310 parity = AB_PARITY_UNDF;
3311 }
3312 goto exit_function;
3313 }
3314 #if ( STEREO_CENTER_BOND4_NORM == 1 )
3315 else {
3316 mult3( at_coord[j], 1.0/z, at_coord[j] );
3317 }
3318 #endif
3319 }
3320 } else
3321 if ( j1 != MAX_NUM_STEREO_ATOM_NEIGH ) {
3322 if ( parity == vABParityUnknown /*AB_PARITY_UNKN*/ ) {
3323 parity = -AB_PARITY_UNDF; /* isotopic composition of H-neighbors contradicts 'unknown' */
3324 }
3325 goto exit_function;
3326 } else /* j1 == MAX_NUM_STEREO_ATOM_NEIGH */
3327 if ( num_z == 0 || are_4at_in_one_plane(at_coord, MIN_SINE) ) {
3328 /* all four neighours in xy plane: undefined geometry. */
3329 if ( num_z > 0 ) {
3330 bAmbiguous |= AMBIGUOUS_STEREO;
3331 }
3332 if ( parity != vABParityUnknown /*AB_PARITY_UNKN*/ ) {
3333 #if (ALWAYS_SET_STEREO_PARITY == 1)
3334 parity = AB_PARITY_EVEN; /* suppose atoms are pre-sorted (testing) */
3335 #else
3336 /* all 4 bonds are in one plain: try to get 0D parities */
3337 if ( AB_PARITY_NONE == (parity = GetStereocenter0DParity( at, cur_at, j1, nSbNeighOrigAtNumb, FlagSC_0D )) ) {
3338 parity = AB_PARITY_UNDF;
3339 } else
3340 if ( ATOM_PARITY_WELL_DEF( parity ) ) {
3341 bAmbiguous &= ~AMBIGUOUS_STEREO; /* 0D parity has resolved the ambiguity */
3342 }
3343 #endif
3344 }
3345 goto exit_function;
3346 }
3347 /***********************************************************
3348 * At this point we have 4 neighboring atoms.
3349 * check for tetrahedral ambiguity in 2D case
3350 */
3351 if ( b2D )
3352 {
3353
3354 n2DTetrahedralAmbiguity = Get2DTetrahedralAmbiguity( at_coord, bAddExplicitNeighbor, (bPointedEdgeStereo & PES_BIT_FIX_SP3_BUG ) );
3355
3356 if ( 0 < n2DTetrahedralAmbiguity )
3357 {
3358 if ( T2D_WARN & n2DTetrahedralAmbiguity ) {
3359 bAmbiguous |= AMBIGUOUS_STEREO;
3360 }
3361 if ( T2D_UNDF & n2DTetrahedralAmbiguity ) {
3362 if ( parity != vABParityUnknown /*AB_PARITY_UNKN*/ ) {
3363 #if (ALWAYS_SET_STEREO_PARITY == 1)
3364 parity = AB_PARITY_EVEN; /* suppose atoms are pre-sorted (testing) */
3365 #else
3366 parity = AB_PARITY_UNDF; /* no parity */
3367 #endif
3368 }
3369 goto exit_function;
3370 }
3371 } else
3372 if ( n2DTetrahedralAmbiguity < 0 ) {
3373 bAmbiguous |= AMBIGUOUS_STEREO_ERROR; /* error */
3374 parity = AB_PARITY_UNDF;
3375 goto exit_function;
3376 }
3377 }
3378
3379 /************************************************************/
3380 /* Move coordinates origin to the neighbor #0 */
3381 for ( j = 1; j < MAX_NUM_STEREO_ATOM_NEIGH; j ++ ) {
3382 diff3(at_coord[j], at_coord[0], at_coord[j]);
3383 }
3384 diff3(at_coord_center, at_coord[0], at_coord_center);
3385
3386 /*
3387 for ( k = 0; k < 3; k++ ) {
3388 for ( j = 1; j < MAX_NUM_STEREO_ATOM_NEIGH; j ++ ) {
3389 at_coord[j][k] -= at_coord[0][k];
3390 }
3391 at_coord_center[k] -= at_coord[0][k];
3392 }
3393 */
3394 /********************************************************
3395 * find the central (cur_at) atom's parity
3396 * (orientation of atoms #1-3 when looking from #0)
3397 ********************************************************/
3398 triple_product = triple_prod_and_min_abs_sine2(&at_coord[1], at_coord_center, bAddExplicitNeighbor, &min_sine, &bAmbiguous);
3399 /*
3400 * check for tetrahedral ambiguity -- leave it out for now
3401 */
3402 if ( fabs(triple_product) > ZERO_FLOAT && (min_sine > MIN_SINE || (fabs(min_sine) > ZERO_FLOAT && (n2DTetrahedralAmbiguity & T2D_OKAY )) ) ) {
3403 /* Even => sorted in correct order, Odd=>transposed */
3404 parity = triple_product > 0.0? AB_PARITY_EVEN : AB_PARITY_ODD;
3405 /* if ( num_explicit_H && at[cur_at].removed_H_parity % 2 ) */
3406 /* odd transposition of the removed implicit H */
3407 /* out_at[cur_at].parity = 3 - out_at[cur_at].parity; */
3408
3409 /* moved; see below */
3410 /* out_at[cur_at].bAmbiguousStereo |= bAmbiguous; */
3411 /* at[cur_at].bAmbiguousStereo |= bAmbiguous; */
3412
3413 /* for 4 attached atoms, moving the implicit H from index=3 to index=0 */
3414 /* can be done in odd number (3) transpositions: (23)(12)(01), which inverts the parity */
3415 if ( j1 == MAX_NUM_STEREO_ATOM_NEIGH-1 ) {
3416 parity = 3 - parity;
3417 }
3418 } else {
3419 #if (ALWAYS_SET_STEREO_PARITY == 1)
3420 parity = AT_PARITY_EVEN; /* suppose atoms are pre-sorted (testing) */
3421 #else
3422 if ( num_z > 0 ) {
3423 bAmbiguous |= AMBIGUOUS_STEREO;
3424 }
3425 parity = AB_PARITY_UNDF; /* no parity: 4 bonds are in one plane. */
3426 #endif
3427 }
3428 exit_function:
3429
3430 if ( parity ) {
3431 out_at[cur_at].bAmbiguousStereo |= bAmbiguous;
3432 at[cur_at].bAmbiguousStereo |= bAmbiguous;
3433 }
3434
3435
3436 /* non-isotopic parity */
3437 if ( at[cur_at].num_H > 1 || parity <= 0 )
3438 ; /* no non-isotopic parity */
3439 else
3440 out_at[cur_at].parity = parity;
3441
3442 /* isotopic parity */
3443 if ( parity == -AB_PARITY_UNDF || parity == -vABParityUnknown /*AB_PARITY_UNKN*/ )
3444 parity = -parity;
3445 if ( parity < 0 )
3446 parity = AB_PARITY_NONE;
3447 out_at[cur_at].parity2 = parity;
3448
3449
3450 parity = PARITY_VAL(out_at[cur_at].parity);
3451 out_at[cur_at].stereo_atom_parity = ATOM_PARITY_WELL_DEF( parity )? AB_PARITY_CALC : parity;
3452 parity = PARITY_VAL(out_at[cur_at].parity2);
3453 out_at[cur_at].stereo_atom_parity2 = ATOM_PARITY_WELL_DEF( parity )? AB_PARITY_CALC : parity;
3454 /*
3455 out_at[cur_at].parity2 = out_at[cur_at].parity; // save for stereo + isotopic canon.
3456 if ( out_at[cur_at].parity ) {
3457 if ( num_explicit_H > 1 || j1 == MAX_NUM_STEREO_ATOM_NEIGH-1 && num_explicit_H ) {
3458 // X H X
3459 // for example, >C< or >C-D
3460 // Y D Y
3461 // parity exists for stereo + isotopic atoms canonicalization only
3462 out_at[cur_at].parity = 0;
3463 }
3464 }
3465 // returning 0 means this can be an adjacent to a stereogenic bond atom
3466 */
3467
3468 return (int)out_at[cur_at].parity2;
3469
3470 }
3471 #undef ADD_EXPLICIT_HYDROGEN_NEIGH
3472 #undef ADD_EXPLICIT_LONE_PAIR_NEIGH
3473
3474 /*************************************************************/
set_stereo_parity(inp_ATOM * at,sp_ATOM * at_output,int num_at,int num_removed_H,int * nMaxNumStereoAtoms,int * nMaxNumStereoBonds,INCHI_MODE nMode,int bPointedEdgeStereo,int vABParityUnknown)3475 int set_stereo_parity( inp_ATOM* at, sp_ATOM* at_output, int num_at, int num_removed_H,
3476 int *nMaxNumStereoAtoms, int *nMaxNumStereoBonds, INCHI_MODE nMode,
3477 int bPointedEdgeStereo, int vABParityUnknown )
3478 {
3479 int num_3D_stereo_atoms=0;
3480 int num_stereo_bonds=0; /* added to fix allene stereo bug reported for FClC=C=CFCl by Burt Leland - 2009-02-05 DT */
3481
3482 int i, is_stereo, num_stereo, max_stereo_atoms=0, max_stereo_bonds=0;
3483 QUEUE *q = NULL;
3484 AT_RANK *nAtomLevel = NULL;
3485 S_CHAR *cSource = NULL;
3486 AT_RANK min_sb_ring_size = 0;
3487
3488 /**********************************************************
3489 *
3490 * Note: this parity reflects only relative positions of
3491 * the atoms-neighbors and their ordering in the
3492 * lists of neighbors.
3493 *
3494 * To obtain the actual parity, the parity of a number
3495 * of neighbors transpositions (to obtain a sorted
3496 * list of numbers assigned to the atoms) should be
3497 * added.
3498 *
3499 **********************************************************/
3500
3501 /*********************************************************************************
3502
3503 An example of parity=1 for stereogenic center, tetrahedral asymmetric atom
3504
3505
3506
3507 (1)
3508 |
3509 |
3510 [C] |
3511 |
3512 (2)------(0)
3513 /
3514 /
3515 /
3516 /
3517 (3)
3518
3519
3520 Notation: (n) is a tetrahedral atom neighbor; n is an index of a neighbor in
3521 the central_at->neighbor[] array : neighbor atom number is central_at->neighbor[n].
3522
3523 (0)-(1), (0)-(2), (0)-(3) are lines connecting atom [C] neighbors to neighbor (0)
3524 (0), (1) and (2) are in the plane
3525 (0)-(3) is directed from the plain to the viewer
3526 [C] is somewhere between (0), (1), (2), (3)
3527 Since (1)-(2)-(3) are in a clockwise order when looking from (0), parity is 2, or even;
3528 otherwise parity would be 1, or odd.
3529
3530 **********************************************************************************
3531
3532 Examples of a stereogenic bond.
3533
3534 Notation: [atom number], (index of a neighbor):
3535 [1] and [2] are atoms connected by the stereogenic bond
3536 numbers in () are indexes of neighbors of [1] or [2].
3537 (12 x 16)z = z-component of [1]-[2] and [1]-[6] cross-product
3538
3539 atom [1] atom [2]
3540 [8] [4] prod01 = (12 x 16)z < 0 prod01 = (21 x 24)z < 0
3541 \ / prod02 = (12 x 18)z > 0 prod02 = (21 x 25)z > 0
3542 (2) (1) 0 transpositions because 0 transpositions because
3543 \ / double bond is in 0 posit. double bond is in 0 position
3544 [1]==(0)(0)==[2] 0 = (prod01 > prod02) 0 = (prod01 > prod02)
3545 / \
3546 (1) (2) result: parity = 2, even result: parity=2, even
3547 / \
3548 [6] [5]
3549
3550
3551
3552 atom [1] atom [2]
3553 [8] [5] prod01 = (12 x 18)z > 0 prod01 = (21 x 24)z > 0
3554 \ / prod02 = (12 x 16)z < 0 prod02 = (21 x 25)z < 0
3555 (0) (2) 2 transpositions to move 1 transposition to move
3556 \ / at [2] from 2 to 0 pos. at [1] from 1 to 0 position
3557 [1]==(2)(1)==[2] 1 = (prod01 > prod02) 1 = (prod01 > prod02)
3558 / \
3559 (1) (0) result: parity = (1+2) result: parity=(1+1)
3560 / \ 2-(1+2)%2 = 1, odd 2-(1+1)%2 = 2, even
3561 [6] [4]
3562
3563
3564 ***********************************************************************************
3565 Note: atoms' numbers [1], [2], [4],... are not used to calculate parity at this
3566 point. They will be used for each numbering in the canonicalization.
3567 Note: parity=3 for a stereo atom means entered undefined bond direction
3568 parity=4 for an atom means parity cannot be determined from the given geometry
3569 ***********************************************************************************/
3570
3571 if ( !at_output || !at ) {
3572 return -1;
3573 }
3574
3575 /* clear stereo descriptors */
3576
3577 for( i = 0; i < num_at; i ++ ) {
3578 at_output[i].parity = 0;
3579 at_output[i].parity2 = 0;
3580 memset(&at_output[i].stereo_bond_neighbor[0], 0, sizeof(at_output[0].stereo_bond_neighbor) );
3581 memset(&at_output[i].stereo_bond_neighbor2[0], 0, sizeof(at_output[0].stereo_bond_neighbor2) );
3582 memset(&at_output[i].stereo_bond_ord[0], 0, sizeof(at_output[0].stereo_bond_ord) );
3583 memset(&at_output[i].stereo_bond_ord2[0], 0, sizeof(at_output[0].stereo_bond_ord2) );
3584 memset(&at_output[i].stereo_bond_z_prod[0], 0, sizeof(at_output[0].stereo_bond_z_prod) );
3585 memset(&at_output[i].stereo_bond_z_prod2[0], 0, sizeof(at_output[0].stereo_bond_z_prod2) );
3586 memset(&at_output[i].stereo_bond_parity[0], 0, sizeof(at_output[0].stereo_bond_parity) );
3587 memset(&at_output[i].stereo_bond_parity2[0], 0, sizeof(at_output[0].stereo_bond_parity2) );
3588 }
3589 /* estimate max numbers of stereo atoms and bonds if isotopic H are added */
3590 if ( nMaxNumStereoAtoms || nMaxNumStereoBonds ) {
3591 for( i = 0, num_stereo = 0; i < num_at; i ++ ) {
3592 int num;
3593 num = can_be_a_stereo_atom_with_isotopic_H( at, i, bPointedEdgeStereo );
3594 if ( num ) {
3595 max_stereo_atoms += num;
3596 } else
3597 if ( (num = can_be_a_stereo_bond_with_isotopic_H( at, i, nMode ) ) ) { /* accept cumulenes */
3598 max_stereo_bonds += num;
3599 }
3600 }
3601 if ( nMaxNumStereoAtoms )
3602 *nMaxNumStereoAtoms = max_stereo_atoms;
3603 if ( nMaxNumStereoBonds )
3604 *nMaxNumStereoBonds = max_stereo_bonds;
3605 }
3606 /* calculate stereo descriptors */
3607 #if ( MIN_SB_RING_SIZE > 0 )
3608 min_sb_ring_size = (AT_RANK)(((nMode & REQ_MODE_MIN_SB_RING_MASK) >> REQ_MODE_MIN_SB_RING_SHFT) & AT_RANK_MASK);
3609 if ( min_sb_ring_size >= 3 ) {
3610 /* create BFS data structure for finding for each stereo bond its min. ring sizes */
3611 q = QueueCreate( num_at+1, sizeof(qInt) );
3612 nAtomLevel = (AT_RANK*)inchi_calloc(sizeof(nAtomLevel[0]),num_at);
3613 cSource = (S_CHAR *)inchi_calloc(sizeof(cSource[0]),num_at);
3614 if ( !q || !cSource || !nAtomLevel ) {
3615 num_3D_stereo_atoms = CT_OUT_OF_RAM;
3616 goto exit_function;
3617 }
3618 } else {
3619 min_sb_ring_size = 2;
3620 }
3621 #endif
3622 /* main cycle: set stereo parities */
3623 for( i = 0, num_stereo = 0; i < num_at; i ++ )
3624 {
3625 is_stereo = set_stereo_atom_parity( at_output, at, i, at+num_at, num_removed_H,
3626 bPointedEdgeStereo, vABParityUnknown ) ;
3627 if ( is_stereo )
3628 {
3629 num_3D_stereo_atoms += ATOM_PARITY_WELL_DEF( is_stereo );
3630 }
3631 else
3632 {
3633 is_stereo = set_stereo_bonds_parity( at_output, at, i, at+num_at, num_removed_H, nMode,
3634 q, nAtomLevel, cSource, min_sb_ring_size,
3635 bPointedEdgeStereo, vABParityUnknown );
3636 if ( RETURNED_ERROR( is_stereo ) ) {
3637 num_3D_stereo_atoms = is_stereo;
3638 break;
3639 }
3640 num_stereo_bonds += (is_stereo != 0); /* added to fix bug reported by Burt Leland - 2009-02-05 DT */
3641 }
3642 num_stereo += (is_stereo != 0);
3643 is_stereo = is_stereo;
3644 }
3645
3646 /* added to fix bug reported by Burt Leland - 2009-02-05 DT */
3647 if ( max_stereo_atoms < num_3D_stereo_atoms && nMaxNumStereoAtoms )
3648 *nMaxNumStereoAtoms = num_3D_stereo_atoms;
3649 if ( max_stereo_bonds < num_stereo_bonds && nMaxNumStereoBonds )
3650 *nMaxNumStereoBonds = num_stereo_bonds;
3651
3652 /*
3653 if ( (nMode & REQ_MODE_SC_IGN_ALL_UU )
3654 REQ_MODE_SC_IGN_ALL_UU
3655 REQ_MODE_SB_IGN_ALL_UU
3656 */
3657
3658 #if ( MIN_SB_RING_SIZE > 0 )
3659 if ( q ) {
3660 q = QueueDelete( q );
3661 }
3662 if ( nAtomLevel )
3663 inchi_free( nAtomLevel );
3664 if ( cSource )
3665 inchi_free( cSource );
3666 exit_function:
3667 #endif
3668
3669
3670 return num_3D_stereo_atoms;
3671 }
3672 /*****************************************************************
3673 * Functions that disconnect bonds
3674 *
3675 *=== During Preprocessing ===
3676 *
3677 * RemoveInpAtBond
3678 * DisconnectMetalSalt (is not aware of bond parities)
3679 * DisconnectAmmoniumSalt
3680 *
3681 *=== Before Normalization ===
3682 *
3683 * remove_terminal_HDT
3684 *
3685 *=== During the Normalization ===
3686 *
3687 * AddOrRemoveExplOrImplH
3688 *
3689 *****************************************************************/
ReconcileAllCmlBondParities(inp_ATOM * at,int num_atoms,int bDisconnected)3690 int ReconcileAllCmlBondParities( inp_ATOM *at, int num_atoms, int bDisconnected )
3691 {
3692 int i, ret = 0;
3693 S_CHAR *visited = (S_CHAR*) inchi_calloc( num_atoms, sizeof(*visited) );
3694 if ( !visited )
3695 return -1; /* out of RAM */
3696 for ( i = 0; i < num_atoms; i ++ ) {
3697 if ( at[i].sb_parity[0] && !visited[i] && !(bDisconnected && is_el_a_metal(at[i].el_number)) ) {
3698 if ( (ret = ReconcileCmlIncidentBondParities( at, i, -1, visited, bDisconnected )) ) {
3699 break; /* error */
3700 }
3701 }
3702 }
3703 inchi_free ( visited );
3704 return ret;
3705 }
3706 /*****************************************************************/
ReconcileCmlIncidentBondParities(inp_ATOM * at,int cur_atom,int prev_atom,S_CHAR * visited,int bDisconnected)3707 int ReconcileCmlIncidentBondParities( inp_ATOM *at, int cur_atom, int prev_atom, S_CHAR *visited, int bDisconnected )
3708 {
3709 /* visited = 0 or parity => atom has not been visited
3710 10 + parity => currently is on the stack + its final parity
3711 20 + parity => has been visited; is not on the stack anymore + its final parity */
3712 int i, j, nxt_atom, ret = 0, len;
3713 int icur2nxt, icur2neigh; /* cur atom neighbors */
3714 int inxt2cur, inxt2neigh; /* next atom neighbors */
3715 int cur_parity, nxt_parity;
3716 int cur_order_parity, nxt_order_parity, cur_sb_parity, nxt_sb_parity, bCurMask, bNxtMask;
3717 /* !(bDisconnected && is_el_a_metal(at[i].el_number) */
3718
3719 if ( at[cur_atom].valence > MAX_NUM_STEREO_BONDS )
3720 return 0; /* ignore */
3721
3722 if ( !at[cur_atom].sb_parity[0] )
3723 return 1; /* wrong call */
3724
3725 if ( visited[cur_atom] >= 10 )
3726 return 2; /* program error */
3727
3728 cur_parity = visited[cur_atom] % 10;
3729
3730 visited[cur_atom] += 10;
3731
3732 for ( i = 0; i < MAX_NUM_STEREO_BONDS && at[cur_atom].sb_parity[i]; i ++ ) {
3733 icur2nxt = (int)at[cur_atom].sb_ord[i];
3734 len = get_opposite_sb_atom( at, cur_atom, icur2nxt, &nxt_atom, &inxt2cur, &j );
3735 if ( !len ) {
3736 return 4; /* could not find the opposite atom: bond parity data error */
3737 }
3738 if ( nxt_atom == prev_atom )
3739 continue;
3740 if ( visited[nxt_atom] >= 20 )
3741 continue; /* back edge, second visit: ignore */
3742 if ( at[nxt_atom].valence > MAX_NUM_STEREO_BONDS )
3743 continue; /* may be treated only after metal disconnection */
3744
3745 if ( bDisconnected && (at[cur_atom].sb_parity[i] & SB_PARITY_FLAG) ) {
3746 cur_sb_parity = (at[cur_atom].sb_parity[i] >> SB_PARITY_SHFT);
3747 bCurMask = 3 << SB_PARITY_SHFT;
3748 } else {
3749 cur_sb_parity = (at[cur_atom].sb_parity[i] & SB_PARITY_MASK);
3750 bCurMask = 3;
3751 }
3752 if ( bDisconnected && (at[nxt_atom].sb_parity[j] & SB_PARITY_FLAG) ) {
3753 nxt_sb_parity = (at[nxt_atom].sb_parity[j] >> SB_PARITY_SHFT);
3754 bNxtMask = 3 << SB_PARITY_SHFT;
3755 } else {
3756 nxt_sb_parity = (at[nxt_atom].sb_parity[j] & SB_PARITY_MASK);
3757 bNxtMask = 3;
3758 }
3759
3760
3761
3762 if ( !ATOM_PARITY_WELL_DEF(cur_sb_parity) ||
3763 !ATOM_PARITY_WELL_DEF(nxt_sb_parity) ) {
3764 if ( cur_sb_parity == nxt_sb_parity ) {
3765 continue;
3766 /*goto move_forward;*/ /* bypass unknown/undefined */
3767 }
3768 return 3; /* sb parities do not match: bond parity data error */
3769 }
3770
3771 icur2neigh = (int)at[cur_atom].sn_ord[i];
3772 inxt2neigh = (int)at[nxt_atom].sn_ord[j];
3773 /* parity of at[cur_atom].neighbor[] premutation to reach this order: { next_atom, neigh_atom, ...} */
3774
3775 /* 1. move next_atom from position=icur2nxt to position=0 =>
3776 * icur2nxt permutations
3777 * 2. move neigh_atom from position=inxt2neigh+(inxt2cur > inxt2neigh) to position=1 =>
3778 * inxt2neigh+(inxt2cur > inxt2neigh)-1 permutations.
3779 * Note if (inxt2cur > inxt2neigh) then move #1 increments neigh_atom position
3780 * Note add 4 because icur2neigh may be negative due to isotopic H removal
3781 */
3782 cur_order_parity = (4+icur2nxt + icur2neigh + (icur2neigh > icur2nxt)) % 2;
3783 /* same for next atom: */
3784 /* parity of at[nxt_atom].neighbor[] premutation to reach this order: { cur_atom, neigh_atom, ...} */
3785 nxt_order_parity = (4+inxt2cur + inxt2neigh + (inxt2neigh > inxt2cur)) % 2;
3786
3787 nxt_parity = visited[nxt_atom] % 10;
3788
3789 if ( !cur_parity ) {
3790 cur_parity = 2 - (cur_order_parity + cur_sb_parity) % 2;
3791 visited[cur_atom] += cur_parity;
3792 } else
3793 if ( cur_parity != 2 - (cur_order_parity + cur_sb_parity) % 2 ) {
3794
3795 /***** reconcile bond parities *****/
3796
3797 /* Each bond parity is split into two values located at the end atoms.
3798 For T (trans) the values are (1,1) or (2,2)
3799 For C (cis) the values are (1,2) or (2,1)
3800 The fact that one pair = another with inverted parities, namely
3801 Inv(1,1) = (2,2) and Inv(1,2) = (2,1), allows one to
3802 simultaneouly invert parities of the current bond end atoms
3803 (at[cur_atom].sb_parity[i], at[nxt_atom].sb_parity[j])
3804 so that the final current atom parity cur_parity
3805 calculated later in stereochemical canonicalization for
3806 each stereobond incident with the current atomis same.
3807 Achieving this is called here RECONCILIATION.
3808 If at the closure of an aromatic circuit the parities of
3809 next atom cannot be reconciled with already calculated then
3810 this function returns 5 (error).
3811 */
3812
3813 at[cur_atom].sb_parity[i] ^= bCurMask;
3814 at[nxt_atom].sb_parity[j] ^= bNxtMask;
3815 cur_sb_parity ^= 3;
3816 nxt_sb_parity ^= 3;
3817 }
3818
3819 if ( !nxt_parity ) {
3820 nxt_parity = 2 - (nxt_order_parity + nxt_sb_parity) % 2;
3821 visited[nxt_atom] += nxt_parity;
3822 } else
3823 if ( nxt_parity != 2 - (nxt_order_parity + nxt_sb_parity) % 2 ) {
3824 return 5; /* algorithm does not work for Mebius-like structures */
3825 }
3826
3827 /* move_forward: */
3828 if ( visited[nxt_atom] < 10 ) {
3829 ret = ReconcileCmlIncidentBondParities( at, nxt_atom, cur_atom, visited, bDisconnected );
3830 if ( ret ) {
3831 break;
3832 }
3833 }
3834 }
3835 visited[cur_atom] += 10; /* all bonds incident to the current atom have
3836 been processed or an error occurred. */
3837 return ret;
3838 }
3839 /*****************************************************************/
get_opposite_sb_atom(inp_ATOM * at,int cur_atom,int icur2nxt,int * pnxt_atom,int * pinxt2cur,int * pinxt_sb_parity_ord)3840 int get_opposite_sb_atom( inp_ATOM *at, int cur_atom, int icur2nxt, int *pnxt_atom, int *pinxt2cur, int *pinxt_sb_parity_ord )
3841 {
3842 AT_NUMB nxt_atom;
3843 int j, len;
3844
3845 len = 0;
3846 while ( len ++ < 20 ) { /* arbitrarily set cumulene length limit to avoid infinite loop */
3847 nxt_atom = at[cur_atom].neighbor[icur2nxt];
3848 for ( j = 0; j < MAX_NUM_STEREO_BONDS && at[nxt_atom].sb_parity[j]; j ++ ) {
3849 if ( cur_atom == at[nxt_atom].neighbor[(int)at[nxt_atom].sb_ord[j]] ) {
3850 /* found the opposite atom */
3851 *pnxt_atom = nxt_atom;
3852 *pinxt2cur = at[nxt_atom].sb_ord[j];
3853 *pinxt_sb_parity_ord = j;
3854 return len;
3855 }
3856 }
3857 if ( j ) {
3858 return 0; /* reached atom(s) with stereobond (sb) parity, the opposite atom has not been found */
3859 }
3860 if ( at[nxt_atom].valence == 2 && 2*BOND_TYPE_DOUBLE == at[nxt_atom].chem_bonds_valence ) {
3861 /* follow cumulene =X= path */
3862 icur2nxt = (at[nxt_atom].neighbor[0] == cur_atom);
3863 cur_atom = nxt_atom;
3864 } else {
3865 return 0; /* neither atom with a sb parity not middle cumulene could be reached */
3866 }
3867 }
3868 return 0; /* too long chain of cumulene was found */
3869 }
3870