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 (11 years, 4 months ago) by gilbertke
File size: 37467 byte(s)
Log Message:
further cleanup and localization of atom data
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4 gilbertke 63
5     #include "asnsym.h"
6 tjod 3 #include "utility.h"
7     #include "bonds_ff.h"
8 gilbertke 63 #include "vibrate.h"
9 tjod 3
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 gilbertke 63
20 tjod 3 EXTERN double hesscut;
21    
22 gilbertke 63 // FILE * fopen_path ( char * , char * , char * ) ;
23 tjod 3
24 gilbertke 63 static void tred2(double **,int , double *, double *);
25     static void tqli(double *,double *,int ,double **);
26 gilbertke 110 static void assign_vibsym(int, char *,char *,int *,double **);
27 gilbertke 103 double get_total_energy(void);
28     void hessian(int, double *,int *, int *, int *,double *);
29 gilbertke 110 void vibrate(int natom,double *atomwt,double *x,double *y,double *z,double *charge);
30 gilbertke 63 //void vibcharge_dipole(double **);
31     //void vibdipole_dipole(double **);
32    
33 gilbertke 110 void vibrate(int natom,double *atomwt,double *x,double *y,double *z,double *charge)
34 tjod 3 {
35 gilbertke 103 int i, k, j, ii, jj, nsym,imix,mm;
36 tjod 3 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 gilbertke 103 double *dipole;
52 tjod 3 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 gilbertke 110 if (atomwt[i] <= 0.0)
62     atomwt[i] = 0.001;
63 tjod 3 }
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 gilbertke 110 mass2[i] = sqrt(atomwt[i]);
85 tjod 3
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 gilbertke 110 matrix[i][i] = hdiag[i]/atomwt[jj];
111 tjod 3 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 gilbertke 110 acoord[i][0] = x[i];
161     acoord[i][1] = y[i];
162     acoord[i][2] = z[i];
163 tjod 3 }
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 gilbertke 110 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 tjod 3
189     }
190    
191     ptgrp(1,vsym,acoord);
192 gilbertke 110 assign_vibsym(natom,ngrp,vsym,&nsym,acoord);
193 tjod 3 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 wdelano 58 /* fprintf(pcmlogfile,"\n==========================================================\n");
257     fprintf(pcmlogfile,"Normal Mode Vibrational Analysis\n\n");
258 tjod 3 if (symmetry.ishape == 1)
259 wdelano 58 fprintf(pcmlogfile,"Molecular Shape: LINEAR\n");
260 tjod 3 else if (symmetry.ishape == 2)
261 wdelano 58 fprintf(pcmlogfile,"Molecular Shape: ASYMMETRIC ROTOR\n");
262 tjod 3 else if (symmetry.ishape == 3)
263 wdelano 58 fprintf(pcmlogfile,"Molecular Shape: PROLATE SYMMETRIC ROTOR\n");
264 tjod 3 else if (symmetry.ishape == 4)
265 wdelano 58 fprintf(pcmlogfile,"Molecular Shape: OBLATE SYMMETRIC ROTOR\n");
266 tjod 3 else if (symmetry.ishape == 5)
267 wdelano 58 fprintf(pcmlogfile,"Molecular Shape: SPHERICAL ROTOR\n");
268 tjod 3 else if (symmetry.ishape == 6)
269 wdelano 58 fprintf(pcmlogfile,"Molecular Shape: PLANAR\n");
270 tjod 3
271 wdelano 58 fprintf(pcmlogfile,"Symmetry Point Group: %s\n",ngrp);
272     fprintf(pcmlogfile,"Symmetry Number: %d\n",symmetry.numsym);
273 tjod 3 if (symmetry.chiral)
274 wdelano 58 fprintf(pcmlogfile,"Chirality: Yes\n");
275 tjod 3 else
276 wdelano 58 fprintf(pcmlogfile,"Chirality: No\n"); */
277 tjod 3
278     if (symmetry.ishape == 1)
279     mm = maxvib - 5;
280     else
281     mm = maxvib - 6;
282    
283 wdelano 58 /* fprintf(pcmlogfile,"==========================================================\n");
284     fprintf(pcmlogfile,"\nFundmental Normal Vibrational Frequencies\n");
285 tjod 3 for (i=0; i < mm; i++)
286 wdelano 58 fprintf(pcmlogfile," %-3d %8.3f %s %8.3f\n",i,eigen[i],symname[vibsym[i]],dipole[i]);
287 tjod 3
288 wdelano 58 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 tjod 3
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 gilbertke 110 mass += atomwt[i];
305 tjod 3 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 gilbertke 103 etot = etrans + erot + evib + emix + get_total_energy();
353 tjod 3
354     // compute enthalpy
355     htrans = etrans + 0.001987*temp;
356     hrot = erot;
357     hvib = evib;
358     hmix = emix;
359 gilbertke 103 htot = htrans + hrot + hvib + hmix + get_total_energy();
360 tjod 3 // 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 gilbertke 103 gtot = gtrans + grot + gvib + gmix + get_total_energy();
366 tjod 3 // 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 wdelano 58 /* 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 tjod 3
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 gilbertke 110 void assign_vibsym(int natom, char *ngrp,char *vsym,int *nsym,double **a)
415 tjod 3 {
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 gilbertke 63 // ======================
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     }