ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/vibrate.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (12 years, 4 months ago) by tjod
File size: 33750 byte(s)
Log Message:
test

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