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