ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/asnsym.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 6 months ago) by tjod
File size: 53717 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
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 }