ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/asnsym.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 11 months ago) by tjod
Original Path: trunk/smi23d/src/mengine/asnsym.c
File size: 53717 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    
7     #define LINEAR 1
8     #define ASYMMETRIC 2
9     #define PROLATE 3
10     #define OBLATE 4
11     #define SPHERICAL 5
12     #define PLANAR 6
13    
14     struct {
15     int ishape,numsym,chiral;
16     double xi,yi,zi;
17     } symmetry;
18    
19     double TOLER = 0.002;
20     double TOL2 = 0.00002;
21    
22     double sign(double,double);
23     void message_alert(char *,char *);
24     int tstc3(double **,double **,int,int,int,double *);
25     int tstc4(double **,double **,int,int,int,double *);
26     int tstc5(double **,double **,int,int,int,double *);
27     void secmom(double **, double *, double *, double *, int *, double (*)[3]);
28     void orients(double **, double **,double *, double *,int,int);
29     void tcopy(double **, double **,int);
30     int findh(double **, double **, int);
31     int findcn(double **,double **,double *, double *, int *);
32     int findc2(double **,double **, double **, int *);
33     int findv(double **,double **, double **, int *);
34     int findsn(double **,double **, double **, int *);
35     int findi(double **,double **);
36     void tform(double (*)[3], double **, double **,int);
37     void reflect(double **, double **,int);
38     int tcomp(double **,double **);
39     void circset(double **, double **,int,int *,int *);
40     void rotates(double **a,double **b,int,double);
41     void hdel(int);
42     void sphere(double **,double **, double **,int *, int *);
43     void sym_invert(double **, double **);
44     void triang(double **,int,int,int,double *, double *,double *,double *,double *,double *);
45     void zero_array(double **,int);
46     void hqrii(double *,int,int,double *,double (*)[3]);
47     void vibsym(int);
48     void ptgrp(int,char *,double **);
49     void sphset(double **, double **,int *, int *,int *);
50    
51     void zero_array(double **a,int num)
52     {
53     int i;
54     for (i=1; i <= num; i++)
55     {
56     a[i][0] = 0.0;
57     a[i][1] = 0.0;
58     a[i][2] = 0.0;
59     }
60     }
61     void triang(double **a,int iat,int jat,int kat,double *alpha, double *beta,
62     double *gamma,double *dij, double *dik, double *djk)
63     {
64     double xi,yi,zi,xj,yj,zj,xk,yk,zk,dotjk,dotik,dotij;
65    
66     xi = a[iat][0];
67     yi = a[iat][1];
68     zi = a[iat][2];
69     xj = a[jat][0];
70     yj = a[jat][1];
71     zj = a[jat][2];
72     xk = a[kat][0];
73     yk = a[kat][1];
74     zk = a[kat][2];
75     *dij = sqrt( (xi-xj)*(xi-xj) + (yi-yj)*(yi-yj) + (zi-zj)*(zi-zj) );
76     *dik = sqrt( (xi-xk)*(xi-xk) + (yi-yk)*(yi-yk) + (zi-zk)*(zi-zk) );
77     *djk = sqrt( (xj-xk)*(xj-xk) + (yj-yk)*(yj-yk) + (zj-zk)*(zj-zk) );
78     dotjk = (xj-xi)*(xk-xi) + (yj-yi)*(yk-yi) + (zj-zi)*(zk-zi);
79     dotik = (xi-xj)*(xk-xj) + (yi-yj)*(yk-yj) + (zi-zj)*(zk-zj);
80     dotij = (xi-xk)*(xj-xk) + (yi-yk)*(yj-yk) + (zi-zk)*(zj-zk);
81     *alpha = acos(dotjk/(*dij**dik));
82     *beta = acos(dotik/(*dij**djk));
83     *gamma = acos(dotij/(*dik**djk));
84    
85     }
86     /* ================================================= */
87     int tstc3(double **a,double **b,int iat,int jat,int kat,double *center)
88     {
89     int itst, iz;
90     double p1[3];
91     double phi3,theta3,fact3;
92     double alpha,beta,gamma,dij,dik,djk;
93     double px,py,pz;
94    
95     p1[0] = p1[1] = p1[2] = 0.0;
96     itst = 0;
97     phi3 = (8.0/3.0)*atan(1.0);
98     theta3 = 0.5*phi3;
99     fact3 = 2.0/3.0;
100    
101     triang(a,iat,jat,kat,&alpha,&beta,&gamma,&dij,&dik,&djk);
102    
103     if (fabs(alpha-theta3) < TOLER && fabs(dij-dik) < TOLER)
104     {
105     px = 0.5*(a[jat][0] + a[kat][0]);
106     py = 0.5*(a[jat][1] + a[kat][1]);
107     pz = 0.5*(a[jat][2] + a[kat][2]);
108     center[0] = a[iat][0] + fact3*(px - a[iat][0]);
109     center[1] = a[iat][1] + fact3*(py - a[iat][1]);
110     center[2] = a[iat][2] + fact3*(pz - a[iat][2]);
111    
112     orients(a,b,p1,center,3,natom);
113     tcopy(b,a,natom);
114     rotates(a,b,3,phi3);
115     iz = tcomp(a,b);
116     return (iz);
117     }
118     return FALSE;
119     }
120     int tstc4(double **a,double **b,int iat,int jat,int kat,double *center)
121     {
122     int itst, iz;
123     double p1[3];
124     double alpha,beta,gamma,dij,dik,djk;
125     double halfpi;
126    
127     p1[0] = p1[1] = p1[2] = 0.0;
128     itst = 0;
129     halfpi = 2.0*atan(1.0);
130    
131     triang(a,iat,jat,kat,&alpha,&beta,&gamma,&dij,&dik,&djk);
132    
133     if (fabs(alpha-halfpi) < TOLER && fabs(dij-dik) < TOLER)
134     {
135     center[0] = 0.5*(a[jat][0] + a[kat][0]);
136     center[1] = 0.5*(a[jat][1] + a[kat][1]);
137     center[2] = 0.5*(a[jat][2] + a[kat][2]);
138     } else if (fabs(beta-halfpi) < TOLER && fabs(dij-djk) < TOLER)
139     {
140     center[0] = 0.5*(a[iat][0] + a[kat][0]);
141     center[1] = 0.5*(a[iat][1] + a[kat][1]);
142     center[2] = 0.5*(a[iat][2] + a[kat][2]);
143     } else if (fabs(gamma-halfpi) < TOLER && fabs(dik-djk) < TOLER)
144     {
145     center[0] = 0.5*(a[jat][0] + a[iat][0]);
146     center[1] = 0.5*(a[jat][1] + a[iat][1]);
147     center[2] = 0.5*(a[jat][2] + a[iat][2]);
148     } else if (fabs(dij-dik) < TOLER && fabs(dij-djk) < TOLER && TOLER && fabs(dik-djk))
149     {
150     center[0] = a[iat][0];
151     center[1] = a[iat][1];
152     center[2] = a[iat][2];
153     } else
154     return FALSE;
155    
156     orients(a,b,p1,center,3,natom);
157     tcopy(b,a,natom);
158     rotates(a,b,3,halfpi);
159     iz = tcomp(a,b);
160     return (iz);
161    
162     }
163     int tstc5(double **a,double **b,int iat,int jat,int kat,double *center)
164     {
165     int itst, iz;
166     double p1[3];
167     double piovr4, phi5,theta5,fact5;
168     double alpha,beta,gamma,dij,dik,djk;
169     double px,py,pz;
170    
171     p1[0] = p1[1] = p1[2] = 0.0;
172     itst = 0;
173     piovr4 = atan(1.0);
174     phi5 = 1.6*piovr4;
175     theta5 = 2.4*piovr4;
176     fact5 = 1.0/(2.0*sin(0.8*piovr4)*sin(0.8*piovr4));
177    
178     triang(a,iat,jat,kat,&alpha,&beta,&gamma,&dij,&dik,&djk);
179    
180     if (fabs(alpha-theta5) < TOLER && fabs(dij-dik) < TOLER)
181     {
182     px = 0.5*(a[jat][0] + a[kat][0]);
183     py = 0.5*(a[jat][1] + a[kat][1]);
184     pz = 0.5*(a[jat][2] + a[kat][2]);
185     center[0] = a[iat][0] + fact5*(px - a[iat][0]);
186     center[1] = a[iat][1] + fact5*(py - a[iat][1]);
187     center[2] = a[iat][2] + fact5*(pz - a[iat][2]);
188     } else if (fabs(beta-theta5) < TOLER && fabs(dij-djk) < TOLER)
189     {
190     px = 0.5*(a[iat][0] + a[kat][0]);
191     py = 0.5*(a[iat][1] + a[kat][1]);
192     pz = 0.5*(a[iat][2] + a[kat][2]);
193     center[0] = a[jat][0] + fact5*(px - a[jat][0]);
194     center[1] = a[jat][1] + fact5*(py - a[jat][1]);
195     center[2] = a[jat][2] + fact5*(pz - a[jat][2]);
196     } else if (fabs(gamma-theta5) < TOLER && fabs(dik-djk) < TOLER)
197     {
198     px = 0.5*(a[jat][0] + a[iat][0]);
199     py = 0.5*(a[jat][1] + a[iat][1]);
200     pz = 0.5*(a[jat][2] + a[iat][2]);
201     center[0] = a[kat][0] + fact5*(px - a[kat][0]);
202     center[1] = a[kat][1] + fact5*(py - a[kat][1]);
203     center[2] = a[kat][2] + fact5*(pz - a[kat][2]);
204     } else
205     return FALSE;
206     orients(a,b,p1,center,3,natom);
207     tcopy(b,a,natom);
208     rotates(a,b,3,phi5);
209     iz = tcomp(a,b);
210     return (iz);
211     }
212     /* ====================================== */
213     void sphset(double **a, double **aset,int *npop, int *nset,int *numset)
214     {
215     int i,j1,j,iset, num;
216     double xtmp,ytmp,ztmp;
217    
218     for (i=1; i <= natom; i++)
219     {
220     aset[i][0] = i;
221     aset[i][1] = atom[i].atomwt;
222     xtmp = a[i][0];
223     ytmp = a[i][1];
224     ztmp = a[i][2];
225     aset[i][2] = sqrt(xtmp*xtmp+ytmp*ytmp+ztmp*ztmp);
226     }
227     iset = 0;
228     for (i=1; i < natom; i++)
229     {
230     num = 0;
231     if (aset[i][0] != 0.0)
232     {
233     iset++;
234     num++;
235     nset[i] = iset;
236     aset[i][0] = 0.0;
237     j1 = i + 1;
238     for (j=j1; j <= natom; j++)
239     {
240     if (aset[j][0] != 0.0)
241     {
242     if (fabs(aset[j][1]-aset[i][1]) <= TOLER &&
243     fabs(aset[j][2]-aset[i][2]) <= TOLER)
244     {
245     num++;
246     nset[j] = iset;
247     aset[j][0] = 0.0;
248     }
249     }
250     }
251     npop[iset] = num;
252     }
253     }
254     *numset = iset;
255     }
256     void sphere(double **a, double **b, double **c, int *nset, int *norder)
257     {
258     int i,j,k,keyset, numset, minnum;
259     int mpop, num3, i2, j1, j2, k1, iat,jat,kat;
260     int iz,iz1,iz2,key,ihop,ixyz;
261     double savez,curz,x,y,z,cmax,theta,tst;
262     int *npop;
263     double halfpi,pi;
264     double p1[3],center[3],save1[3],save2[3];
265    
266     halfpi = 2.0*atan(1.0);
267     pi = 2.0*halfpi;
268     p1[0] = p1[1] = p1[2] = 0.0;
269    
270     keyset = 0;
271     npop = ivector(0,natom+1);
272     sphset(a,c,npop,nset,&numset);
273    
274     minnum = natom;
275     *norder = 0;
276     for (i=1; i <= numset; i++)
277     {
278     j = npop[i];
279     if (npop[i] <= minnum && npop[i] > 3)
280     {
281     minnum = npop[i];
282     keyset = i;
283     }
284     }
285    
286     key = natom;
287     mpop = 0;
288     for (i=1; i <= natom; i++)
289     {
290     if (nset[i] == keyset)
291     {
292     mpop++;
293     if (key < i)
294     key = i;
295     nset[mpop] = i;
296     }
297     }
298    
299     if (mpop < 3)
300     {
301     *norder = 0;
302     free_ivector(npop,0,natom+1);
303     return;
304     }
305     num3 = 0;
306     i2 = mpop-2;
307     j2 = mpop-1;
308     for (i=1; i <= i2; i++)
309     {
310     iat = nset[i];
311     j1 = i + 1;
312     for (j=j1; j <= j2; j++)
313     {
314     jat = nset[j];
315     k1 = j+1;
316     for (k=k1; k <= mpop; k++)
317     {
318     kat = nset[k];
319     if (*norder <= 3)
320     {
321     tcopy(a,c,natom);
322     iz = tstc5(c,b,iat,jat,kat,center);
323     if (iz == TRUE)
324     {
325     *norder = 5;
326     savez = fabs(c[key][2]);
327     save1[0] = center[0];
328     save1[1] = center[1];
329     save1[2] = center[2];
330     }else
331     {
332     tcopy(a,c,natom);
333     zero_array(b,natom);
334     iz1 = tstc4(c,b,iat,jat,kat,center);
335     if (iz1 == TRUE)
336     {
337     *norder = 4;
338     tcopy(c,a,natom);
339     } else
340     {
341     if (num3 != 2)
342     {
343     tcopy(a,c,natom);
344     iz2 = tstc3(c,b,iat,jat,kat,center);
345     if (iz2 == TRUE)
346     {
347     ihop = num3+1;
348     if (ihop == 1)
349     {
350     *norder = 3;
351     num3 = 1;
352     save1[0] = center[0];
353     save1[1] = center[1];
354     save1[2] = center[2];
355     } else if (ihop == 2)
356     {
357     num3 = 2;
358     save2[0] = center[0];
359     save2[1] = center[1];
360     save2[2] = center[2];
361     }
362     }
363     }
364     }
365     }
366     } else if (*norder == 5)
367     {
368     tcopy(a,c,natom);
369     iz = tstc5(c,b,iat,jat,kat,center);
370     if (iz == TRUE)
371     {
372     curz = fabs(c[key][2]);
373     if (fabs(curz-savez) > TOLER && (savez < curz ))
374     {
375     savez = curz;
376     save1[0] = center[0];
377     save1[1] = center[1];
378     save1[2] = center[2];
379     }
380     }
381     } else if (*norder == 4)
382     {
383     tcopy(a,c,natom);
384     iz = tstc4(c,b,iat,jat,kat,center);
385     if (iz == TRUE)
386     {
387     if (fabs(center[2]) < TOLER)
388     {
389     x = center[0];
390     y = center[1];
391     theta = halfpi;
392     if (fabs(y) > TOLER) theta = atan(x/y);
393     rotates(c,a,3,theta);
394     goto L_5;
395     }
396     }
397     }
398     }
399     }
400     }
401     L_5:
402     if (*norder == 0)
403     {
404     center[0] = 0.5*(a[iat][0] + a[jat][0]);
405     center[1] = 0.5*(a[iat][1] + a[jat][1]);
406     center[2] = 0.5*(a[iat][2] + a[jat][2]);
407     orients(a,b,p1,center,3,natom);
408     x = b[iat][0];
409     y = b[iat][1];
410     theta = halfpi;
411     if (fabs(y) > TOLER) theta = atan(x/y);
412     rotates(b,a,3,theta);
413     } else if (*norder == 3)
414     {
415     center[0] = 0.5*(save1[0] + save2[0]);
416     center[1] = 0.5*(save1[1] + save2[1]);
417     center[2] = 0.5*(save1[2] + save2[2]);
418     orients(a,b,p1,center,3,natom);
419     tcopy(b,a,natom);
420     b[0][0] = save1[0];
421     b[0][1] = save1[1];
422     b[0][2] = save1[2];
423     b[1][0] = save2[0];
424     b[1][1] = save2[1];
425     b[1][2] = save2[2];
426     orients(b,c,p1,center,3,2);
427     tcopy(c,b,2);
428     save1[0] = b[0][0];
429     save1[1] = b[0][1];
430     save1[2] = b[0][2];
431     save2[0] = b[1][0];
432     save2[1] = b[1][1];
433     save2[2] = b[1][2];
434     x = 0.5*(save1[0]-save2[1]);
435     y = 0.5*(save1[1]+save2[0]);
436     theta = halfpi;
437     if (fabs(y) > TOLER) theta = atan(x/y);
438     rotates(a,b,3,theta);
439     tcopy(b,a,natom);
440    
441     } else if (*norder == 4)
442     {
443     x = a[key][0];
444     y = a[key][1];
445     z = a[key][2];
446     cmax = fabs(z);
447     ixyz = 3;
448     for (i=1; i <= 2; i++)
449     {
450     tst = fabs(a[key][i]);
451     if (fabs(tst-cmax) > TOLER && (cmax <= tst))
452     {
453     ixyz = i;
454     cmax = tst;
455     }
456     }
457     if (ixyz < 3)
458     {
459     if (ixyz == 2)
460     ixyz = 1;
461     else
462     ixyz = 2;
463     rotates(a,b,ixyz,halfpi);
464     tcopy(b,a,natom);
465     }
466     if (a[key][3] < 0.0)
467     {
468     rotates(a,b,1,pi);
469     tcopy(b,a,natom);
470     }
471     } else if (*norder == 5)
472     {
473     orients(a,b,p1,save1,3,natom);
474     tcopy(b,a,natom);
475     if (a[key][3] < 0.0)
476     {
477     rotates(a,b,1,pi);
478     tcopy(b,a,natom);
479     }
480     }
481    
482    
483     free_ivector(npop,0,natom+1);
484     }
485     void sym_invert(double **a, double **b)
486     {
487     int i,j;
488     double t[3][3];
489     for (i=0; i < 3; i++)
490     {
491     for (j=0; j <3; j++)
492     t[i][j] = 0.0;
493     }
494     t[0][0] = -1.0;
495     t[1][1] = -1.0;
496     t[2][2] = -1.0;
497     tform(t,a,b,natom);
498     }
499     void reflect(double **a, double **b, int iaxis)
500     {
501     int i,j;
502     double t[3][3];
503     for (i=0; i < 3; i++)
504     {
505     for (j=0; j <3; j++)
506     t[i][j] = 0.0;
507     t[i][i] = 1.0;
508     }
509     t[iaxis-1][iaxis-1] = -1.0;
510     tform(t,a,b,natom);
511     }
512     /* ==================================== */
513     int tcomp(double **a, double **b)
514     {
515     int i,j, k;
516     int *nrs;
517     double dsum, x1, x2;
518    
519     nrs = ivector(0,natom+1);
520     for (i=0; i <= natom; i++)
521     nrs[i] = 0;
522    
523     for (i=1; i <= natom; i++)
524     {
525     for (j=1; j <= natom; j++)
526     {
527     if (fabs(atom[i].atomwt-atom[j].atomwt) <= TOLER)
528     {
529     if (nrs[i] == 0)
530     {
531     dsum = 0;
532     for (k=0; k < 3; k++)
533     {
534     x1 = a[i][k];
535     x2 = b[j][k];
536     dsum += (a[i][k]-b[j][k])*(a[i][k]-b[j][k]);
537     }
538     if (dsum <= TOLER)
539     {
540     nrs[i] = j;
541     goto L_5;
542     }
543     }
544     }
545     }
546     free_ivector(nrs,0,natom+1);
547     return FALSE;
548     L_5:
549     continue;
550     }
551     free_ivector(nrs,0,natom+1);
552     return TRUE;
553     }
554     void tcopy(double **a, double **b, int num)
555     {
556     int i,j;
557     for (i=1; i <= num; i++)
558     {
559     for (j=0; j < 3; j++)
560     b[i][j] = a[i][j];
561     }
562     }
563     void tform(double t[][3], double **a, double **b, int num)
564     {
565     int i;
566     for (i=1; i <= num; i++)
567     {
568     b[i][0] = t[0][0]*a[i][0] + t[0][1]*a[i][1] + t[0][2]*a[i][2];
569     b[i][1] = t[1][0]*a[i][0] + t[1][1]*a[i][1] + t[1][2]*a[i][2];
570     b[i][2] = t[2][0]*a[i][0] + t[2][1]*a[i][1] + t[2][2]*a[i][2];
571     }
572     }
573     /* ==================================== */
574     void circset(double **a, double **aset,int ixyz,int *nset,int *numset)
575     {
576     int i,j, i1, i2, i3, iset,j1;
577     double q2, q3, an,ap,ad;
578    
579     i1 = ixyz;
580     i2 = ixyz%3+1;
581     i3 = i2%3 +1;
582     i1--;
583     i2--;
584     i3--;
585     for (i=1; i <= natom; i++)
586     {
587     aset[i][0] = atom[i].atomwt;
588     aset[i][1] = a[i][i1];
589     q2 = a[i][i2];
590     q3 = a[i][i3];
591     aset[i][2] = sqrt(q2*q2+q3*q3);
592     }
593     // set 0 on axis atoms
594     for (i=1; i <= natom; i++)
595     {
596     if (fabs(aset[i][2]) <= TOLER)
597     {
598     nset[i] = 0;
599     aset[i][0] = 0.0;
600     }
601     }
602     // remaining sets
603     iset = 0;
604     for (i=1; i < natom; i++) // loop is natom-1
605     {
606     if (aset[i][0] != 0.0)
607     {
608     iset++;
609     nset[i] = iset;
610     an = aset[i][0];
611     ap = aset[i][1];
612     ad = aset[i][2];
613     aset[i][0] = 0.0;
614     j1 = i+1;
615     for (j=j1; j <= natom; j++)
616     {
617     if (fabs(aset[j][0]-an) <= TOL2 &&
618     fabs(aset[j][1]-ap) <= TOLER &&
619     fabs(aset[j][2]-ad) <= TOLER)
620     {
621     nset[j] = iset;
622     aset[j][0] = 0.0;
623     }
624     }
625     }
626     }
627     *numset = iset;
628     for (i=1; i <= natom; i++)
629     aset[i][0] = atom[i].atomwt;
630     }
631     /* =============================================== */
632     void rotates(double **a,double **b,int iaxis,double theta)
633     {
634     int i,j,i1,i2,i3;
635     double s,c;
636     double t[3][3];
637    
638     i1 = iaxis;
639     i2 = iaxis%3 +1;
640     i3 = i2%3 + 1;
641     i1--;
642     i2--;
643     i3--;
644     s = sin(theta);
645     c = cos(theta);
646     for (i=0; i < 3; i++)
647     {
648     for (j=0; j < 3; j++)
649     t[i][j] = 0.0;
650     }
651     t[i1][i1] = 1.0;
652     t[i2][i2] = c;
653     t[i2][i3] = s;
654     t[i3][i2] = -s;
655     t[i3][i3] = c;
656     tform(t,a,b,natom);
657     }
658     /* =============================================== */
659     void orients(double **a, double **b, double *p1, double *p2, int nset, int num)
660     {
661     int i,j, i1, i2, i3;
662     double trvec[3], t[3][3];
663     double v1, v2, v3, vnorm;
664     double alph, beta, gamm, v2v2, v3v3, v2233;
665     double tol = 1.0e-17;
666    
667     trvec[0] = p1[0];
668     trvec[1] = p1[1];
669     trvec[2] = p1[2];
670    
671     for (i=1; i <= num; i++)
672     {
673     for (j=0; j < 3; j++)
674     a[i][j] -= trvec[j];
675     }
676    
677     i1 = nset;
678     i2 = i1%3 + 1;
679     i3 = i2%3 + 1;
680     i1--;
681     i2--;
682     i3--;
683     v1 = p2[i1] - p1[i1];
684     v2 = p2[i2] - p1[i2];
685     v3 = p2[i3] - p1[i3];
686     vnorm = sqrt(v1*v1 + v2*v2 + v3*v3);
687     if (vnorm == 0.0)
688     return;
689     if ( fabs(v1) < TOLER && fabs(v2) < TOLER && fabs(v3) < TOLER)
690     return;
691     if (fabs(v2) < tol && fabs(v3) < tol)
692     {
693     if (v1 > 0.0)
694     {
695     t[i1][i1] = 1.0;
696     t[i2][i2] = 1.0;
697     t[i3][i3] = 1.0;
698     t[i1][i2] = t[i2][i1] = 0.0;
699     t[i1][i3] = t[i3][i1] = 0.0;
700     t[i3][i2] = t[i2][i3] = 0.0;
701     } else
702     {
703     t[i1][i1] = -1.0;
704     t[i2][i2] = 1.0;
705     t[i3][i3] = 1.0;
706     t[i1][i2] = t[i2][i1] = 0.0;
707     t[i1][i3] = t[i3][i1] = 0.0;
708     t[i3][i2] = t[i2][i3] = 0.0;
709     }
710     } else
711     {
712     alph = v1/vnorm;
713     beta = v2/vnorm;
714     gamm = v3/vnorm;
715     v2v2 = v2*v2;
716     v3v3 = v3*v3;
717     v2233 = 1.0/(v2v2+v3v3);
718     t[i1][i1] = alph;
719     t[i1][i2] = beta;
720     t[i1][i3] = gamm;
721     t[i2][i1] = -t[i1][i2];
722     t[i3][i1] = -t[i1][i3];
723     t[i2][i3] = v2*v3*(alph-1.0)*v2233;
724     t[i3][i2] = t[i2][i3];
725     t[i2][i2] = (v2v2*alph +v3v3)*v2233;
726     t[i3][i3] = (v3v3*alph +v2v2)*v2233;
727     }
728     tform(t,a,b,num);
729     for (i=1; i <= num; i++)
730     {
731     for (j=0; j < 3; j++)
732     b[i][j] += trvec[j];
733     }
734     }
735     /* ================================== */
736     int findi(double **a, double **b)
737     {
738     int iz;
739     sym_invert(a,b);
740     iz = tcomp(a,b);
741     return (iz);
742     }
743     int findh(double **a, double **b, int nset)
744     {
745     int iz;
746     reflect(a,b,nset);
747     iz = tcomp(a,b);
748     return (iz);
749     }
750     int findcn(double **a, double **b,double *p1, double *p2, int *maxcn)
751     {
752     int i,j,k,l, iz, iaxis;
753     int i1, i2;
754     double halfpi, pi,theta;
755     double t[3][3];
756    
757     halfpi = 2.0*atan(1.0);
758     pi = 2.0*halfpi;
759     iaxis = 0;
760     for (i=6; i >= 2; i--)
761     {
762     theta = 2.0*pi/(float)i;
763     for (j=0; j < 3; j++)
764     {
765     for (k=0; k < 3; k++)
766     {
767     for (l=0; l < 3; l++)
768     {
769     t[k][l] = 0.0;
770     }
771     }
772     i1 = j+1;
773     i2 = j+2;
774     if (i1 > 2) i1 -= 3;
775     if (i2 > 2) i2 -= 3;
776     t[i1][i1] = cos(theta);
777     t[i2][i2] = cos(theta);
778     t[i1][i2] = sin(theta);
779     t[i2][i1] = -sin(theta);
780     t[j][j] = 1.0;
781     tform(t,a,b,natom);
782     iz = tcomp(a,b);
783     if (iz == TRUE)
784     {
785     *maxcn = i;
786     for (k=0; k < 3; k++)
787     {
788     p1[k] = 0.0;
789     p2[k] = 0.0;
790     if (k == j) p2[k] = 1.0;
791     }
792     orients(b,a,p1,p2,3,natom);
793     return TRUE;
794     }
795     }
796     }
797     *maxcn = 1;
798     return FALSE;
799     }
800     int findsn(double **a, double **b, double **c, int *maxcn)
801     {
802     int iz;
803     double pi,theta;
804     pi = 4.0*atan(1.0);
805     theta = pi/(*maxcn);
806     rotates(a,b,3,theta);
807     reflect(b,c,3);
808     iz = tcomp(a,c);
809     return iz;
810     }
811     int findc2(double **a, double **b, double **aset, int *nset)
812     {
813     int i,j,k, j1, numset, iset, jset, iz;
814     double halfpi, pi, x,y,theta, xi,yi;
815     double proi, ani,disi;
816    
817     halfpi = 2.0*atan(1.0);
818     pi = 2.0*halfpi;
819    
820     circset(a,aset,3,nset,&numset);
821     for (i=1; i <= numset; i++)
822     {
823     for (j=1; j <= natom; j++)
824     {
825     if(nset[j] == i) // atom belongs to set
826     {
827     if ( fabs(aset[j][1]) < TOLER )
828     {
829     j1 = j+1;
830     for (k=j1; k <= natom; k++)
831     {
832     if (nset[k] == i)
833     {
834     x = (a[j][0]+a[k][0])*0.5;
835     y = (a[j][1]+a[k][1])*0.5;
836     theta = halfpi;
837     if (fabs(y) > TOLER) theta = -atan(x/y);
838     rotates(a,b,3,theta);
839     rotates(b,aset,2,pi);
840     iz = tcomp(b,aset);
841     if (iz == TRUE)
842     {
843     tcopy(b,a,natom);
844     return TRUE;
845     }
846     }
847     }
848     }
849     }
850     }
851     }
852     //
853     circset(a,aset,3,nset,&numset);
854     iset = 1;
855     for (i=1; i <= natom; i++)
856     if (nset[i] == iset) goto L_10;
857     return FALSE;
858     //
859     L_10:
860     proi = aset[i][1];
861     ani = aset[i][0];
862     disi = aset[i][2];
863     xi = a[i][0];
864     yi = a[i][1];
865     j1 = iset+1;
866     for (j=j1; j <= numset; j++)
867     {
868     for (k=1; k <= natom; k++)
869     {
870     if (nset[k] == j)
871     {
872     jset = j;
873     goto L_20;
874     }
875     }
876     return FALSE;
877     L_20:
878     if ( fabs(proi+aset[k][1]) > TOLER || fabs(ani-aset[k][0]) > TOL2 ||
879     fabs(disi-aset[k][2]) > TOLER )
880     goto L_30;
881    
882     for (k=1; k <= natom; k++)
883     {
884     if (nset[k] == jset)
885     {
886     x = (xi+a[k][0])*0.5;
887     y = (yi+a[k][1])*0.5;
888     theta = halfpi;
889     if (fabs(y) > TOL2) theta = -atan(x/y);
890     rotates(a,b,3,theta);
891     rotates(b,aset,2,pi);
892     iz = tcomp(b,aset);
893     if (iz == TRUE)
894     {
895     tcopy(b,a,natom);
896     return TRUE;
897     }
898     }
899     }
900     L_30:
901     continue;
902     }
903     return FALSE;
904     }
905     int findv(double **a, double **b, double **aset, int *nset)
906     {
907     int i,j, j1, numset, iset, jset, iz;
908     double halfpi, pi, x,y,theta;
909    
910     halfpi = 2.0*atan(1.0);
911     pi = 2.0*halfpi;
912    
913     circset(a,aset,3,nset,&numset);
914     iset = 1;
915     for (i=1; i < natom; i++)
916     {
917     if (nset[i] == iset)
918     {
919     jset = i;
920     goto L_5;
921     }
922     }
923     return FALSE;
924     L_5:
925     j1 = jset + 1;
926     for (j=j1; j <= natom; j++)
927     {
928     if (nset[j] == iset)
929     {
930     x = (a[i][0]+a[j][0])*0.5;
931     y = (a[i][1]+a[j][1])*0.5;
932     if (fabs(x) <= TOLER && fabs(y) <= TOLER)
933     {
934     x = 0.5*a[j][0];
935     y = 0.5*a[j][1];
936     }
937     theta = halfpi;
938     if (fabs(y) > TOLER) theta = -atan(x/y);
939     rotates(a,b,3,theta);
940     reflect(b,aset,1);
941     iz = tcomp(b,aset);
942     if (iz == TRUE)
943     {
944     tcopy(b,a,natom);
945     return TRUE;
946     }
947     }
948     }
949     return FALSE;
950     }
951     /* ======================================================= */
952     // compute moments of inertia
953     void secmom(double **xyz, double *xi,double *yi, double *zi, int *ishape,double t[3][3])
954     {
955     int i,j,k,ii;
956     float b[3];
957     double a[6], eigval[3], eigvec[3][3], atemp;
958     double summ, wti, xx,yy,zz;
959    
960     for (i=0; i < 3; i++)
961     {
962     b[i] = 0.0;
963     eigval[i] = 0.0;
964     for (j=0; j < 3; j++)
965     eigvec[i][j] = 0.0;
966     }
967     for (i=0; i < 6; i++)
968     a[i] = 0.0;
969    
970     summ = 0.0;
971     for (i=1; i <= natom; i++)
972     {
973     wti = atom[i].atomwt;
974     for (j=0; j < 3; j++)
975     b[j] += xyz[i][j]*wti;
976     summ += wti;
977     }
978     for (i=0; i < 3; i++)
979     b[i] /= summ;
980     for (i=1; i <= natom; i++)
981     {
982     for (j=0; j < 3; j++)
983     xyz[i][j] -= b[j];
984     }
985     // compute moments of inertia , sort and place smallest moment on x axis
986     for (i=1; i <= natom; i++)
987     {
988     wti = atom[i].atomwt;
989     *xi = xyz[i][0];
990     *yi = xyz[i][1];
991     *zi = xyz[i][2];
992     xx = *xi* *xi;
993     yy = *yi* *yi;
994     zz = *zi* *zi;
995     a[0] += wti*(yy+zz);
996     a[2] += wti*(zz+xx);
997     a[5] += wti*(xx+yy);
998     a[1] -= wti**xi**yi;
999     a[3] -= wti**xi**zi;
1000     a[4] -= wti**yi**zi;
1001     }
1002     summ = fabs(a[1]) + fabs(a[3]) + fabs(a[4]);
1003     if (summ <= 1.0e-15)
1004     {
1005     eigval[0] = a[0];
1006     eigval[1] = a[2];
1007     eigval[2] = a[5];
1008     eigvec[0][0] = 1.0;
1009     eigvec[1][1] = 1.0;
1010     eigvec[2][2] = 1.0;
1011     } else
1012     {
1013     hqrii(a,3,3,eigval,eigvec);
1014     }
1015    
1016     for (i=0; i < 2; i++)
1017     {
1018     ii = i + 1;
1019     for (j=ii; j < 3; j++)
1020     {
1021     if (eigval[i] > eigval[j])
1022     {
1023     atemp = eigval[i];
1024     eigval[i] = eigval[j];
1025     eigval[j] = atemp;
1026     for (k=0; k < 3; k++)
1027     {
1028     atemp = eigvec[k][i];
1029     eigvec[k][i] = eigvec[k][j];
1030     eigvec[k][j] = atemp;
1031     }
1032     }
1033     }
1034     }
1035     *xi = eigval[0] / 6.0228;
1036     *yi = eigval[1] / 6.0228;
1037     *zi = eigval[2] / 6.0228;
1038     // now rotate the molecule
1039     for (i=1; i <= natom; i++)
1040     {
1041     for (j=0; j < 3; j++)
1042     {
1043     b[j] = 0.0;
1044     for (k=0; k < 3; k++)
1045     b[j] += eigvec[j][k]*xyz[i][k];
1046     }
1047     xyz[i][0] = b[0];
1048     xyz[i][1] = b[1];
1049     xyz[i][2] = b[2];
1050     }
1051     // assign shape
1052     *ishape = 0;
1053     if (*xi <= 0.01)
1054     *ishape = 1;
1055     else
1056     {
1057     if (fabs(*xi+*yi-*zi) <= 0.01)
1058     *ishape = 6;
1059     else
1060     {
1061     if (fabs( (*zi/(*yi))-(*yi/(*xi))) <= 0.01)
1062     *ishape = 5;
1063     else
1064     {
1065     atemp = (2.0/(*yi) - 1.0/(*zi) - 1.0/(*xi))/ (1.0/(*xi) - 1.0/(*zi));
1066     if (atemp <= -0.999)
1067     *ishape = 3;
1068     else
1069     {
1070     if (atemp >= 0.999)
1071     *ishape = 4;
1072     else
1073     *ishape = 2;
1074     }
1075     }
1076     }
1077     }
1078     }
1079     // assign point group, symmetry number and chirality
1080     void vibsym(int mode)
1081     {
1082     int i;
1083     double **a;
1084     char ngrp[4];
1085    
1086     a = dmatrix (0, natom+1, 0,3);
1087     for (i=1; i <= natom; i++)
1088     {
1089     a[i][0] = atom[i].x;
1090     a[i][1] = atom[i].y;
1091     a[i][2] = atom[i].z;
1092     }
1093    
1094     if (mode == 0)
1095     {
1096     strcpy(ngrp,"");
1097     ptgrp(mode, ngrp,a);
1098     }
1099     free_dmatrix(a ,0, natom+1, 0,3);
1100     }
1101    
1102     void ptgrp(int mode,char *ngrp,double **a)
1103     {
1104     int i, j, k, ishape, numsym, iaxis, maxcn, xxx;
1105     int iz, iz1, iz2, iopt;
1106     int chiral, *nset;
1107     double **b, **c;
1108     double t[3][3], p1[3], p2[3];
1109     double xi, yi, zi;
1110     float halfpi, pi;
1111     char string[256];
1112    
1113     halfpi = 2.0*atan(1.0);
1114     pi = 2.0*halfpi;
1115     chiral = FALSE;
1116     strcpy(ngrp,"Cs");
1117    
1118    
1119     b = dmatrix (0, natom+1, 0,3);
1120     c = dmatrix (0, natom+1, 0,3);
1121     nset = ivector(0,natom+1);
1122    
1123     secmom(a, &xi, &yi, &zi, &ishape, t);
1124    
1125     if (ishape == 1) // linear
1126     {
1127     p1[0] = 0.0;
1128     p1[1] = 0.0;
1129     p1[2] = 0.0;
1130     p2[0] = 1.0;
1131     p2[1] = 0.0;
1132     p2[1] = 0.0;
1133     orients(a,b,p1,p2,3,natom);
1134     tcopy(b,a,natom);
1135     iaxis = findh(a,b,3);
1136     if (iaxis == TRUE)
1137     {
1138     strcpy(ngrp,"D*h");
1139     numsym = 2;
1140     } else
1141     {
1142     strcpy(ngrp,"C*v");
1143     numsym = 1;
1144     }
1145     } else if (ishape == 2) // asymmetric
1146     /*------------------------*
1147     ASyMMETRiC TOP MOLECULES
1148     *------------------------*
1149    
1150     These molecules can have no axes of order greater than 2. Thus
1151     the possible point groups are: D2H, D2, C2V, C2H, C2, Ci, CS,
1152     and C1. FiNDCN routine finds the principal axes and align the
1153     principal axes with the z-axis. */
1154     {
1155     iaxis = findcn(a,b,p1,p2,&maxcn);
1156     if (maxcn == 1)
1157     {
1158     j = findh(a,b,3);
1159     if (j == TRUE)
1160     {
1161     strcpy(ngrp,"Cs");
1162     numsym = 1;
1163     }else
1164     {
1165     i = findh(a,b,2);
1166     if (i == TRUE)
1167     {
1168     strcpy(ngrp,"Cs");
1169     numsym = 1;
1170     } else
1171     {
1172     k = findh(a,b,1);
1173     if (k == TRUE)
1174     {
1175     strcpy(ngrp,"Cs");
1176     numsym = 1;
1177     } else
1178     {
1179     xxx = findi(a,b);
1180     if (xxx == TRUE)
1181     {
1182     strcpy(ngrp,"Ci");
1183     numsym = 1;
1184     } else
1185     {
1186     strcpy(ngrp,"C1");
1187     numsym = 1;
1188     chiral = TRUE;
1189     }
1190     }
1191     }
1192     }
1193     } else // maxcn == 2
1194     {
1195     iaxis = findc2(a,b,c,nset);
1196     if (iaxis == FALSE)
1197     {
1198     i = findv(a,b,c,nset);
1199     if (i == TRUE)
1200     {
1201     strcpy(ngrp,"C2v");
1202     numsym = 2;
1203     } else
1204     {
1205     j = findh(a,b,3);
1206     if (j == TRUE)
1207     {
1208     strcpy(ngrp,"C2h");
1209     numsym = 2;
1210     } else
1211     {
1212     strcpy(ngrp,"C2");
1213     numsym = 2;
1214     chiral = TRUE;
1215     }
1216     }
1217     } else
1218     {
1219     i = findh(a,b,3);
1220     if (i == TRUE)
1221     {
1222     strcpy(ngrp,"D2h");
1223     numsym = 4;
1224     } else
1225     {
1226     j = findv(a,b,c,nset);
1227     if (j == TRUE)
1228     {
1229     strcpy(ngrp,"D2d");
1230     numsym = 4;
1231     } else
1232     {
1233     strcpy(ngrp,"D2");
1234     numsym = 4;
1235     chiral = TRUE;
1236     }
1237     }
1238     }
1239     }
1240     } else if (ishape == 3 || ishape == 4) //prolate
1241     {
1242     /*-----------------------*
1243     SyMMETRiC TOP MOLECULES
1244     *-----------------------*
1245    
1246     These molecules can belong to any axial point group, thus only
1247     the cubic point groups (T, Td, Th, O, Oh, i, ih) are impossible.
1248     However, except in rare cases the unique axis is a rotation axis
1249     of order 3 or greater. */
1250    
1251     iz = findcn(a,b,p1,p2,&maxcn);
1252     iz = findsn(a,b,c,&maxcn);
1253     if (iz == TRUE)
1254     {
1255     iz1 = findv(a,b,c,nset);
1256     if (iz1 == FALSE)
1257     {
1258     sprintf(ngrp,"S%1d",2*maxcn);
1259     numsym = maxcn;
1260     goto L_2000;
1261     }
1262     }
1263     iz = findc2(a,b,c,nset);
1264     if (iz == FALSE)
1265     {
1266     iz1 = findv(a,b,c,nset);
1267     if (iz1 == TRUE)
1268     {
1269     sprintf(ngrp,"C%1dv",maxcn);
1270     numsym = maxcn;
1271     } else
1272     {
1273     iz2 = findh(a,b,3);
1274     if (iz2 == TRUE)
1275     {
1276     sprintf(ngrp,"C%1dh",maxcn);
1277     numsym = maxcn;
1278     } else
1279     {
1280     sprintf(ngrp,"C%1d",maxcn);
1281     numsym = maxcn;
1282     chiral = TRUE;
1283     }
1284     }
1285     }else
1286     {
1287     iz1 = findh(a,b,3);
1288     if (iz1 == TRUE)
1289     {
1290     sprintf(ngrp,"D%1dh",maxcn);
1291     numsym = 2*maxcn;
1292     } else
1293     {
1294     iz2 = findv(a,b,c,nset);
1295     if (iz2 == TRUE)
1296     {
1297     sprintf(ngrp,"D%1dd",maxcn);
1298     numsym = 2*maxcn;
1299     } else
1300     {
1301     sprintf(ngrp,"D%1d",maxcn);
1302     numsym = 2*maxcn;
1303     chiral = TRUE;
1304     }
1305     }
1306     }
1307     // SPHERiCAL TOP MOLECULES
1308     } else if (ishape == 5) // spherical
1309     {
1310     sphere(a,b,c,nset,&maxcn);
1311     iopt = maxcn - 2;
1312     if (iopt == 1)
1313     {
1314     iz = findi(a,b);
1315     if (iz == TRUE)
1316     {
1317     sprintf(ngrp,"Th");
1318     numsym = 12;
1319     } else
1320     {
1321     iz1 = findv(a,b,c,nset);
1322     if (iz1 == TRUE)
1323     {
1324     sprintf(ngrp,"Td");
1325     numsym = 12;
1326     } else
1327     {
1328     sprintf(ngrp,"T");
1329     numsym = 12;
1330     chiral = TRUE;
1331     }
1332     }
1333     } else if (iopt == 2)
1334     {
1335     iz = findi(a,b);
1336     if (iz == TRUE)
1337     {
1338     sprintf(ngrp,"Oh");
1339     numsym = 24;
1340     } else
1341     {
1342     sprintf(ngrp,"O");
1343     numsym = 24;
1344     chiral = TRUE;
1345     }
1346     } else if (iopt == 3)
1347     {
1348     iz = findi(a,b);
1349     if (iz == TRUE)
1350     {
1351     sprintf(ngrp,"Ih");
1352     numsym = 60;
1353     } else
1354     {
1355     sprintf(ngrp,"I");
1356     numsym = 60;
1357     chiral = TRUE;
1358     }
1359     }
1360     } else if (ishape == 6) // planar
1361     {
1362     iz = findh(a,b,1);
1363     iz1 = findh(a,b,2);
1364     iz2 = findcn(a,b,p1,p2,&maxcn);
1365     if (maxcn == 2)
1366     {
1367     if (iz == TRUE && iz1 == TRUE)
1368     {
1369     sprintf(ngrp,"D%1dh",maxcn);
1370     numsym = maxcn*2;
1371     } else
1372     {
1373     if (iz == TRUE || iz1 == TRUE)
1374     {
1375     sprintf(ngrp,"C%1dv",maxcn);
1376     numsym = maxcn;
1377     }else
1378     {
1379     i = findh(a,b,3);
1380     if (i == TRUE)
1381     {
1382     sprintf(ngrp,"C%1dh",maxcn);
1383     numsym = maxcn;
1384     } else
1385     {
1386     sprintf(ngrp,"C%1d",maxcn);
1387     numsym = maxcn;
1388     chiral = TRUE;
1389     }
1390     }
1391     }
1392     } else
1393     {
1394     iz = findsn(a,b,c,&maxcn);
1395     if (iz == TRUE)
1396     {
1397     iz1 = findv(a,b,c,nset);
1398     if (iz1 == FALSE)
1399     {
1400     sprintf(ngrp,"S%1d",2*maxcn);
1401     numsym = maxcn;
1402     goto L_2000;
1403     }
1404     }
1405     iz = findc2(a,b,c,nset);
1406     if (iz == FALSE)
1407     {
1408     iz1 = findv(a,b,c,nset);
1409     if (iz1 == TRUE)
1410     {
1411     sprintf(ngrp,"C%1dv",maxcn);
1412     numsym = maxcn;
1413     } else
1414     {
1415     iz2 = findh(a,b,3);
1416     if (iz2 == TRUE)
1417     {
1418     sprintf(ngrp,"C%1dh",maxcn);
1419     numsym = maxcn;
1420     } else
1421     {
1422     sprintf(ngrp,"C%1d",maxcn);
1423     numsym = maxcn;
1424     chiral = TRUE;
1425     }
1426     }
1427     } else
1428     {
1429     iz1 = findh(a,b,3);
1430     if (iz1 == TRUE)
1431     {
1432     sprintf(ngrp,"D%1dh",maxcn);
1433     numsym = 2*maxcn;
1434     } else
1435     {
1436     iz2 = findv(a,b,c,nset);
1437     if (iz2 == TRUE)
1438     {
1439     sprintf(ngrp,"D%1dd",maxcn);
1440     numsym = 2*maxcn;
1441     } else
1442     {
1443     sprintf(ngrp,"D%1d",maxcn);
1444     numsym = 2*maxcn;
1445     chiral = TRUE;
1446     }
1447     }
1448     }
1449     }
1450     }
1451     L_2000:
1452     free_dmatrix(b ,0, natom+1, 0,3);
1453     free_dmatrix(c ,0, natom+1, 0,3);
1454     free_ivector(nset, 0, natom+1);
1455    
1456     if (mode == 2)
1457     {
1458     symmetry.chiral = chiral;
1459     symmetry.numsym = numsym;
1460     symmetry.ishape = ishape;
1461     symmetry.xi = xi;
1462     symmetry.yi = yi;
1463     symmetry.zi = zi;
1464     }
1465    
1466     if (mode == 0)
1467     {
1468     sprintf(string,"Symmetry group: %s\nMoments of Inertia\nIx: %7.2f\nIy: %7.2f\nIz: %7.2f",ngrp,xi,yi,zi);
1469     message_alert(string,"Symmetry Information");
1470     }
1471     }
1472     /* ========================================== */
1473     void hqrii(double *a,int n,int m,double *e,double v[3][3])
1474     {
1475     #define V(I_,J_) (*(v+(I_)*(n)+(J_)))
1476     int _do0, _do1, _do10, _do11, _do12, _do13, _do14, _do15,
1477     _do16, _do17, _do18, _do19, _do2, _do20, _do21, _do22, _do3,
1478     _do4, _do5, _do6, _do7, _do8, _do9, i, i_, ig, ii, im1, iord,
1479     ip1, irank, itere, itere_, j, j_, jrank, k, k_, kk, kk_, kp1,
1480     kpiv, krank, l, ll, nm1, nm2;
1481     double c, del, ee, eps, eps1, eps2, eps3, ff, fn, gersch, h, r,
1482     ra, rn, s, seps, sinv, sorter, sum, summ, t, u, vn, w[5][5],
1483     ww, z, zero;
1484    
1485     double *const A = &a[0] - 1;
1486     double *const E = &e[0] - 1;
1487    
1488     /*************************************************************
1489     *
1490     * HQRII IS A DIAGONALISATION ROUTINE, WRITTEN BY YOSHITAKA BEPPU OF
1491     * NAGOYA UNIVERSITY, JAPAN.
1492     * FOR DETAILS SEE 'COMPUTERS & CHEMISTRY' VOL.6 1982. PAGE 000.
1493     *
1494     * ON INPUT A = MATRIX TO BE DIAGONALISED (PACKED CANONICAL)
1495     * N = SIZE OF MATRIX TO BE DIAGONALISED.
1496     * M = NUMBER OF EIGENVECTORS NEEDED.
1497     * E = ARRAY OF SIZE AT LEAST N
1498     * V = ARRAY OF SIZE AT LEAST NMAX*M
1499     *
1500     * ON OUTPUT E = EIGENVALUES
1501     * V = EIGENVECTORS IN ARRAY OF SIZE NMAX*M
1502     *
1503     ************************************************************************ */
1504    
1505     /* EPS3 AND EPS ARE MACHINE-PRECISION DEPENDENT
1506     * */
1507     eps3 = 1.e-30;
1508     zero = 0.e0;
1509     ll = (n*(n + 1))/2 + 1;
1510     eps = 1.e-8;
1511     iord = -1;
1512     nm1 = n - 1;
1513     if( n == 2 )
1514     goto L_90;
1515     nm2 = n - 2;
1516     krank = 0;
1517     /* HOUSEHOLDER TRANSFORMATION */
1518     for( k = 1, _do0 = nm2; k <= _do0; k++ )
1519     {
1520     k_ = k - 1;
1521     kp1 = k + 1;
1522     krank += k;
1523     w[k_][1] = A[krank];
1524     sum = 0.;
1525     jrank = krank;
1526     for( j = kp1, _do1 = n; j <= _do1; j++ )
1527     {
1528     j_ = j - 1;
1529     w[j_][1] = A[jrank + k];
1530     jrank += j;
1531     sum += (w[j_][1])*(w[j_][1]);
1532     }
1533     s = sign( sqrt( sum ), w[kp1 - 1][1] );
1534     w[k_][0] = -s;
1535     w[kp1 - 1][1] += s;
1536     A[k + krank] = w[kp1 - 1][1];
1537     h = w[kp1 - 1][1]*s;
1538     if( fabs( h ) < eps3 )
1539     goto L_80;
1540     summ = 0.e0;
1541     irank = krank;
1542     for( i = kp1, _do2 = n; i <= _do2; i++ )
1543     {
1544     i_ = i - 1;
1545     sum = 0.e0;
1546     for( j = kp1, _do3 = i; j <= _do3; j++ )
1547     {
1548     j_ = j - 1;
1549     sum += A[j + irank]*w[j_][1];
1550     }
1551     if( i >= n )
1552     goto L_40;
1553     ip1 = i + 1;
1554     jrank = i*(i + 3)/2;
1555     for( j = ip1, _do4 = n; j <= _do4; j++ )
1556     {
1557     j_ = j - 1;
1558     sum += A[jrank]*w[j_][1];
1559     jrank += j;
1560     }
1561     L_40:
1562     w[i_][0] = sum/h;
1563     irank += i;
1564     summ += w[i_][0]*w[i_][1];
1565     }
1566     u = summ*0.5e0/h;
1567     jrank = krank;
1568     for( j = kp1, _do5 = n; j <= _do5; j++ )
1569     {
1570     j_ = j - 1;
1571     w[j_][0] = w[j_][1]*u - w[j_][0];
1572     for( i = kp1, _do6 = j; i <= _do6; i++ )
1573     {
1574     i_ = i - 1;
1575     A[i + jrank] += w[i_][0]*w[j_][1] + w[j_][0]*w[i_][1];
1576     }
1577     jrank += j;
1578     }
1579     L_80:
1580     A[krank] = h;
1581     }
1582     L_90:
1583     w[nm1 - 1][1] = A[(nm1*(nm1 + 1))/2];
1584     w[n - 1][1] = A[(n*(n + 1))/2];
1585     w[nm1 - 1][0] = A[nm1 + (n*(n - 1))/2];
1586     w[n - 1][0] = 0.e0;
1587     gersch = fabs( w[0][1] ) + fabs( w[0][0] );
1588     for( i = 1, _do7 = nm1; i <= _do7; i++ )
1589     {
1590     i_ = i - 1;
1591     if ( fabs( w[i_ + 1][1] ) + fabs( w[i_][0] ) + fabs( w[i_ + 1][0] ) > gersch)
1592     gersch = fabs( w[i_ + 1][1] ) + fabs( w[i_][0] ) + fabs( w[i_ + 1][0] );
1593     }
1594     del = eps*gersch;
1595     for( i = 1, _do8 = n; i <= _do8; i++ )
1596     {
1597     i_ = i - 1;
1598     w[i_][2] = w[i_][0];
1599     E[i] = w[i_][1];
1600     v[m-1][i_] = E[i];
1601     }
1602     if( fabs( del ) < eps3 )
1603     goto L_220;
1604     /* QR-METHOD WITH ORIGIN SHIFT */
1605     k = n;
1606     L_120:
1607     l = k;
1608     L_130:
1609     if( fabs( w[l - 2][2] ) < del )
1610     goto L_140;
1611     l -= 1;
1612     if( l > 1 )
1613     goto L_130;
1614     L_140:
1615     if( l == k )
1616     goto L_170;
1617     ww = (E[k - 1] + E[k])*0.5e0;
1618     r = E[k] - ww;
1619     z = sign( sqrt( (w[k - 2][2])*(w[k - 2][2]) + r*r ), r ) + ww;
1620     ee = E[l] - z;
1621     E[l] = ee;
1622     ff = w[l - 1][2];
1623     r = sqrt( ee*ee + ff*ff );
1624     j = l;
1625     goto L_160;
1626     L_150:
1627     r = sqrt( (E[j])*(E[j]) + (w[j - 1][2])*(w[j - 1][2]) );
1628     w[j - 2][2] = s*r;
1629     ee = E[j]*c;
1630     ff = w[j - 1][2]*c;
1631     L_160:
1632     c = E[j]/r;
1633     s = w[j - 1][2]/r;
1634     ww = E[j + 1] - z;
1635     E[j] = (ff*c + ww*s)*s + ee + z;
1636     E[j + 1] = c*ww - s*ff;
1637     j += 1;
1638     if( j < k )
1639     goto L_150;
1640     w[k - 2][2] = E[k]*s;
1641     E[k] = E[k]*c + z;
1642     goto L_120;
1643     L_170:
1644     k -= 1;
1645     if( k > 1 )
1646     goto L_120;
1647     /* * * * * * * * * * * * *
1648     *
1649     * AT THIS POINT THE ARRAY 'E' CONTAINS THE UN-ORDERED EIGENVALUES
1650     *
1651     * * * * * * * * * * * * *
1652     * STRAIGHT SELECTION SORT OF EIGENVALUES */
1653     sorter = 1.e0;
1654     if( iord < 0 )
1655     sorter = -1.e0;
1656     j = n;
1657     L_180:
1658     l = 1;
1659     ii = 1;
1660     ll = 1;
1661     for( i = 2, _do9 = j; i <= _do9; i++ )
1662     {
1663     i_ = i - 1;
1664     if( (E[i] - E[l])*sorter > 0.e0 )
1665     goto L_190;
1666     l = i;
1667     goto L_200;
1668     L_190:
1669     ii = i;
1670     ll = l;
1671     L_200:
1672     ;
1673     }
1674     if( ii == ll )
1675     goto L_210;
1676     ww = E[ll];
1677     E[ll] = E[ii];
1678     E[ii] = ww;
1679     L_210:
1680     j = ii - 1;
1681     if( j >= 2 )
1682     goto L_180;
1683     L_220:
1684     if( !m )
1685     return;
1686     /***************
1687     * ORDERING OF EIGENVALUES COMPLETE.
1688     ***************
1689     * INVERSE-ITERATION FOR EIGENVECTORS */
1690     fn = (float)( n );
1691     eps1 = 1.e-5;
1692     seps = sqrt( eps );
1693     eps2 = 0.05e0;
1694     rn = 0.e0;
1695     ra = eps*0.6180339887485e0;
1696     /* 0.618... IS THE FIBONACCI NUMBER (-1+SQRT(5))/2. */
1697     ig = 1;
1698     for( i = 1, _do10 = m; i <= _do10; i++ )
1699     {
1700     i_ = i - 1;
1701     im1 = i - 1;
1702     for( j = 1, _do11 = n; j <= _do11; j++ )
1703     {
1704     j_ = j - 1;
1705     w[j_][2] = 0.e0;
1706     w[j_][3] = w[j_][0];
1707     w[j_][4] = v[m-1][j_] - E[i];
1708     rn += ra;
1709     if( rn >= eps )
1710     rn -= eps;
1711     v[i_][j_] = rn;
1712     }
1713     for( j = 1, _do12 = nm1; j <= _do12; j++ )
1714     {
1715     j_ = j - 1;
1716     if( fabs( w[j_][4] ) >= fabs( w[j_][0] ) )
1717     goto L_240;
1718     w[j_][1] = -w[j_][4]/w[j_][0];
1719     w[j_][4] = w[j_][0];
1720     t = w[j_ + 1][4];
1721     w[j_ + 1][4] = w[j_][3];
1722     w[j_][3] = t;
1723     w[j_][2] = w[j_ + 1][3];
1724     if( fabs( w[j_][2] ) < eps3 )
1725     w[j_][2] = del;
1726     w[j_ + 1][3] = 0.e0;
1727     goto L_250;
1728     L_240:
1729     if( fabs( w[j_][4] ) < eps3 )
1730     w[j_][4] = del;
1731     w[j_][1] = -w[j_][0]/w[j_][4];
1732     L_250:
1733     w[j_ + 1][3] += w[j_][2]*w[j_][1];
1734     w[j_ + 1][4] += w[j_][3]*w[j_][1];
1735     }
1736     if( fabs( w[n - 1][4] ) < eps3 )
1737     w[n - 1][4] = del;
1738     for( itere = 1; itere <= 5; itere++ )
1739     {
1740     itere_ = itere - 1;
1741     if( itere == 1 )
1742     goto L_280;
1743     for( j = 1, _do13 = nm1; j <= _do13; j++ )
1744     {
1745     j_ = j - 1;
1746     if( fabs( w[j_][2] ) < eps3 )
1747     goto L_270;
1748     t = v[i_][j_];
1749     v[i_][j_] = v[i_][j_+1];
1750     v[i_][j_+1] = t;
1751     L_270:
1752     v[i_][j_+1] += v[i_][j_]*w[j_][1];
1753     }
1754     L_280:
1755     v[i_][n-1] /= w[n - 1][4];
1756     v[i_][nm1 - 1] = (v[i_][nm1 - 1] - v[i_][n-1]*w[nm1 - 1][3])/w[nm1 - 1][4];
1757     vn = 1.0e-20;
1758     if (fabs( v[i_][n-1] ) > vn)
1759     vn = fabs( v[i_][n-1] );
1760     if ( fabs( v[i_][nm1 - 1] ) > vn)
1761     vn = fabs( v[i_][nm1 - 1] );
1762     if( n == 2 )
1763     goto L_300;
1764     k = nm2;
1765     L_290:
1766     v[i_][k - 1] = (v[i_][k - 1] - v[i_][k]*w[k - 1][3] - v[i_][k + 1]*
1767     w[k - 1][2])/w[k - 1][4];
1768     if (vn < 1.0e-20)
1769     vn = 1.0e-20;
1770     if (fabs( v[i_][k - 1] ) > vn)
1771     vn = fabs( v[i_][n-1] );
1772     k -= 1;
1773     if( k >= 1 )
1774     goto L_290;
1775     L_300:
1776     s = eps1/vn;
1777     for( j = 1, _do14 = n; j <= _do14; j++ )
1778     {
1779     j_ = j - 1;
1780     v[i_][j_] *= s;
1781     }
1782     if( itere > 1 && vn > 1 )
1783     goto L_330;
1784     }
1785     /* TRANSFORMATION OF EIGENVECTORS */
1786     L_330:
1787     if( n == 2 )
1788     goto L_380;
1789     krank = nm2*(n + 1)/2;
1790     kpiv = nm2*nm1/2;
1791     for( k = nm2; k >= 1; k-- )
1792     {
1793     k_ = k - 1;
1794     kp1 = k + 1;
1795     if( fabs( A[kpiv] ) <= eps3 )
1796     goto L_360;
1797     sum = 0.e0;
1798     for( kk = kp1, _do15 = n; kk <= _do15; kk++ )
1799     {
1800     kk_ = kk - 1;
1801     sum += A[krank]*v[i_][kk_];
1802     krank += kk;
1803     }
1804     s = -sum/A[kpiv];
1805     for( kk = n, _do16 = kp1; kk >= _do16; kk-- )
1806     {
1807     kk_ = kk - 1;
1808     krank -= kk;
1809     v[i_][kk_] += A[krank]*s;
1810     }
1811     L_360:
1812     kpiv -= k;
1813     krank -= kp1;
1814     }
1815     L_380:
1816     for( j = ig, _do17 = i; j <= _do17; j++ )
1817     {
1818     j_ = j - 1;
1819     if( fabs( E[j] - E[i] ) < eps2 )
1820     goto L_400;
1821     }
1822     j = i;
1823     L_400:
1824     ig = j;
1825     if( ig == i )
1826     goto L_430;
1827     /* RE-ORTHOGONALISATION */
1828     for( k = ig, _do18 = im1; k <= _do18; k++ )
1829     {
1830     k_ = k - 1;
1831     sum = 0.e0;
1832     for( j = 1, _do19 = n; j <= _do19; j++ )
1833     {
1834     j_ = j - 1;
1835     sum += v[k_][j_]*v[i_][j_];
1836     }
1837     s = -sum;
1838     for( j = 1, _do20 = n; j <= _do20; j++ )
1839     {
1840     j_ = j - 1;
1841     v[i_][j_] += v[k_][j_]*s;
1842     }
1843     }
1844     /* NORMALISATION */
1845     L_430:
1846     sum = 1.e-24;
1847     for( j = 1, _do21 = n; j <= _do21; j++ )
1848     {
1849     j_ = j - 1;
1850     sum += (v[i_][j_])*(v[i_][j_]);
1851     }
1852     sinv = 1.e0/sqrt( sum );
1853     for( j = 1, _do22 = n; j <= _do22; j++ )
1854     {
1855     j_ = j - 1;
1856     v[i_][j_] *= sinv;
1857     }
1858     }
1859     return;
1860     #undef V
1861     }