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