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 (11 years, 4 months ago) by tjod
File size: 33750 byte(s)
Log Message:
test

Line File contents
1 #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 }