ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/fnd_ring.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 19129 byte(s)
Log Message:
test

Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5 #include "utility.h"
6 #include "gmmx.h"
7
8 #define MAX_RINGSIZE 30
9
10 void message_alert(char *,char *);
11 void find_rings(void);
12 void fnd_rings1(int *, int **);
13 void del_atom(int, int **);
14 void log_ring1(int, int *);
15 int icompare(int, int *, int *);
16 void gmmxsort(int, int *);
17 int in_ring(int,int);
18 void add_all_rings(int,int *);
19 void trim(int ,int *,int *);
20 void check_N2(int,int *,int *);
21 void getUnion(int *, int *,int *);
22 int getRing(int,int,int *, int **,int *);
23 void delete_bond(int,int, int **);
24 void breakBond(int,int **);
25 void repack(int,int **);
26 void addq(int,int *,int);
27 int deleteq(int *,int);
28 void add_pair(int,int,int,int);
29 void find_rings1(int *,int **);
30 int get_intersection(int *, int *);
31
32 void find_rings()
33 {
34 int i,j,k,l,kk, jji;
35 int found_terminal, ntemp;
36 int ia,ib,icount,found,tmp_size, tmp_ring[50];
37 int *iexist, **itemiat;
38
39
40 ntemp = natom;
41 iexist = ivector(1, ntemp+1);
42 itemiat = imatrix(1,ntemp+1, 0,MAXIAT);
43 gmmxring.nrings = 0;
44 for (i=0; i < 200; i++)
45 {
46 gmmxring.nring_atoms[i] = 0;
47 gmmxring.fused[i] = -1;
48 for (j=0; j < 30; j++)
49 gmmxring.ring_atoms[i][j] = 0;
50 }
51
52 for (i=1; i <= natom; i++)
53 {
54 iexist[i] = atom[i].type;
55 for(j=0; j < MAXIAT; j++)
56 itemiat[i][j] = atom[i].iat[j];
57 }
58
59 L_1:
60 found_terminal = FALSE;
61 for (i=1; i <= natom; i++)
62 {
63 if (iexist[i] != 0)
64 {
65 jji = 0;
66 for(j=0; j < MAXIAT; j++)
67 {
68 if (itemiat[i][j] != 0)
69 jji++;
70 }
71 if (jji <= 1)
72 {
73 del_atom(i,itemiat);
74 iexist[i] = 0;
75 found_terminal = TRUE;
76 }
77 }
78 }
79 if (found_terminal == TRUE)
80 goto L_1;
81 find_rings1(iexist, itemiat);
82 free_ivector(iexist, 1, ntemp+1);
83 free_imatrix(itemiat, 1, ntemp+1, 0, MAXIAT);
84 // sort on ring size
85 L_2:
86 found = FALSE;
87 for (i=0; i < gmmxring.nrings-1; i++)
88 {
89 for (j=i+1; j < gmmxring.nrings; j++)
90 {
91 if (gmmxring.nring_atoms[j] < gmmxring.nring_atoms[i])
92 {
93 tmp_size = gmmxring.nring_atoms[i];
94 for (k=0; k < tmp_size; k++)
95 tmp_ring[k] = gmmxring.ring_atoms[i][k];
96
97 gmmxring.nring_atoms[i] = gmmxring.nring_atoms[j];
98 for (k=0; k < gmmxring.nring_atoms[j]; k++)
99 gmmxring.ring_atoms[i][k] = gmmxring.ring_atoms[j][k];
100
101 gmmxring.nring_atoms[j] = tmp_size;
102 for (k=0; k < tmp_size; k++)
103 gmmxring.ring_atoms[j][k] = tmp_ring[k];
104 found = TRUE;
105 }
106 }
107 }
108 if (found == TRUE) goto L_2;
109 // find fused rings
110 // now find branch atom pairs and what rings they occur in
111 gmmxring.npairs = 0;
112 for (i=0; i < 200; i++)
113 {
114 gmmxring.pair[i][0] = -1;
115 gmmxring.pair[i][1] = -1;
116 for (j=0; j < 10; j++)
117 gmmxring.in_rings[i][j] = -1;
118 }
119 for (i=0; i < gmmxring.nbranch; i++)
120 {
121 icount = 0;
122 ia = gmmxring.branch_atoms[i];
123 for (j=0; j < gmmxring.nrings; j++)
124 {
125 if (in_ring(ia,j))
126 {
127 gmmxring.in_rings[i][icount] = j;
128 icount++;
129 }
130 }
131 // if icount == 1 we should remove - not a branch atom
132 }
133 // found which rings each branch atom belongs to
134 for (i=0; i < gmmxring.nbranch; i++)
135 {
136 for (k=0; gmmxring.in_rings[i][k] != -1; k++)
137 {
138 ia = gmmxring.in_rings[i][k];
139 for (kk=k+1; gmmxring.in_rings[i][kk] != -1; kk++)
140 {
141 ib = gmmxring.in_rings[i][kk];
142 for (j=i+1; j < gmmxring.nbranch; j++)
143 {
144 for (l=0; gmmxring.in_rings[j][l] != -1; l++)
145 {
146 if (ia == gmmxring.in_rings[j][l] && ib == gmmxring.in_rings[j][l+1])
147 add_pair(i,j,ia,ib);
148 }
149 }
150 }
151 }
152 }
153 // now group fused systems
154 gmmxring.nfused = 0;
155 for (i=0; i < gmmxring.npairs; i++)
156 {
157 found = FALSE;
158 ia = gmmxring.pair[i][2];
159 ib = gmmxring.pair[i][3];
160 for (j=i+1; j < gmmxring.npairs; j++)
161 {
162 if (ia == gmmxring.pair[j][2] || ia == gmmxring.pair[j][3] || ib == gmmxring.pair[j][2] || ib == gmmxring.pair[j][3])
163 {
164 if (gmmxring.fused[ia] == -1)
165 {
166 found = TRUE;
167 gmmxring.fused[ia] = gmmxring.nfused;
168 gmmxring.fused[ib] = gmmxring.nfused;
169 gmmxring.fused[gmmxring.pair[j][2]] = gmmxring.nfused;
170 gmmxring.fused[gmmxring.pair[j][3]] = gmmxring.nfused;
171 gmmxring.nfused++;
172 } else
173 {
174 found = TRUE;
175 gmmxring.fused[ib] = gmmxring.fused[ia];
176 gmmxring.fused[gmmxring.pair[j][2]] = gmmxring.fused[ia];
177 gmmxring.fused[gmmxring.pair[j][3]] = gmmxring.fused[ia];
178 }
179 }
180 }
181 if (found == FALSE && gmmxring.fused[ia] == -1)
182 {
183 gmmxring.fused[ia] = gmmxring.nfused;
184 gmmxring.fused[ib] = gmmxring.nfused;
185 gmmxring.nfused++;
186 }
187 }
188 }
189 /* ==================================== */
190 void del_atom(int i, int **itemiat)
191 {
192 int j, itiat, k;
193
194 for (j=0; j < MAXIAT; j++)
195 {
196 if (itemiat[i][j] != 0)
197 {
198 itiat = itemiat[i][j];
199 for (k=0; k < MAXIAT; k++)
200 {
201 if (itemiat[itiat][k] == i)
202 itemiat[itiat][k] = 0;
203 }
204 }
205 }
206 for (j=0; j < MAXIAT; j++)
207 itemiat[i][j] = 0;
208 }
209 /* ==================================== */
210 void add_all_rings(int j,int *tmp)
211 {
212 int i,ia,k;
213 for (i=0; i < 10; i++)
214 {
215 if (gmmxring.in_rings[j][i] != -1)
216 {
217 ia = gmmxring.in_rings[j][i];
218 for (k=0; k < 10; k++)
219 {
220 if (tmp[k] == ia)
221 break;
222 else if (tmp[k] == -1)
223 {
224 tmp[k] = ia;
225 break;
226 }
227 }
228 }
229 }
230 }
231 /* ======================================= */
232 int icompare(int iatom, int *iarray1, int *iarray2)
233 {
234 int i;
235
236 gmmxsort(iatom, iarray1);
237 gmmxsort(iatom, iarray2);
238 for (i=0; i < iatom; i++)
239 {
240 if (iarray1[i] != iarray2[i])
241 return FALSE;
242 }
243 return TRUE;
244 }
245 /* ======================================= */
246 void gmmxsort(int iatom, int *iarray)
247 {
248 int i, j, itemp;
249
250 for (i=0; i < iatom-1; i++)
251 {
252 for(j=i+1; j < iatom; j++)
253 {
254 if (iarray[i] > iarray[j])
255 {
256 itemp = iarray[i];
257 iarray[i] = iarray[j];
258 iarray[j] = itemp;
259 }
260 }
261 }
262 }
263 /* =================================================== */
264 int in_ring(int ia,int k)
265 {
266 int i;
267 for (i=0; i < gmmxring.nring_atoms[k]; i++)
268 {
269 if (ia == gmmxring.ring_atoms[k][i])
270 return TRUE;
271 }
272 return FALSE;
273 }
274 /* =================================================== */
275 void add_pair(int ia,int ib,int k,int l)
276 {
277 int tmpa, tmpb;
278
279 if (ia < ib)
280 {
281 tmpa = ia;
282 tmpb = ib;
283 } else
284 {
285 tmpa = ib;
286 tmpb = ia;
287 }
288
289 // get here pair is not found create new pair
290 gmmxring.pair[gmmxring.npairs][0] = tmpa;
291 gmmxring.pair[gmmxring.npairs][1] = tmpb;
292 gmmxring.pair[gmmxring.npairs][2] = k;
293 gmmxring.pair[gmmxring.npairs][3] = l;
294 gmmxring.npairs++;
295 return;
296 }
297 // ============================================
298 void find_rings1(int *iexist,int **itemiat)
299 {
300 int i,j,k,ia,ib;
301 int rsize,ncount,numbond,nring,ringsize;
302 int *newatom,*atnum,*jj;
303 // int atbond[50][6];
304 int **atbond;
305 int *ring;
306 int n2count,ntrimset;
307 int *nodeN2,*trimSet,*nodesDone;
308 int smallestdegree,smallest;
309 int nodesToBreak;
310 int newring;
311
312 gmmxring.nrings = 0;
313 gmmxring.nbranch = 0;
314 // count size of ring system
315 ncount = 0;
316 for (i=1; i <= natom; i++)
317 if (iexist[i] != 0)
318 ncount++;
319
320 if (ncount < 2) return;
321 // allocate space
322 rsize = ncount;
323 newatom = ivector(0,natom+1);
324 atnum = ivector(0,rsize);
325 jj = ivector(0,rsize);
326 ring = ivector(0,rsize);
327 nodeN2 = ivector(0,rsize);
328 trimSet = ivector(0,rsize);
329 nodesDone = ivector(0,rsize);
330 atbond = imatrix(0,rsize+2,0,6);
331
332 ncount = 0;
333 for (i=1; i <= natom; i++)
334 {
335 if (iexist[i] != 0)
336 {
337 atnum[ncount] = i;
338 newatom[i] = ncount;
339 ncount++;
340 }
341 }
342 // count attachments and get bridgehead atoms
343 gmmxring.nbranch = 0;
344 numbond = 0;
345 for (i=0; i < rsize; i++)
346 {
347 nodeN2[i] = -1;
348 for (j=0; j < 6;j++)
349 atbond[i][j] = -1;
350 }
351 for (i=0; i < rsize; i++)
352 {
353 ia = atnum[i];
354 ncount = 0;
355 n2count = 0;
356 for (j=0; j < MAXIAT; j++)
357 {
358 if (itemiat[ia][j] != 0)
359 {
360 atbond[i][ncount] = newatom[itemiat[ia][j]];
361 ncount++;
362 if (ia < itemiat[ia][j])
363 numbond++;
364 }
365 }
366 jj[i] = ncount;
367 if (ncount >= 3)
368 {
369 gmmxring.branch_atoms[gmmxring.nbranch] = ia;
370 gmmxring.nbranch++;
371 } else if (ncount == 2)
372 {
373 nodeN2[n2count] = i;
374 n2count++;
375 }
376
377 }
378 // number of rings
379 nring = numbond - rsize + 1;
380 // loop over N2
381 ntrimset = 0;
382 for (i=0; i < rsize; i++)
383 trimSet[i] = -1;
384 do {
385
386 // check current connectivities
387 smallestdegree = 7;
388 smallest = -1;
389 for (i=0; i < rsize; i++)
390 {
391 ncount = 0;
392 for (j=0; j < 6; j++)
393 {
394 if (atbond[i][j] != -1)
395 ncount++;
396 }
397 jj[i] = ncount;
398
399 if (jj[i] == 1)
400 {
401 trim(i,&ntrimset,trimSet);
402 if (ntrimset > rsize)
403 message_alert("Error in TrimSet","Error");
404 for (j=0; j < 6; j++)
405 {
406 if (atbond[i][j] != -1)
407 {
408 ia = atbond[i][j];
409 for (k=0; j < 6; k++)
410 {
411 if (atbond[ia][k] == i)
412 {
413 atbond[i][j] = -1;
414 atbond[ia][k] = -1;
415 break;
416 }
417 }
418 }
419 }
420 } else if (jj[i] == 2)
421 {
422 check_N2(i,&n2count,nodeN2);
423 if (n2count > rsize)
424 message_alert("Error in check_N2","Error");
425 }
426
427 if (jj[i] < smallestdegree && jj[i] > 0)
428 {
429 smallest = i;
430 smallestdegree = jj[i];
431 }
432 }
433 if (smallest == -1)
434 break;
435 for(i=0; i < rsize; i++)
436 ring[i] = -1;
437 if (smallestdegree == 2)
438 {
439 nodesToBreak = 0;
440 for (i=0; i < n2count; i++)
441 {
442 newring = getRing(nodeN2[i],rsize,jj,atbond,ring);
443 if (newring)
444 {
445 for (j=0; ring[j] != -1 && j < rsize; j++)
446 ring[j] = atnum[ring[j]];
447 ringsize = j;
448 if (ringsize > rsize)
449 message_alert("Error in getRing","Error");
450 log_ring1(ringsize,ring);
451 }
452 nodesDone[nodesToBreak] = nodeN2[i];
453 nodesToBreak++;
454 }
455 for (i=0; i < nodesToBreak; i++)
456 {
457 breakBond(nodesDone[i],atbond);
458 nodeN2[i] = -1;
459 }
460 n2count = 0;
461 } else if (smallestdegree == 3)
462 {
463 newring = getRing(smallest,rsize,jj,atbond,ring);
464 if (newring)
465 {
466 ia = ring[0];
467 ib = ring[1];
468
469 for (j=0; ring[j] != -1 && j < rsize; j++)
470 ring[j] = atnum[ring[j]];
471 ringsize = j;
472 log_ring1(ringsize,ring);
473 delete_bond(ia,ib,atbond);
474 }
475 }
476 repack(rsize,atbond);
477 } while (ntrimset < rsize);
478 // free space
479 free_ivector(newatom ,0,natom+1);
480 free_ivector(atnum ,0,rsize);
481 free_ivector(jj ,0,rsize);
482 free_ivector(ring ,0,rsize);
483 free_ivector(nodeN2 ,0,rsize);
484 free_ivector(trimSet ,0,rsize);
485 free_ivector(nodesDone ,0,rsize);
486 free_imatrix(atbond ,0,rsize+2,0,6);
487 }
488 // ============= queque stuff =======================
489 #define MAX_QUE_SIZE 100
490 int queue[MAX_QUE_SIZE];
491
492 void addq(int front,int *rear,int item)
493 {
494 *rear = (*rear+1)%MAX_QUE_SIZE;
495 if (front == *rear)
496 {
497 message_alert("Error in Queue","Error");
498 return;
499 }
500 queue[*rear] = item;
501 }
502 int deleteq(int *front,int rear)
503 {
504 if (*front == rear) return(-1);
505 *front = (*front+1)%MAX_QUE_SIZE;
506 return queue[*front];
507 }
508 // ===================================================
509 int getRing(int ia,int rsize,int *jj, int **atbond,int *ring)
510 {
511 int i,j,ib,ic,intersect,itmp;
512 int front,rear;
513 int **path;
514
515 path = imatrix(0,rsize,0,rsize);
516 for (i=0; i < rsize; i++)
517 {
518 for (j=0; j < rsize; j++)
519 path[i][j] = -1;
520 }
521 front = rear = -1;
522 // get all neighbors of root node, put them on queque
523 // and update path
524 for (i=0; i < jj[ia]; i++)
525 {
526 itmp = atbond[ia][i];
527 if (jj[atbond[ia][i]] > 0)
528 {
529 addq(front,&rear,atbond[ia][i]);
530 path[ia][0] = ia;
531 path[atbond[ia][i]][0] = ia;
532 path[atbond[ia][i]][1] = atbond[ia][i];
533 }
534 }
535 while (front >= -1)
536 {
537 ib = deleteq(&front,rear); // get new node from queue
538 if (ib == -1) // empty queue
539 {
540 free_imatrix(path,0,rsize,0,rsize);
541 return FALSE;
542 }
543 for (i=0; i < jj[ib]; i++)
544 {
545 ic = atbond[ib][i];
546 if (path[ic][0] != -1) // already visited
547 {
548 for (j=0; path[ib][j] != -1; j++)
549 ;
550 if (path[ib][j-2] != ic) // found collision
551 {
552 intersect = get_intersection(path[ib],path[ic]);
553 if (intersect == 1)
554 {
555 for(j=0; j < rsize; j++)
556 ring[j] = -1;
557 getUnion(path[ib],path[ic],ring);
558 free_imatrix(path,0,rsize,0,rsize);
559 return TRUE;
560 }
561 }
562 } else // put on queue
563 {
564 addq(front,&rear,ic);
565 for (j=0;path[ib][j] != -1; j++)
566 path[ic][j] = path[ib][j];
567 path[ic][j] = ic;
568 }
569 }
570
571 }
572 free_imatrix(path,0,rsize,0,rsize);
573 return FALSE;
574 }
575 // ======================================
576 int get_intersection(int *path1, int *path2)
577 {
578 int i,j,ia,icount;
579 icount = 0;
580 for (i=0; path1[i] != -1; i++)
581 {
582 ia = path1[i];
583 for (j=0; path2[j] != -1; j++)
584 {
585 if (path2[j] == ia)
586 icount++;
587 }
588 }
589 return icount;
590 }
591 // ======================================
592 void getUnion(int *path1, int *path2,int *ring)
593 {
594 int i,j,k,icount;
595 icount = 0;
596
597 for (i=0; path1[i] != -1; i++)
598 {
599 ring[icount] = path1[i];
600 icount++;
601 }
602 for (j=0; path2[j] != -1; j++)
603 ;
604 for (k=j-1; k > 0; k--)
605 {
606 ring[icount] = path2[k];
607 icount++;
608 }
609 }
610 // ======================================
611 void log_ring1(int iatom, int *path)
612 {
613 int i, isame, j;
614 int iarray1[50], iarray2[50];
615
616 if (iatom > MAX_RINGSIZE)
617 {
618 message_alert("Ring too large in find ring","Error");
619 return;
620 }
621 for (i=0; i < gmmxring.nrings; i++)
622 {
623 if (gmmxring.nring_atoms[i] == iatom)
624 {
625 isame = FALSE;
626 for (j=0; j < 50; j++)
627 {
628 iarray1[j] = 0;
629 iarray2[j] = 0;
630 }
631 for(j=0; j < iatom; j++)
632 iarray1[j] = path[j];
633 for(j=0; j < gmmxring.nring_atoms[i]; j++)
634 iarray2[j] = gmmxring.ring_atoms[i][j];
635
636 isame = icompare(iatom,iarray1, iarray2);
637 if (isame == TRUE)
638 return;
639 }
640 }
641 gmmxring.nring_atoms[gmmxring.nrings] = iatom;
642 for (i=0; i < iatom; i++)
643 gmmxring.ring_atoms[gmmxring.nrings][i] = path[i];
644 gmmxring.nrings++;
645 }
646 // =================================
647 void delete_bond(int ia,int ib, int **atbond)
648 {
649 int i;
650
651 for (i=0; i < 6; i++)
652 {
653 if (atbond[ia][i] == ib)
654 atbond[ia][i] = -1;
655 if (atbond[ib][i] == ia)
656 atbond[ib][i] = -1;
657 }
658 }
659 // =================================
660 void check_N2(int ia,int *n2count,int *nodes)
661 {
662 int i;
663
664 for (i=0; i < *n2count; i++)
665 {
666 if (nodes[i] == ia) // already there
667 return;
668 }
669 nodes[*n2count] = ia;
670 *n2count += 1;
671 }
672 // =================================
673 void trim(int ia,int *j,int *trimset)
674 {
675 int i;
676
677 for (i=0; i < *j; i++)
678 {
679 if (trimset[i] == ia)
680 return;
681 }
682 trimset[*j] = ia;
683 *j += 1;
684 }
685 // =================================
686 void breakBond(int ia,int **atbond)
687 {
688 int i,j,ib;
689
690 for (i=0; i < 6; i++)
691 {
692 if (atbond[ia][i] != -1)
693 {
694 ib = atbond[ia][i];
695 for (j=0; j < 6; j++)
696 {
697 if (atbond[ib][j] == ia)
698 {
699 atbond[ia][i] = -1;
700 atbond[ib][j] = -1;
701 return;
702 }
703 }
704 }
705 }
706 }
707 // =====================================
708 void repack(int rsize,int **atbond)
709 {
710 int i,j,k;
711
712 for (i=0; i < rsize; i++)
713 {
714 for (j=0; j < 5; j++)
715 {
716 if (atbond[i][j] == -1)
717 {
718 for (k=j+1; k < 6; k++)
719 atbond[i][k-1] = atbond[i][k];
720 atbond[i][5] = -1;
721 }
722 }
723 }
724 }
725