ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/vibrate.c
Revision: 110
Committed: Thu Mar 12 01:43:50 2009 UTC (12 years, 1 month ago) by gilbertke
File size: 37467 byte(s)
Log Message:
further cleanup and localization of atom data
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4
5 #include "asnsym.h"
6 #include "utility.h"
7 #include "bonds_ff.h"
8 #include "vibrate.h"
9
10 char *symname[] = {"A","A'","A''","Ag","Au",
11 "A1","A2","A1'","A2'","A1''","A2''",
12 "A1g","A1u","A2g","A2u","B",
13 "B1","B2","B3","Bg","Bu","B1g",
14 "B1u","B2g","B2u","B3g","B3u","E",
15 "E1","E2","E3","E'","E''",
16 "Eg","Eu","E1'","E1''","E2'",
17 "E2''","E1g","E1u","E2g","E2u","T","Tg","Tu" };
18
19
20 EXTERN double hesscut;
21
22 // FILE * fopen_path ( char * , char * , char * ) ;
23
24 static void tred2(double **,int , double *, double *);
25 static void tqli(double *,double *,int ,double **);
26 static void assign_vibsym(int, char *,char *,int *,double **);
27 double get_total_energy(void);
28 void hessian(int, double *,int *, int *, int *,double *);
29 void vibrate(int natom,double *atomwt,double *x,double *y,double *z,double *charge);
30 //void vibcharge_dipole(double **);
31 //void vibdipole_dipole(double **);
32
33 void vibrate(int natom,double *atomwt,double *x,double *y,double *z,double *charge)
34 {
35 int i, k, j, ii, jj, nsym,imix,mm;
36 int maxvar, maxhess, ihess, nfreq, maxvib;
37 int *hinit, *hstop, *hindex;
38 int *vibsym;
39 double mass, efact,cnosym,cimix,temp,qptf,wdt,efwt,cpfact,efwt2;
40 double xyzi,ssym,entrov1,entrov2,ezpe;
41 double etrans,htrans,strans,gtrans,cptrans;
42 double erot,hrot,srot,grot,cprot;
43 double evib,hvib,svib,gvib,cpvib;
44 double emix,hmix,smix,gmix,cpmix;
45 double etot,htot,stot,gtot,cptot;
46 double *hdiag, *h;
47 double **matrix, *a;
48 double **acoord;
49 double *eigen, **vects, *mass2;
50 double factor, convert, lightspd, vnorm;
51 double *dipole;
52 double dqx,dqy,dqz,dqxx,dqyy,dqzz,dqxy,dqxz,dqyx,dqyz,dqzx,dqzy;
53 char ngrp[4], vsym[4];
54
55 convert = 4.184e2;
56 lightspd = 2.99792458e-2;
57 factor = sqrt(convert)/(2.0*PI*lightspd);
58
59 for (i=1; i <= natom; i++)
60 {
61 if (atomwt[i] <= 0.0)
62 atomwt[i] = 0.001;
63 }
64
65 maxvar = 3*natom;
66 maxhess = (3*natom*(3*natom-1))/2;
67 maxvib = 3*natom;
68
69 vibsym = ivector(0,maxvar);
70 hinit = ivector(0,maxvar);
71 hstop = ivector(0,maxvar);
72 hindex = ivector(0,maxhess);
73 hdiag = dvector(0,maxvar);
74 h = dvector(0,maxhess);
75 matrix = dmatrix(0, maxvib, 0, maxvib);
76 a = dvector(0, maxvib);
77 mass2 = dvector(0, natom+1);
78 eigen = dvector(0, 3*natom);
79 vects = dmatrix(0, 3*natom, 0, 3*natom);
80 acoord = dmatrix (0, natom+1, 0,3);
81 dipole = dvector(0,maxvar);
82
83 for (i=1; i <= natom; i++)
84 mass2[i] = sqrt(atomwt[i]);
85
86 hesscut = 0.0;
87 hessian(maxhess, h, hinit, hstop, hindex, hdiag);
88 ihess = 0;
89 for (i=0; i < maxvar; i++)
90 {
91 matrix[i][i] = hdiag[i];
92 ii = i;
93 for (k = hinit[i]; k < hstop[i]; k++)
94 {
95 ii++;
96 matrix[i][ii] = h[k];
97 matrix[ii][i] = h[k];
98 }
99 }
100 nfreq = maxvar;
101
102 tred2(matrix, nfreq, eigen, a);
103 tqli(eigen,a,nfreq,matrix);
104
105 // now do mass weighted
106 ihess = 0;
107 for (i=0; i < maxvar; i++)
108 {
109 jj = i/3 + 1;
110 matrix[i][i] = hdiag[i]/atomwt[jj];
111 ii = i;
112 for (k = hinit[i]; k < hstop[i]; k++)
113 {
114 ii++;
115 j = ii/3 + 1;
116 matrix[i][ii] = h[k]/(mass2[jj]*mass2[j]);
117 matrix[ii][i] = matrix[i][ii];
118 }
119 }
120
121 nfreq = maxvar;
122 tred2(matrix, nfreq, eigen, a);
123 tqli(eigen,a,nfreq,matrix);
124
125 // vibfile = fopen_path(Savebox.path,Savebox.fname,"w");
126
127 /* if (vibfile == NULL)
128 {
129 message_alert("Error opening vibration output file.","Vibration Setup");
130 } */
131
132 // fprintf(vibfile,"\nCurrent Structure\n");
133
134 // for (i=1; i <= natom; i++)
135 // fprintf(vibfile,"%d %d %f %f %f\n",i,atom[i].atomnum,atom[i].x, atom[i].y, atom[i].z);
136
137
138 // fprintf(vibfile,"\n\nMass Weighted Eigenvectors\n");
139
140 for (i=0; i < maxvar; i++)
141 eigen[i] = factor*sqrt(fabs(eigen[i]));
142
143 for (i=0; i < maxvar; i++)
144 {
145 vnorm = 0.0;
146 for (j=0; j < maxvar; j++)
147 {
148 k = j/3 + 1;
149 matrix[j][i] /= mass2[k];
150 vnorm += matrix[j][i]*matrix[j][i];
151 }
152 vnorm = sqrt(vnorm);
153 for (j=0; j < maxvar; j++)
154 matrix[j][i] /= vnorm;
155 }
156 ///
157 // get overall molecular symmetry
158 for (i=1; i <= natom; i++)
159 {
160 acoord[i][0] = x[i];
161 acoord[i][1] = y[i];
162 acoord[i][2] = z[i];
163 }
164 strcpy(ngrp," ");
165 ptgrp(2,ngrp,acoord);
166 // compute dipole moment
167
168 for (i=0; i < maxvar; i++)
169 {
170 dqx = 0.0;
171 dqy = 0.0;
172 dqz = 0.0;
173 dqxx = dqyy = dqzz = 0.0;
174 dqxy = dqxz = 0.0;
175 dqyx = dqyz = 0.0;
176 dqzx = dqzy = 0.0;
177
178 for (j=0; j < maxvar; j+= 3)
179 {
180 strcpy(vsym," ");
181 k = j/3 + 1;
182 acoord[k][0] = x[k] + matrix[j][i];
183 acoord[k][1] = y[k] + matrix[j+1][i];
184 acoord[k][2] = z[k] + matrix[j+2][i];
185 dqx += charge[k]*matrix[j][i];
186 dqy += charge[k]*matrix[j+1][i];
187 dqz += charge[k]*matrix[j+2][i];
188
189 }
190
191 ptgrp(1,vsym,acoord);
192 assign_vibsym(natom,ngrp,vsym,&nsym,acoord);
193 vibsym[i] = nsym;
194
195 }
196
197 /* icycle = 3*natom/5;
198 for (i=0; i < (5*icycle); i += 5)
199 {
200 fprintf(vibfile,"\nFreq: %8.3f %8.3f %8.3f %8.3f %8.3f\n",eigen[i],eigen[i+1],eigen[i+2],eigen[i+3],eigen[i+4]);
201 fprintf(vibfile,"Sym: %s %s %s %s %s\n",symname[vibsym[i]],symname[vibsym[i+1]],symname[vibsym[i+2]],
202 symname[vibsym[i+3]],symname[vibsym[i+4]]);
203 fprintf(vibfile,"Int: %8.4f %8.4f %8.4f %8.4f %8.4f\n",dipole[i],dipole[i+1],dipole[i+2],dipole[i+3],dipole[i+4]);
204 fprintf(vibfile,"\n");
205
206 for (j=0; j < maxvar; j++)
207 fprintf(vibfile,"%-d : %8.3f %8.3f %8.3f %8.3f %8.3f\n",j,matrix[j][i],matrix[j][i+1],matrix[j][i+2],
208 matrix[j][i+3],matrix[j][i+4]);
209 }
210 if (icycle*5 < 3*natom)
211 {
212 iz = (3*natom-icycle*5);
213 i = icycle*5;
214 if (iz == 1)
215 {
216 fprintf(vibfile,"\nFreq: %8.3f\n",eigen[i]);
217 fprintf(vibfile,"Sym: %s\n",symname[vibsym[i]]);
218 fprintf(vibfile,"Int: %8.4f\n",dipole[i]);
219 fprintf(vibfile,"\n");
220
221 for (j=0; j < maxvar; j++)
222 fprintf(vibfile,"%-d : %8.3f %8.3f\n",j,matrix[j][i],matrix[j][i+1]);
223 } else if (iz == 2)
224 {
225 fprintf(vibfile,"\nFreq: %8.3f %8.3f\n",eigen[i],eigen[i+1]);
226 fprintf(vibfile,"Sym: %s %s\n",symname[vibsym[i]],symname[vibsym[i+1]]);
227 fprintf(vibfile,"Int: %8.4f %8.4f\n",dipole[i],dipole[i+1]);
228 fprintf(vibfile,"\n");
229
230 for (j=0; j < maxvar; j++)
231 fprintf(vibfile,"%-d : %8.3f %8.3f\n",j,matrix[j][i],matrix[j][i+1]);
232 } else if (iz == 3)
233 {
234 fprintf(vibfile,"\nFreq: %8.3f %8.3f %8.3f\n",eigen[i],eigen[i+1],eigen[i+2]);
235 fprintf(vibfile,"Sym: %s %s %s\n",symname[vibsym[i]],symname[vibsym[i+1]],symname[vibsym[i+2]]);
236 fprintf(vibfile,"Int: %8.4f %8.4f %8.4f\n",dipole[i],dipole[i+1],dipole[i+2]);
237 fprintf(vibfile,"\n");
238
239 for (j=0; j < maxvar; j++)
240 fprintf(vibfile,"%-d : %8.3f %8.3f %8.3f\n",j,matrix[j][i],matrix[j][i+1],matrix[j][i+2]);
241 } else if (iz == 4)
242 {
243 fprintf(vibfile,"\nFreq: %8.3f %8.3f %8.3f %8.3f\n",eigen[i],eigen[i+1],eigen[i+2],eigen[i+3]);
244 fprintf(vibfile,"Sym: %s %s %s %s\n",symname[vibsym[i]],symname[vibsym[i+1]],symname[vibsym[i+2]],
245 symname[vibsym[i+3]]);
246 fprintf(vibfile,"Int: %8.4f %8.4f %8.4f %8.4f\n",dipole[i],dipole[i+1],dipole[i+2],dipole[i+3]);
247 fprintf(vibfile,"\n");
248
249 for (j=0; j < maxvar; j++)
250 fprintf(vibfile,"%-d : %8.3f %8.3f %8.3f %8.3f\n",j,matrix[j][i],matrix[j][i+1],matrix[j][i+2],
251 matrix[j][i+3]);
252 }
253 }
254 fclose(vibfile); */
255
256 /* fprintf(pcmlogfile,"\n==========================================================\n");
257 fprintf(pcmlogfile,"Normal Mode Vibrational Analysis\n\n");
258 if (symmetry.ishape == 1)
259 fprintf(pcmlogfile,"Molecular Shape: LINEAR\n");
260 else if (symmetry.ishape == 2)
261 fprintf(pcmlogfile,"Molecular Shape: ASYMMETRIC ROTOR\n");
262 else if (symmetry.ishape == 3)
263 fprintf(pcmlogfile,"Molecular Shape: PROLATE SYMMETRIC ROTOR\n");
264 else if (symmetry.ishape == 4)
265 fprintf(pcmlogfile,"Molecular Shape: OBLATE SYMMETRIC ROTOR\n");
266 else if (symmetry.ishape == 5)
267 fprintf(pcmlogfile,"Molecular Shape: SPHERICAL ROTOR\n");
268 else if (symmetry.ishape == 6)
269 fprintf(pcmlogfile,"Molecular Shape: PLANAR\n");
270
271 fprintf(pcmlogfile,"Symmetry Point Group: %s\n",ngrp);
272 fprintf(pcmlogfile,"Symmetry Number: %d\n",symmetry.numsym);
273 if (symmetry.chiral)
274 fprintf(pcmlogfile,"Chirality: Yes\n");
275 else
276 fprintf(pcmlogfile,"Chirality: No\n"); */
277
278 if (symmetry.ishape == 1)
279 mm = maxvib - 5;
280 else
281 mm = maxvib - 6;
282
283 /* fprintf(pcmlogfile,"==========================================================\n");
284 fprintf(pcmlogfile,"\nFundmental Normal Vibrational Frequencies\n");
285 for (i=0; i < mm; i++)
286 fprintf(pcmlogfile," %-3d %8.3f %s %8.3f\n",i,eigen[i],symname[vibsym[i]],dipole[i]);
287
288 fprintf(pcmlogfile,"\nMoments of Inertia\n");
289 fprintf(pcmlogfile," IX = %8.3f IY = %8.3f IZ = %8.3f\n",symmetry.xi,symmetry.yi,symmetry.zi); */
290
291 strcpy(vibdata.ptgrp,ngrp);
292 vibdata.mom_ix = symmetry.xi;
293 vibdata.mom_iy = symmetry.yi;
294 vibdata.mom_iz = symmetry.zi;
295
296 if (symmetry.chiral)
297 imix = 2;
298 else
299 imix = 1;
300 // set values
301 temp = 298.15;
302 mass = 0.0;
303 for (i=1; i <= natom; i++)
304 mass += atomwt[i];
305 efact = 1.987;
306 cnosym = symmetry.numsym;
307 cimix = (double)imix;
308 etrans=htrans=strans=gtrans=cptrans=0.0;
309 erot=hrot=srot=grot=cprot=0.0;
310 evib=hvib=svib=gvib=cpvib=0.0;
311 emix=hmix=smix=gmix=cpmix=0.0;
312 etot=htot=stot=gtot=cptot=0.0;
313
314 // compute entropy, internal energy
315 strans = 2.2868*(5.0*log10(temp)+3.0*log10(mass))-2.3135;
316 etrans = 0.001987*1.5*temp;
317 if (strans < 0.0) strans = 0.0;
318 if (symmetry.ishape == 1)
319 {
320 xyzi = sqrt(symmetry.xi*symmetry.xi+symmetry.yi*symmetry.yi+symmetry.zi*symmetry.zi);
321 srot = 4.5736*(log10(temp)+log10(xyzi)-0.6049)+ efact;
322 erot = 0.001987*temp;
323 }else
324 {
325 xyzi = symmetry.xi*symmetry.yi*symmetry.zi;
326 srot = 2.2868*(3.0*log10(temp)+log10(xyzi)-1.3176)+ 2.9805;
327 erot = 0.001987*temp*1.5;
328 }
329 if (srot < 0.0) srot = 0.0;
330 ssym = -2.2868*2.0*log10(cnosym);
331 srot += ssym;
332
333 qptf = 0.0;
334 ezpe = 0.0;
335 entrov1 = 0.0;
336 entrov2 = 0.0;
337 for (i=0; i < mm; i++)
338 {
339 wdt = fabs(eigen[i])/temp;
340 efwt = exp(-1.43893*wdt);
341 qptf += efwt;
342 entrov1 += eigen[i]*efwt/(1.0-efwt);
343 entrov2 += log(1.0-efwt);
344
345 evib += eigen[i]*0.0028591*efwt/(1.0-efwt);
346 ezpe += 0.5*eigen[i]*0.0028591;
347 }
348 svib = efact*(1.438903*entrov1/temp - entrov2);
349 smix = 2.2868*2.0*log10(cimix);
350 stot = strans + srot + svib + smix;
351 evib += ezpe;
352 etot = etrans + erot + evib + emix + get_total_energy();
353
354 // compute enthalpy
355 htrans = etrans + 0.001987*temp;
356 hrot = erot;
357 hvib = evib;
358 hmix = emix;
359 htot = htrans + hrot + hvib + hmix + get_total_energy();
360 // compute free energy
361 gtrans = htrans - temp*strans*0.001;
362 grot = hrot - temp*srot*0.001;
363 gvib = hvib - temp*svib*0.001;
364 gmix = hmix - temp*smix*0.001;
365 gtot = gtrans + grot + gvib + gmix + get_total_energy();
366 // compute heat capacity
367 cptrans = 1.987 + 1.987*1.5;
368 if (symmetry.ishape == 1)
369 cprot = 1.987;
370 else
371 cprot = 1.987*1.5;
372 cpfact = 4.111/(temp*temp);
373 for (i=0; i < mm; i++)
374 {
375 wdt = fabs(eigen[i])/temp;
376 efwt = exp(-1.43893*wdt);
377 efwt2 = 1.0-efwt;
378 cpvib += eigen[i]*eigen[i]*cpfact*efwt/(efwt2*efwt2);
379 }
380 cptot = cptrans + cprot + cpvib + cpmix;
381 //
382 vibdata.etot = etot;
383 vibdata.htot = htot;
384 vibdata.stot = stot;
385 vibdata.gtot = gtot;
386 vibdata.cptot = cptot;
387 // print out
388 /* fprintf(pcmlogfile,"\nTemperatur: %8.3f K \n",temp);
389 fprintf(pcmlogfile,"\nZero Point Energy: %8.3f kcal/mol \n",ezpe);
390 fprintf(pcmlogfile,"\n Energy Enthalpy Entropy Free Energy Heat Capacity\n");
391 fprintf(pcmlogfile,"\n kcal/mol kcal/mol eu kcal/mol cal/mol/deg\n");
392 fprintf(pcmlogfile,"Translational: %10.3f %10.3f %10.3f %10.3f %10.3f\n",etrans,htrans,strans,gtrans,cptrans);
393 fprintf(pcmlogfile,"Rotational: %10.3f %10.3f %10.3f %10.3f %10.3f\n",erot,hrot,srot,grot,cprot);
394 fprintf(pcmlogfile,"Vibrational: %10.3f %10.3f %10.3f %10.3f %10.3f\n",evib,hvib,svib,gvib,cpvib);
395 fprintf(pcmlogfile,"Mixing: %10.3f %10.3f %10.3f %10.3f %10.3f\n",emix,hmix,smix,gmix,cpmix);
396 fprintf(pcmlogfile,"Total: %10.3f %10.3f %10.3f %10.3f %10.3f\n",etot,htot,stot,gtot,cptot); */
397
398 //
399 free_dvector(eigen, 0, 3*natom);
400 free_dmatrix(vects, 0, 3*natom, 0, 3*natom);
401 free_ivector(vibsym ,0,maxvar);
402 free_ivector(hinit ,0,maxvar);
403 free_ivector(hstop ,0,maxvar);
404 free_ivector(hindex ,0,maxhess);
405 free_dvector(hdiag ,0,maxvar);
406 free_dvector(h ,0,maxhess);
407
408 free_dmatrix(matrix ,0, maxvib, 0, maxvib);
409 free_dvector(a ,0, maxvib);
410 free_dvector(mass2, 0 , natom+1);
411 free_dvector(dipole,0,maxvar);
412 }
413 /* =============================================== */
414 void assign_vibsym(int natom, char *ngrp,char *vsym,int *nsym,double **a)
415 {
416 int iz;
417 double **b;
418
419 b = dmatrix (0, natom+1, 0,3);
420
421 *nsym = 0;
422 if (strcmp(ngrp,"Cs") == 0)
423 {
424 if (strcmp(vsym,"Cs") == 0)
425 *nsym = 1; // A'
426 else
427 *nsym = 2; // a''
428 } else if (strcmp(ngrp,"Ci") == 0)
429 {
430 if (strcmp(vsym,"Ci") == 0)
431 *nsym = 3; // ag
432 else
433 *nsym = 4; // au
434 } else if (strcmp(ngrp,"C1") == 0)
435 {
436 *nsym = 0; // a
437 } else if (strcmp(ngrp,"C2") == 0)
438 {
439 if (strcmp(vsym,"C2") == 0)
440 *nsym = 0; // a
441 else
442 *nsym = 15; // b
443 } else if (strcmp(ngrp,"C3") == 0)
444 {
445 if (strcmp(vsym,"C3") == 0)
446 *nsym = 0; // a
447 else
448 *nsym = 27; // e
449 } else if (strcmp(ngrp,"C4") == 0)
450 {
451 if (strcmp(vsym,"C4") == 0)
452 *nsym = 0; // a
453 else
454 {
455 if (strcmp(vsym,"C2") == 0)
456 *nsym = 15; // b
457 else
458 *nsym = 27; // e
459 }
460 } else if (strcmp(ngrp,"C5") == 0)
461 {
462 if (strcmp(vsym,"C5") == 0)
463 *nsym = 0; // a
464 else
465 *nsym = 28; // e?
466 } else if (strcmp(ngrp,"C6") == 0)
467 {
468 if (strcmp(vsym,"C6") == 0)
469 *nsym = 0; // a
470 else
471 {
472 if (strcmp(vsym,"C3") == 0)
473 *nsym = 15; // b
474 else
475 *nsym = 27; // e
476 }
477 } else if (strcmp(ngrp,"C2v") == 0)
478 {
479 if (strcmp(vsym,"C2v") == 0)
480 *nsym = 5; // a1
481 else
482 {
483 if (strcmp(vsym,"C2") == 0)
484 *nsym = 6; //a2
485 else
486 {
487 iz = findh(a,b,3);
488 if (iz)
489 *nsym = 16; // b1
490 else
491 *nsym = 17; // b2
492 }
493 }
494 } else if (strcmp(ngrp,"C3v") == 0)
495 {
496 if (strcmp(vsym,"C3v") == 0)
497 *nsym = 5; // a1
498 else
499 {
500 if (strcmp(vsym,"C3") == 0)
501 *nsym = 6; // a2
502 else
503 *nsym = 27; // e
504 }
505 } else if (strcmp(ngrp,"C4v") == 0)
506 {
507 if (strcmp(vsym,"C4v") == 0)
508 *nsym = 5; // a1
509 else
510 {
511 if (strcmp(vsym,"C4") == 0)
512 *nsym = 6; // a2
513 else
514 {
515 if (strcmp(vsym,"C2v") == 0)
516 *nsym = 16; //b?
517 else
518 *nsym = 27; // e
519 }
520 }
521 } else if (strcmp(ngrp,"C5v") == 0)
522 {
523 if (strcmp(vsym,"C5v") == 0)
524 *nsym = 5; // a1
525 else
526 {
527 if (strcmp(vsym,"C5") == 0)
528 *nsym = 6; // a2
529 else
530 *nsym = 28; // e?
531 }
532 } else if (strcmp(ngrp,"C6v") == 0)
533 {
534 if (strcmp(vsym,"C6v") == 0)
535 *nsym = 5; // a1
536 else
537 {
538 if (strcmp(vsym,"C6") == 0)
539 *nsym = 6; // a2
540 else
541 {
542 if (strcmp(vsym,"C3v") == 0)
543 *nsym = 16; // b?
544 else
545 *nsym = 28; // e?
546 }
547 }
548 } else if (strcmp(ngrp,"C2h") == 0)
549 {
550 if (strcmp(vsym,"C2h") == 0)
551 *nsym = 3; // ag
552 else
553 {
554 if (strcmp(vsym,"C2") == 0)
555 *nsym = 4; // au
556 else
557 {
558 if (strcmp(vsym,"Ci") == 0)
559 *nsym = 19; // bg
560 else
561 *nsym = 20; // bu
562 }
563 }
564 } else if (strcmp(ngrp,"C3h") == 0)
565 {
566 if (strcmp(vsym,"C3h") == 0)
567 *nsym = 1; // a'
568 else
569 {
570 if (strcmp(vsym,"C3") == 0)
571 *nsym = 2; // a"
572 else
573 {
574 if (strcmp(vsym,"Cs") == 0)
575 *nsym = 31; // e'
576 else
577 *nsym = 32; // e"
578 }
579 }
580 } else if (strcmp(ngrp,"C4h") == 0)
581 {
582 if (strcmp(vsym,"C4h") == 0)
583 *nsym = 3; // ag
584 else
585 {
586 if (strcmp(vsym,"C2h") == 0)
587 *nsym = 19; // Bg
588 else
589 {
590 if (strcmp(vsym,"C4") == 0)
591 *nsym = 4; // au
592 else
593 {
594 if (strcmp(vsym,"S4") == 0)
595 *nsym = 20; // bu
596 else
597 {
598 iz = findi(a,b);
599 if (iz)
600 *nsym = 33 ; // eg
601 else
602 *nsym = 34 ; // eu
603 }
604 }
605 }
606 }
607 } else if (strcmp(ngrp,"C5h") == 0)
608 {
609 if (strcmp(vsym,"C5h") == 0)
610 *nsym = 1; // a'
611 else
612 {
613 if (strcmp(vsym,"C5") == 0)
614 *nsym = 2; // a"
615 else
616 {
617 if (strcmp(vsym,"Cs") == 0)
618 *nsym = 35; // e?'
619 else
620 *nsym = 36; // e?"
621 }
622 }
623 } else if (strcmp(ngrp,"C6h") == 0)
624 {
625 if (strcmp(vsym,"C6h") == 0)
626 *nsym = 3; // ag
627 else
628 {
629 if (strcmp(vsym,"C6") == 0)
630 *nsym = 4; // au
631 else
632 {
633 if (strcmp(vsym,"S6") == 0)
634 *nsym = 19; // bg
635 else
636 {
637 if (strcmp(vsym,"C3h") == 0)
638 *nsym = 20; // bu
639 else
640 {
641 if (strcmp(vsym,"C2") == 0)
642 {
643 iz = findi(a,b);
644 if (iz)
645 *nsym = 41 ; // e2g
646 else
647 *nsym = 42 ; // e2u
648 } else
649 {
650 iz = findi(a,b);
651 if (iz)
652 *nsym = 39; // e1g
653 else
654 *nsym = 40; // e1u
655 }
656 }
657 }
658 }
659 }
660 } else if (strcmp(ngrp,"D2") == 0)
661 {
662 if (strcmp(vsym,"D2") == 0)
663 *nsym = 0; // a
664 else
665 {
666 if (strcmp(vsym,"C2") == 0)
667 *nsym = 16; // b1 or b2
668 }
669 } else if (strcmp(ngrp,"D3") == 0)
670 {
671 if (strcmp(vsym,"D3") == 0)
672 *nsym = 5; // a1
673 else
674 {
675 if (strcmp(vsym,"C3") == 0)
676 *nsym = 6; // a2
677 else
678 *nsym = 27; // e
679 }
680 } else if (strcmp(ngrp,"D4") == 0)
681 {
682 if (strcmp(vsym,"D4") == 0)
683 *nsym = 5; // a1
684 else
685 {
686 if (strcmp(vsym,"C4") == 0)
687 *nsym = 6; // a2
688 else
689 {
690 if (strcmp(vsym,"C1") == 0)
691 *nsym = 27; // e
692 else
693 {
694 if (strcmp(vsym,"D2") == 0)
695 *nsym = 16; // b1
696 else
697 *nsym = 17; // b2
698 }
699 }
700 }
701 } else if (strcmp(ngrp,"D5") == 0)
702 {
703 if (strcmp(vsym,"D5") == 0)
704 *nsym = 5; // a1
705 else
706 {
707 if (strcmp(vsym,"C5") == 0)
708 *nsym = 6; // a2
709 else
710 *nsym = 28; // e?
711 }
712 } else if (strcmp(ngrp,"D6") == 0)
713 {
714 if (strcmp(vsym,"D6") == 0)
715 *nsym = 5; // a1
716 else
717 {
718 if (strcmp(vsym,"C6") == 0)
719 *nsym = 6; // a2
720 else
721 {
722 if (strcmp(vsym,"D3") == 0)
723 *nsym = 16; // b1
724 else
725 {
726 if (strcmp(vsym,"C3") == 0)
727 *nsym = 17; // b2
728 else
729 {
730 if (strcmp(vsym,"C2") == 0)
731 *nsym = 29; // e2
732 else
733 *nsym = 28; // e1
734 }
735 }
736 }
737 }
738 } else if (strcmp(ngrp,"D2h") == 0)
739 {
740 if (strcmp(vsym,"D2h") == 0)
741 *nsym = 3; // ag
742 else
743 {
744 if (strcmp(vsym,"D2") == 0)
745 *nsym = 4; // au
746 else
747 {
748 iz = findi(a,b);
749 if (iz)
750 *nsym = 21; // b1g
751 else
752 *nsym = 22; // b1u?
753 }
754 }
755 } else if (strcmp(ngrp,"D3h") == 0)
756 {
757 if (strcmp(vsym,"D3h") == 0)
758 *nsym = 5; // a1
759 else
760 {
761 if (strcmp(vsym,"C3h") == 0)
762 *nsym = 6; // a2
763 else
764 {
765 if (strcmp(vsym,"Cs") == 0)
766 *nsym = 31; // e'
767 else
768 {
769 if (strcmp(vsym,"D3") == 0)
770 *nsym = 5; // a1
771 else
772 {
773 if (strcmp(vsym,"C3v") == 0)
774 *nsym = 10; // a2"
775 else
776 *nsym = 32; // e"
777 }
778 }
779 }
780 }
781 } else if (strcmp(ngrp,"D4h") == 0)
782 {
783 if (strcmp(vsym,"D4h") == 0)
784 *nsym = 11; // a1g
785 else
786 {
787 if (strcmp(vsym,"C4h") == 0)
788 *nsym = 13; // a2g
789 else
790 {
791 if (strcmp(vsym,"D4") == 0)
792 *nsym = 12; // a1u
793 else
794 {
795 if (strcmp(vsym,"C4v") == 0)
796 *nsym = 14; // a2u
797 else
798 {
799 iz = findi(a,b);
800 if (iz)
801 {
802 if (strcmp(vsym,"D2h") == 0)
803 *nsym = 21; // b1g
804 else
805 *nsym = 33; // eg
806 } else
807 {
808 if (strcmp(vsym,"D2d") == 0)
809 *nsym = 22; // b?u
810 else
811 *nsym = 34; //eu
812 }
813 }
814 }
815 }
816 }
817 } else if (strcmp(ngrp,"D5h") == 0)
818 {
819 if (strcmp(vsym,"D5h") == 0)
820 *nsym = 7; // a1'
821 else
822 {
823 if (strcmp(vsym,"C5h") == 0)
824 *nsym = 8; // a2'
825 else
826 {
827 if (strcmp(vsym,"D5") == 0)
828 *nsym = 9; // a1''
829 else
830 {
831 if (strcmp(vsym,"C5v") == 0)
832 *nsym = 10; // a2''
833 else
834 *nsym = 36; // e"
835 }
836 }
837 }
838 } else if (strcmp(ngrp,"D6h") == 0)
839 {
840 if (strcmp(vsym,"D6h") == 0)
841 *nsym = 11; // a1g
842 else
843 {
844 if (strcmp(vsym,"C6h") == 0)
845 *nsym = 13; // a2g
846 else
847 {
848 if (strcmp(vsym,"D3d") == 0)
849 *nsym = 21; // b?g
850 else
851 {
852 if (strcmp(vsym,"D6") == 0)
853 *nsym = 12; // a1u
854 else
855 {
856 if (strcmp(vsym,"C6v") == 0)
857 *nsym = 14; // a2u
858 else
859 {
860 if (strcmp(vsym,"D3h") == 0)
861 *nsym = 22; // b?u
862 else
863 {
864 if ( strcmp(vsym,"C2") == 0 || strcmp(vsym,"D2") == 0 ||
865 strcmp(vsym,"C2h") == 0 )
866 {
867 iz = findi(a,b);
868 if (iz)
869 *nsym = 41; // e2g
870 else
871 *nsym = 42; // e2u
872 } else
873 {
874 iz = findi(a,b);
875 if (iz)
876 *nsym = 39; // e1g
877 else
878 *nsym = 40; // e1u
879 }
880 }
881 }
882 }
883 }
884 }
885 }
886 } else if (strcmp(ngrp,"D2d") == 0)
887 {
888 if (strcmp(vsym,"D2d") == 0 )
889 *nsym = 5; // a1
890 else
891 {
892 if (strcmp(vsym,"S4") == 0)
893 *nsym = 6; // a2
894 else
895 {
896 if (strcmp(vsym,"D2") == 0)
897 *nsym = 16; // b1
898 else
899 {
900 if (strcmp(vsym,"C2v") == 0)
901 *nsym = 17; // b2
902 else
903 *nsym = 27; // e
904 }
905 }
906 }
907 } else if (strcmp(ngrp,"D3d") == 0)
908 {
909 iz = findi(a,b);
910 if (iz)
911 {
912 if (strcmp(vsym,"D3d") == 0 )
913 *nsym = 11; // a1g
914 else
915 {
916 if (strcmp(vsym,"S6") == 0 )
917 *nsym = 13; // a2g
918 else
919 *nsym = 33; // eg
920 }
921 } else
922 {
923 if (strcmp(vsym,"D3") == 0 )
924 *nsym = 12; // a1u
925 else
926 {
927 if (strcmp(vsym,"C3v") == 0 )
928 *nsym = 14; // a2u
929 else
930 *nsym = 34; // eu
931 }
932 }
933 } else if (strcmp(ngrp,"D4d") == 0)
934 {
935 if (strcmp(vsym,"D4d") == 0 )
936 *nsym = 5; // a1
937 else
938 {
939 if (strcmp(vsym,"S8") == 0)
940 *nsym = 6; // a2
941 else
942 {
943 if (strcmp(vsym,"D4") == 0)
944 *nsym = 16; // b1
945 else
946 {
947 if (strcmp(vsym,"C4v") == 0)
948 *nsym = 17; // b2
949 else
950 *nsym = 27; // e
951 }
952 }
953 }
954 } else if (strcmp(ngrp,"D5d") == 0)
955 {
956 iz = findi(a,b);
957 if (iz)
958 {
959 if (strcmp(vsym,"D5d") == 0 )
960 *nsym = 11; // a1g
961 else
962 {
963 if (strcmp(vsym,"S10") == 0 )
964 *nsym = 13; // a2g
965 else
966 *nsym = 33; // eg
967 }
968 } else
969 {
970 if (strcmp(vsym,"D5") == 0 )
971 *nsym = 12; // a1u
972 else
973 {
974 if (strcmp(vsym,"C5v") == 0 )
975 *nsym = 14; // a2u
976 else
977 *nsym = 34; // eu
978 }
979 }
980 } else if (strcmp(ngrp,"D6d") == 0)
981 {
982 if (strcmp(vsym,"D6d") == 0 )
983 *nsym = 5; // a1
984 else
985 {
986 if (strcmp(vsym,"S12") == 0)
987 *nsym = 6; // a2
988 else
989 {
990 if (strcmp(vsym,"D6") == 0)
991 *nsym = 16; // b1
992 else
993 {
994 if (strcmp(vsym,"C6v") == 0)
995 *nsym = 17; // b2
996 else
997 *nsym = 27; // e
998 }
999 }
1000 }
1001 } else if (strcmp(ngrp,"S4") == 0)
1002 {
1003 if (strcmp(vsym,"S4") == 0 )
1004 *nsym = 0; // a
1005 else
1006 {
1007 if (strcmp(vsym,"C2") == 0)
1008 *nsym = 15; // b
1009 else
1010 *nsym = 27; // e
1011 }
1012 } else if (strcmp(ngrp,"S6") == 0)
1013 {
1014 iz = findi(a,b);
1015 if (iz)
1016 {
1017 if (strcmp(vsym,"S6") == 0 )
1018 *nsym = 3; // ag
1019 else
1020 *nsym = 33; // eg
1021 } else
1022 {
1023 if (strcmp(vsym,"C3") == 0 )
1024 *nsym = 4; // au
1025 else
1026 *nsym = 34; // eu
1027 }
1028 } else if (strcmp(ngrp,"S8") == 0)
1029 {
1030 if (strcmp(vsym,"S8") == 0 )
1031 *nsym = 0; // a
1032 else
1033 {
1034 if (strcmp(vsym,"C4") == 0)
1035 *nsym = 15; // b
1036 else
1037 *nsym = 27; // e
1038 }
1039 } else if (strcmp(ngrp,"T") == 0)
1040 {
1041 if (strcmp(vsym,"T") == 0 || strcmp(vsym,"C3") == 0)
1042 *nsym = 0; // a
1043 else
1044 {
1045 if (strcmp(vsym,"D2") == 0 || strcmp(vsym,"C2") == 0)
1046 *nsym = 27; // e
1047 else
1048 *nsym = 43; // t
1049 }
1050 } else if (strcmp(ngrp,"Th") == 0)
1051 {
1052 if (strcmp(vsym,"Th") == 0)
1053 *nsym = 3; // ag
1054 else
1055 {
1056 if (strcmp(vsym,"T") == 0)
1057 *nsym = 4; // au
1058 else
1059 {
1060 iz = findi(a,b);
1061 if (iz)
1062 {
1063 if (strcmp(vsym,"D2h") == 0 || strcmp(vsym,"D2") == 0)
1064 *nsym = 33; // Eg
1065 else
1066 *nsym = 44; // Tg
1067 } else
1068 {
1069 if (strcmp(vsym,"D2") == 0 || strcmp(vsym,"C2") == 0)
1070 *nsym = 34; // eu
1071 else
1072 *nsym = 45; // Tu
1073 }
1074 }
1075 }
1076 } else if (strcmp(ngrp,"Td") == 0)
1077 {
1078 if (strcmp(vsym,"Td") == 0)
1079 *nsym = 5; // a1
1080 else
1081 {
1082 if (strcmp(vsym,"T") == 0)
1083 *nsym = 6; // a2
1084 else
1085 {
1086 if (strcmp(vsym,"D2") == 0 || strcmp(vsym,"D2d") == 0)
1087 *nsym = 27; // e
1088 else
1089 *nsym = 43; // T
1090 }
1091 }
1092 } else if (strcmp(ngrp,"O") == 0)
1093 {
1094 if (strcmp(vsym,"O") == 0)
1095 *nsym = 5; // a1
1096 else
1097 {
1098 if (strcmp(vsym,"T") == 0)
1099 *nsym = 6; // a2
1100 else
1101 {
1102 if (strcmp(vsym,"D2") == 0)
1103 *nsym = 27; // e
1104 else
1105 *nsym = 43; // T
1106 }
1107 }
1108 }
1109 free_dmatrix(b ,0, natom+1, 0,3);
1110 }
1111 // ======================
1112 void tred2(double **a,int n, double *d, double *e)
1113 {
1114
1115 int i,j,k,l, nn;
1116 double scale, hh, h, g, f;
1117
1118 nn = n-1;
1119 for (i= nn; i >= 1; i--)
1120 {
1121 l = i - 1;
1122 h = scale = 0.0;
1123 if (l > 0)
1124 {
1125 for (k=0; k <= l; k++)
1126 scale += fabs(a[i][k]);
1127 if (scale == 0.0)
1128 e[i] = a[i][l];
1129 else
1130 {
1131 for (k=0; k <= l; k++)
1132 {
1133 a[i][k] /= scale;
1134 h += a[i][k]*a[i][k];
1135 }
1136 f = a[i][l];
1137 g = f>0 ? -sqrt(h) : sqrt(h);
1138 e[i] = scale*g;
1139 h -= f*g;
1140 a[i][l] = f-g;
1141 f = 0.0;
1142 for (j=0; j <= l; j++)
1143 {
1144 a[j][i] = a[i][j]/h;
1145 g = 0.0;
1146 for (k=0; k <= j; k++)
1147 g += a[j][k]*a[i][k];
1148 for (k = j+1; k <= l; k++)
1149 g += a[k][j]*a[i][k];
1150 e[j] = g/h;
1151 f += e[j]*a[i][j];
1152 }
1153 hh = f/(h+h);
1154 for (j=0; j <= l; j++)
1155 {
1156 f = a[i][j];
1157 e[j] = g = e[j] - hh*f;
1158 for (k=0; k <= j; k++)
1159 a[j][k] -= (f*e[k]+g*a[i][k]);
1160 }
1161 }
1162 }else
1163 e[i] = a[i][l];
1164 d[i] = h;
1165 }
1166 d[0] = 0.0;
1167 e[0] = 0.0;
1168 for (i=0; i < n; i++)
1169 {
1170 l = i - 1;
1171 if (i > 0)
1172 {
1173 for (j=0; j <= l; j++)
1174 {
1175 g = 0.0;
1176 for (k=0; k <= l; k++)
1177 g += a[i][k]*a[k][j];
1178 for (k=0; k <= l; k++)
1179 a[k][j] -= g*a[k][i];
1180 }
1181 }
1182 d[i] = a[i][i];
1183 a[i][i] = 1.0;
1184 for (j=0; j <= l; j++)
1185 {
1186 a[j][i] = a[i][j] = 0.0;
1187 }
1188 }
1189 }
1190 /* =========================================== */
1191 void tqli(double *d, double *e, int n, double **z)
1192 {
1193 int m,l,iter,i,k,j;
1194 double s,r,p,g,f,dd,c,b;
1195
1196 for (i=1; i < n; i++)
1197 e[i-1] = e[i];
1198 e[n-1] = 0.0;
1199
1200 for (l=0; l < n; l++)
1201 {
1202 iter = 0;
1203 do
1204 {
1205 for (m=l; m < n-1; m++)
1206 {
1207 dd = fabs(d[m]) + fabs(d[m+1]);
1208 if (fabs(e[m])+dd == dd) break;
1209 }
1210 if (m!= l)
1211 {
1212 if (iter++ == 30)
1213 {
1214 fprintf(pcmlogfile,"Error in tqli\n");
1215 message_alert("Error in Tqli","Error");
1216 return;
1217 }
1218 g = (d[l+1]-d[l])/(2.0*e[l]);
1219 r = sqrt((g*g)+1.0);
1220 g = d[m]-d[l]+e[l]/(g+SIGN(r,g));
1221 s = c = 1.0;
1222 p = 0.0;
1223 for (i=m-1; i >= l; i--)
1224 {
1225 f = s*e[i];
1226 b = c*e[i];
1227 if (fabs(f) >= fabs(g))
1228 {
1229 c = g/f;
1230 r = sqrt((c*c)+1.0);
1231 e[i+1] = f*r;
1232 c *= (s=1.0/r);
1233 } else
1234 {
1235 s = f/g;
1236 r = sqrt((s*s)+1.0);
1237 e[i+1] = g*r;
1238 s *= (c=1.0/r);
1239 }
1240 g = d[i+1]-p;
1241 r = (d[i]-g)*s+2.0*c*b;
1242 p = s*r;
1243 d[i+1] = g+p;
1244 g = c*r-b;
1245 for (k=0; k < n; k++)
1246 {
1247 f = z[k][i+1];
1248 z[k][i+1] = s*z[k][i]+c*f;
1249 z[k][i] = c*z[k][i]-s*f;
1250 }
1251 }
1252 d[l] = d[l]-p;
1253 e[l] = g;
1254 e[m] = 0.0;
1255 }
1256 } while (m!= l);
1257 }
1258 // sort d and z
1259 for (i=0; i < n; i++)
1260 {
1261 p = d[k=i];
1262 for (j=i+1; j < n; j++)
1263 if (d[j] >= p) p = d[k=j];
1264 if (k != i)
1265 {
1266 d[k] = d[i];
1267 d[i] = p;
1268 for (j=0; j < n; j++)
1269 {
1270 p = z[j][i];
1271 z[j][i] = z[j][k];
1272 z[j][k] = p;
1273 }
1274 }
1275 }
1276 }