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, 4 months ago) by wdelano
File size: 37460 byte(s)
Log Message: 
branch created
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "asnsym.h"
7 #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 #include "vibrate.h"
14
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
25 EXTERN double hesscut;
26
27 // FILE * fopen_path ( char * , char * , char * ) ;
28
29 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 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 /* fprintf(pcmlogfile,"\n==========================================================\n");
263 fprintf(pcmlogfile,"Normal Mode Vibrational Analysis\n\n");
264 if (symmetry.ishape == 1)
265 fprintf(pcmlogfile,"Molecular Shape: LINEAR\n");
266 else if (symmetry.ishape == 2)
267 fprintf(pcmlogfile,"Molecular Shape: ASYMMETRIC ROTOR\n");
268 else if (symmetry.ishape == 3)
269 fprintf(pcmlogfile,"Molecular Shape: PROLATE SYMMETRIC ROTOR\n");
270 else if (symmetry.ishape == 4)
271 fprintf(pcmlogfile,"Molecular Shape: OBLATE SYMMETRIC ROTOR\n");
272 else if (symmetry.ishape == 5)
273 fprintf(pcmlogfile,"Molecular Shape: SPHERICAL ROTOR\n");
274 else if (symmetry.ishape == 6)
275 fprintf(pcmlogfile,"Molecular Shape: PLANAR\n");
276
277 fprintf(pcmlogfile,"Symmetry Point Group: %s\n",ngrp);
278 fprintf(pcmlogfile,"Symmetry Number: %d\n",symmetry.numsym);
279 if (symmetry.chiral)
280 fprintf(pcmlogfile,"Chirality: Yes\n");
281 else
282 fprintf(pcmlogfile,"Chirality: No\n"); */
283
284 if (symmetry.ishape == 1)
285 mm = maxvib - 5;
286 else
287 mm = maxvib - 6;
288
289 /* fprintf(pcmlogfile,"==========================================================\n");
290 fprintf(pcmlogfile,"\nFundmental Normal Vibrational Frequencies\n");
291 for (i=0; i < mm; i++)
292 fprintf(pcmlogfile," %-3d %8.3f %s %8.3f\n",i,eigen[i],symname[vibsym[i]],dipole[i]);
293
294 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
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 /* 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
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 // ======================
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 }