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