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, 11 months ago) by tjod
File size: 19129 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #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