1 /*
2 * (c) 2004 Iowa State University
3 * see the LICENSE file in the top level directory
4 */
5
6 /* BasisSet.cpp
7
8 Basis Set classes
9
10 Seperated from other files 4/1998 BMB
11 */
12
13 #include "Globals.h"
14 #include "GlobalExceptions.h"
15 #include "BFiles.h"
16 #include "MyTypes.h"
17 #include "BasisSet.h"
18 #include <math.h>
19 #include <string.h>
20 #include <cctype>
21
22 // Basis set related routines
23
24 // Constructor and destructor
BasisSet(long nAtoms,long nShells)25 BasisSet::BasisSet(long nAtoms, long nShells) {
26 Shells.reserve(nShells);
27 NumShells = NumFuncs = 0;
28 MapLength = nAtoms;
29 goodCharges = false;
30 BasisMap.reserve(2*nAtoms);
31 NuclearCharge.reserve(nAtoms);
32 for (long i=0; i<nAtoms; i++) {
33 BasisMap.push_back(0);
34 BasisMap.push_back(-1); //clear the basismap
35 NuclearCharge.push_back(-1);
36 }
37 }
~BasisSet(void)38 BasisSet::~BasisSet(void) {
39 }
40
GetNumBasisFuncs(bool UseSphericalHarmonics) const41 long BasisSet::GetNumBasisFuncs(bool UseSphericalHarmonics) const {
42 long result;
43 if (!UseSphericalHarmonics)
44 result = NumFuncs;
45 else {
46 result=0;
47 for (long atom=0; atom<MapLength; atom++) {
48 for (long i=BasisMap[2*atom]; i<=BasisMap[2*atom+1]; i++)
49 result += Shells[i].GetNumFuncs(UseSphericalHarmonics);
50 }
51 }
52 return result;
53 }
GetShellIndexArray(long * IndexArray) const54 void BasisSet::GetShellIndexArray(long * IndexArray) const {
55 long i, ShellStart, ShellEnd, j, index=0, shell=0;
56 for (i=0; i<NumFuncs; i++) IndexArray[i] = -1;
57 for (i=0; i<MapLength; i++) {
58 ShellStart = BasisMap[2*i];
59 ShellEnd = BasisMap[2*i+1];
60 for (j=ShellStart; j<=ShellEnd; j++) {
61 IndexArray[shell] = index;
62 index += Shells[j].GetNumFuncs(false);
63 shell++;
64 }
65 }
66 }
Normalize(bool,bool)67 void BasisSet::Normalize(bool /*InputNormed*/, bool /*NormOutput*/) {
68 //Set up the normalization constants: (2/pi)^3/4 * 2^n/n!!
69 double SNorm = sqrt(sqrt(pow(2.0/kPi, 3)));
70 double PNorm = SNorm*2;
71 double DNorm = PNorm*2/sqrt(3.0);
72 double FNorm = DNorm*2/sqrt(5.0);
73 double GNorm = FNorm*2/sqrt(7.0);
74 double HNorm = GNorm*2/sqrt(9.0);
75 double INorm = HNorm*2/sqrt(11.0);
76 for (long ishell=0; ishell<NumShells; ishell++) {
77 float fExp;
78 long NumPrims;
79 short ShellType = Shells[ishell].ShellType;
80 std::vector<float> & NormCoef = Shells[ishell].NormCoef;
81 std::vector<float> & InputCoef = Shells[ishell].InputCoef;
82 NumPrims = Shells[ishell].NumPrims;
83 for (long iprim=0; iprim<NumPrims; iprim++) {
84 fExp = Shells[ishell].Exponent[iprim];
85 double ExpNorm = sqrt(sqrt(fExp));
86 double ExpNorm2 = ExpNorm*ExpNorm;
87 ExpNorm *= ExpNorm2;
88 //Normalization constant is: [fExp^(n+3/2) 2^(2n+3/2)/(pi^3/2 n!!)]^1/2
89 //This factor is correct for functions of the form xx, xxx, etc but is off
90 //for the xy, xxy, etc functions which must be multiplied by sqrt3, sqrt5 etc later
91 switch (ShellType) {
92 case LShell:
93 NormCoef[iprim+NumPrims] = (float) (PNorm*ExpNorm*ExpNorm2*InputCoef[iprim+NumPrims]);
94 case SShell:
95 NormCoef[iprim] = (float) (SNorm*ExpNorm*InputCoef[iprim]);
96 break;
97 case PShell:
98 NormCoef[iprim] = (float) (PNorm*ExpNorm*ExpNorm2*InputCoef[iprim]);
99 break;
100 case DShell:
101 NormCoef[iprim] = (float) (DNorm*fExp*ExpNorm*InputCoef[iprim]);
102 break;
103 case FShell:
104 NormCoef[iprim] = (float) (FNorm*fExp*ExpNorm*ExpNorm2*InputCoef[iprim]);
105 break;
106 case GShell:
107 NormCoef[iprim] = (float) (GNorm*fExp*fExp*ExpNorm*InputCoef[iprim]);
108 break;
109 //These two need to be checked!!!!
110 case HShell:
111 NormCoef[iprim] = (float) (HNorm*fExp*fExp*ExpNorm*ExpNorm2*InputCoef[iprim]);
112 break;
113 case IShell:
114 NormCoef[iprim] = (float) (INorm*fExp*fExp*fExp*ExpNorm*InputCoef[iprim]);
115 break;
116 }
117 }
118 }
119 }
WriteBasis(BufferFile * File,long AtomNum) const120 void BasisSet::WriteBasis(BufferFile * File, long AtomNum) const {
121 double PI = acos(-1.0);
122 double PI32 = PI * sqrt(PI);
123 if ((AtomNum < MapLength)&&(AtomNum>=0)) {
124 long StartShell = BasisMap[2*AtomNum];
125 long EndShell = BasisMap[2*AtomNum + 1];
126 char ShellLabel, Line[kMaxLineLength];
127 long ishell, iprim;
128 for (ishell=StartShell; ishell<=EndShell; ishell++) {
129 switch (Shells[ishell].ShellType) {
130 case LShell:
131 ShellLabel = 'L';
132 break;
133 case SShell:
134 ShellLabel = 'S';
135 break;
136 case PShell:
137 ShellLabel = 'P';
138 break;
139 case DShell:
140 ShellLabel = 'D';
141 break;
142 case FShell:
143 ShellLabel = 'F';
144 break;
145 case GShell:
146 ShellLabel = 'G';
147 break;
148 case HShell:
149 ShellLabel = 'H';
150 break;
151 case IShell:
152 ShellLabel = 'I';
153 break;
154 }
155 sprintf(Line, " %c %d", ShellLabel, Shells[ishell].NumPrims);
156 File->WriteLine(Line, true);
157 for (iprim=0; iprim<Shells[ishell].NumPrims; iprim++) {
158 double EE = Shells[ishell].Exponent[iprim] * Shells[ishell].Exponent[iprim];
159 double FACS = PI32/(EE*sqrt(EE));
160 double FACP = 0.5*FACS/EE;
161 double FACD = 0.75*FACS/(EE*EE);
162 double FACF = 1.875*FACS/(EE*EE*EE);
163 double FACG = 6.5625*FACS/(EE*EE*EE*EE);
164 double FACH = 29.53125*FACS/(EE*EE*EE*EE*EE);
165 double FACI = 162.421875*FACS/(EE*EE*EE*EE*EE*EE);
166 float backconv, backl;
167 switch (Shells[ishell].ShellType) {
168 case LShell:
169 backl = (float) sqrt(FACP);
170 case SShell:
171 backconv = (float) sqrt(FACS);
172 break;
173 case PShell:
174 backconv = (float) sqrt(FACP);
175 break;
176 case DShell:
177 backconv = (float) sqrt(FACD);
178 break;
179 case FShell:
180 backconv = (float) sqrt(FACF);
181 break;
182 case GShell:
183 backconv = (float) sqrt(FACG);
184 break;
185 case HShell:
186 backconv = (float) sqrt(FACH);
187 break;
188 case IShell:
189 backconv = (float) sqrt(FACI);
190 break;
191 }
192 sprintf(Line, " %ld %f %f", iprim+1, Shells[ishell].Exponent[iprim],
193 Shells[ishell].NormCoef[iprim]*backconv);
194 File->WriteLine(Line, false);
195 if (Shells[ishell].ShellType == LShell) {
196 sprintf(Line, " %f", Shells[ishell].NormCoef[iprim+Shells[ishell].NumPrims]*backl);
197 File->WriteLine(Line, false);
198 }
199 File->WriteLine("", true);
200 }
201 }
202 } else {
203 File->WriteLine("Error! Basis not defined for this atom!", true);
204 }
205 File->WriteLine("",true);
206 }
207
BasisShell(void)208 BasisShell::BasisShell(void) {
209 NumPrims = 0;
210 }
~BasisShell(void)211 BasisShell::~BasisShell(void) {
212 }
GetNumFuncs(bool UseSphericalHarmonics) const213 long BasisShell::GetNumFuncs(bool UseSphericalHarmonics) const {
214 //returns the # of functions for the shell based on the shell type
215 long ret;
216 long type = ShellType;
217 if (UseSphericalHarmonics) type += 10;
218 switch (type) {
219 case LShell: //L shell
220 case SHLShell:
221 ret = 4;
222 break;
223 case SShell: //S shell
224 case SHSShell:
225 ret = 1;
226 break;
227 case PShell: //P shell
228 case SHPShell:
229 ret = 3;
230 break;
231 case DShell: //D shell
232 ret = 6;
233 break;
234 case FShell: //F shell
235 ret = 10;
236 break;
237 case GShell: //G shell
238 ret = 15;
239 break;
240 case HShell:
241 ret = 21;
242 break;
243 case IShell:
244 ret = 28;
245 break;
246 case SHDShell: //Spherical harmonic D
247 ret = 5;
248 break;
249 case SHFShell: //Spherical harmonic F
250 ret = 7;
251 break;
252 case SHGShell: //Spherical harmonic G
253 ret = 9;
254 break;
255 case SHHShell: //Spherical harmonic H
256 ret = 11;
257 break;
258 case SHIShell: //Spherical harmonic I
259 ret = 13;
260 break;
261 }
262 return ret;
263 }
264 //Return the start of the function in the angular momenta list (MINI)
GetAngularStart(bool UseSphericalHarmonics) const265 long BasisShell::GetAngularStart(bool UseSphericalHarmonics) const {
266 long ret;
267 long type = ShellType;
268 if (UseSphericalHarmonics) type += 10;
269 switch (type) {
270 case LShell: //L shell
271 case SHLShell:
272 ret = 0;
273 break;
274 case SShell: //S shell
275 case SHSShell:
276 ret = 0;
277 break;
278 case PShell: //P shell
279 case SHPShell:
280 ret = 1;
281 break;
282 case DShell: //D shell
283 case SHDShell:
284 ret = 4;
285 break;
286 case FShell: //F shell
287 ret = 10;
288 break;
289 case GShell: //G shell
290 ret = 20;
291 break;
292 case HShell:
293 ret = 35;
294 break;
295 case IShell:
296 ret = 56;
297 break;
298 case SHFShell:
299 ret = 9;
300 break;
301 case SHGShell:
302 ret = 16;
303 break;
304 case SHHShell:
305 ret = 25;
306 break;
307 case SHIShell:
308 ret = 36;
309 break;
310 }
311 return ret;
312 }
GetLabel(char * label,short FuncNum,bool UseSphericalHarmonics) const313 void BasisShell::GetLabel(char * label, short FuncNum, bool UseSphericalHarmonics) const {
314 long type = ShellType;
315 if (UseSphericalHarmonics) type += 10;
316 switch (type) {
317 case LShell:
318 case SShell:
319 case SHLShell:
320 case SHSShell:
321 if (FuncNum==0) { //s or s part of an L shell
322 label[0]='S';
323 label[1]=0;
324 break;
325 } else FuncNum--;
326 case PShell:
327 case SHPShell:
328 label[0] = 'P';
329 switch (FuncNum) {
330 case 0: //Px
331 label[1]='x';
332 break;
333 case 1: //Py
334 label[1]='y';
335 break;
336 case 2: //Pz
337 label[1]='z';
338 break;
339 }
340 label[2]=0;
341 break;
342 case DShell: //D shell
343 label[0]='D';
344 switch (FuncNum) {
345 case 0: //Dxx
346 label[1]='x';
347 label[2]='2';
348 break;
349 case 1: //Dyy
350 label[1]='y';
351 label[2]='2';
352 break;
353 case 2: //Dzz
354 label[1]='z';
355 label[2]='2';
356 break;
357 case 3: //Dxy
358 label[1]='x';
359 label[2]='y';
360 break;
361 case 4: //Dxz
362 label[1]='x';
363 label[2]='z';
364 break;
365 case 5: //Dyz
366 label[1]='y';
367 label[2]='z';
368 break;
369 }
370 label[3]=0;
371 break;
372 case FShell: //F shell
373 label[0]='F';
374 switch (FuncNum) {
375 case 0: //Fxxx
376 label[1]='x';
377 label[2]='3';
378 label[3]=0;
379 break;
380 case 1: //Fyyy
381 label[1]='y';
382 label[2]='3';
383 label[3]=0;
384 break;
385 case 2: //Fzzz
386 label[1]='z';
387 label[2]='3';
388 label[3]=0;
389 break;
390 case 3: //Fxxy
391 label[1]='x';
392 label[2]='2';
393 label[3]='y';
394 break;
395 case 4: //Fxxz
396 label[1]='x';
397 label[2]='2';
398 label[3]='z';
399 break;
400 case 5: //Fxyy
401 label[1]='x';
402 label[2]='y';
403 label[3]='2';
404 break;
405 case 6: //Fyyz
406 label[1]='y';
407 label[2]='2';
408 label[3]='z';
409 break;
410 case 7: //Fxzz
411 label[1]='x';
412 label[2]='z';
413 label[3]='2';
414 break;
415 case 8: //Fyzz
416 label[1]='y';
417 label[2]='z';
418 label[3]='2';
419 break;
420 case 9: //Fxyz
421 label[1]='x';
422 label[2]='y';
423 label[3]='z';
424 break;
425 }
426 label[4]=0;
427 break;
428 case GShell: //G shell
429 label[0]='G';
430 switch (FuncNum) {
431 case 0: //Gx4
432 label[1]='x';
433 label[2]='4';
434 label[3]=0;
435 break;
436 case 1: //Gy4
437 label[1]='y';
438 label[2]='4';
439 label[3]=0;
440 break;
441 case 2: //Gz4
442 label[1]='z';
443 label[2]='4';
444 label[3]=0;
445 break;
446 case 3: //Gx3y
447 label[1]='x';
448 label[2]='3';
449 label[3]='y';
450 label[4]=0;
451 break;
452 case 4: //Gx3z
453 label[1]='x';
454 label[2]='3';
455 label[3]='z';
456 label[4]=0;
457 break;
458 case 5: //Gxy3
459 label[1]='x';
460 label[2]='y';
461 label[3]='3';
462 label[4]=0;
463 break;
464 case 6: //Gy3z
465 label[1]='y';
466 label[2]='3';
467 label[3]='z';
468 label[4]=0;
469 break;
470 case 7: //Gxz3
471 label[1]='x';
472 label[2]='z';
473 label[3]='3';
474 label[4]=0;
475 break;
476 case 8: //Gyz3
477 label[1]='y';
478 label[2]='z';
479 label[3]='3';
480 label[4]=0;
481 break;
482 case 9: //Gx2y2
483 label[1]='x';
484 label[2]='2';
485 label[3]='y';
486 label[4]='2';
487 break;
488 case 10://Gx2z2
489 label[1]='x';
490 label[2]='2';
491 label[3]='z';
492 label[4]='2';
493 break;
494 case 11://Gy2z2
495 label[1]='y';
496 label[2]='2';
497 label[3]='z';
498 label[4]='2';
499 break;
500 case 12://Gx2yz
501 label[1]='x';
502 label[2]='2';
503 label[3]='y';
504 label[4]='z';
505 break;
506 case 13://Gxy2z
507 label[1]='x';
508 label[2]='y';
509 label[3]='2';
510 label[4]='z';
511 break;
512 case 14://Gxyz2
513 label[1]='x';
514 label[2]='y';
515 label[3]='z';
516 label[4]='2';
517 break;
518 }
519 label[5]=0;
520 break;
521 case HShell:
522 switch (FuncNum) {
523 case 0:
524 strcpy(label, "Hx5");
525 break;
526 case 1:
527 strcpy(label, "Hy5");
528 break;
529 case 2:
530 strcpy(label, "Hz5");
531 break;
532 case 3:
533 strcpy(label, "Hx4y");
534 break;
535 case 4:
536 strcpy(label, "Hx4z");
537 break;
538 case 5:
539 strcpy(label, "Hxy4");
540 break;
541 case 6:
542 strcpy(label, "Hy4z");
543 break;
544 case 7:
545 strcpy(label, "Hxz4");
546 break;
547 case 8:
548 strcpy(label, "Hyz4");
549 break;
550 case 9:
551 strcpy(label, "Hx3y2");
552 break;
553 case 10:
554 strcpy(label, "Hx3z2");
555 break;
556 case 11:
557 strcpy(label, "Hx2y3");
558 break;
559 case 12:
560 strcpy(label, "Hy3z2");
561 break;
562 case 13:
563 strcpy(label, "Hx2z3");
564 break;
565 case 14:
566 strcpy(label, "Hy2z3");
567 break;
568 case 15:
569 strcpy(label, "Hx3yz");
570 break;
571 case 16:
572 strcpy(label, "Hxy3z");
573 break;
574 case 17:
575 strcpy(label, "Hxyz3");
576 break;
577 case 18:
578 strcpy(label, "Hx2y2z");
579 break;
580 case 19:
581 strcpy(label, "Hx2yz2");
582 break;
583 case 20:
584 strcpy(label, "Hxy2z2");
585 break;
586 }
587 break;
588 case IShell:
589 switch (FuncNum) {
590 case 0:
591 strcpy(label, "Ix6");
592 break;
593 case 1:
594 strcpy(label, "Iy6");
595 break;
596 case 2:
597 strcpy(label, "Iz6");
598 break;
599 case 3:
600 strcpy(label, "Ix5y");
601 break;
602 case 4:
603 strcpy(label, "Ix5z");
604 break;
605 case 5:
606 strcpy(label, "Ixy5");
607 break;
608 case 6:
609 strcpy(label, "Iy5z");
610 break;
611 case 7:
612 strcpy(label, "Ixz5");
613 break;
614 case 8:
615 strcpy(label, "Iyz5");
616 break;
617 case 9:
618 strcpy(label, "Ix4y2");
619 break;
620 case 10:
621 strcpy(label, "Ix4z2");
622 break;
623 case 11:
624 strcpy(label, "Ix2y4");
625 break;
626 case 12:
627 strcpy(label, "Iy4z2");
628 break;
629 case 13:
630 strcpy(label, "Ix2z4");
631 break;
632 case 14:
633 strcpy(label, "Iy2z4");
634 break;
635 case 15:
636 strcpy(label, "Ix4yz");
637 break;
638 case 16:
639 strcpy(label, "Ixy4z");
640 break;
641 case 17:
642 strcpy(label, "Ixyz4");
643 break;
644 case 18:
645 strcpy(label, "Ix3y3");
646 break;
647 case 19:
648 strcpy(label, "Ix3z3");
649 break;
650 case 20:
651 strcpy(label, "Iy3z3");
652 break;
653 case 21:
654 strcpy(label, "Ix3y2z");
655 break;
656 case 22:
657 strcpy(label, "Ix3yz2");
658 break;
659 case 23:
660 strcpy(label, "Ix2y3z");
661 break;
662 case 24:
663 strcpy(label, "Ixy3z2");
664 break;
665 case 25:
666 strcpy(label, "Ix2yz3");
667 break;
668 case 26:
669 strcpy(label, "Ixy2z3");
670 break;
671 case 27:
672 strcpy(label, "Ix2y2z2");
673 break;
674 }
675 break;
676 case SHDShell:
677 label[0]='D';
678 switch (FuncNum) {
679 case 0:
680 strcpy((char *) &(label[1]), "xy");
681 break;
682 case 1:
683 strcpy((char *) &(label[1]), "yz");
684 break;
685 case 2:
686 strcpy((char *) &(label[1]), "3z2-r2");
687 break;
688 case 3:
689 strcpy((char *) &(label[1]), "xz");
690 break;
691 case 4:
692 strcpy((char *) &(label[1]), "x2-y2");
693 break;
694 }
695 break;
696 case SHFShell:
697 label[0]='F';
698 switch (FuncNum) {
699 case 0:
700 strcpy((char *) &(label[1]), "x3-3xy2");
701 break;
702 case 1:
703 strcpy((char *) &(label[1]), "x2z-y2z");
704 break;
705 case 2:
706 strcpy((char *) &(label[1]), "x(5z2-r2)");
707 break;
708 case 3:
709 strcpy((char *) &(label[1]), "z(5z2-3r2)");
710 break;
711 case 4:
712 strcpy((char *) &(label[1]), "y(5z2-r2)");
713 break;
714 case 5:
715 strcpy((char *) &(label[1]), "xyz");
716 break;
717 case 6:
718 strcpy((char *) &(label[1]), "3x2y-y3");
719 break;
720 }
721 break;
722 case SHGShell:
723 label[0]='G';
724 switch (FuncNum) {
725 case 0:
726 strcpy((char *) &(label[1]), "x4+y4-6x2y2");
727 break;
728 case 1:
729 strcpy((char *) &(label[1]), "xz(x2-3y2)");
730 break;
731 case 2:
732 strcpy((char *) &(label[1]), "(x2-y2)(7z2-r2)");
733 break;
734 case 3:
735 strcpy((char *) &(label[1]), "xz(7z2-3r2)");
736 break;
737 case 4:
738 strcpy((char *) &(label[1]), "35z4-30z2r2+3r4");
739 break;
740 case 5:
741 strcpy((char *) &(label[1]), "yz(7z2-3r2)");
742 break;
743 case 6:
744 strcpy((char *) &(label[1]), "xy(7z2-r2)");
745 break;
746 case 7:
747 strcpy((char *) &(label[1]), "yz(3x2-y2)");
748 break;
749 case 8:
750 strcpy((char *) &(label[1]), "xy(x2-y2)");
751 break;
752 }
753 break;
754 case SHHShell:
755 label[0]='H';
756 switch (FuncNum) {
757 case 0:
758 strcpy((char *) &(label[1]), "x5-10x3y2+5xy4");
759 break;
760 case 1:
761 strcpy((char *) &(label[1]), "z(x4+y4-6x2y2)");
762 break;
763 case 2:
764 strcpy((char *) &(label[1]), "(x3-3xy2)(9z2-r2)");
765 break;
766 case 3:
767 strcpy((char *) &(label[1]), "(x2-y2)(3z3-zr2)");
768 break;
769 case 4:
770 strcpy((char *) &(label[1]), "21xz4-14xz2r2+xr4");
771 break;
772 case 5:
773 strcpy((char *) &(label[1]), "63z5-70z3r2+15zr4");
774 break;
775 case 6:
776 strcpy((char *) &(label[1]), "21yz4-14yz2r2+yr4");
777 break;
778 case 7:
779 strcpy((char *) &(label[1]), "xyz(3z3-zr2)");
780 break;
781 case 8:
782 strcpy((char *) &(label[1]), "(3x2y-y3)(9z2-r2)");
783 break;
784 case 9:
785 strcpy((char *) &(label[1]), "z(x3y-xy3)");
786 break;
787 case 10:
788 strcpy((char *) &(label[1]), "5x4y-10x2y3+5y5");
789 break;
790 }
791 break;
792 case SHIShell:
793 label[0]='I';
794 switch (FuncNum) {
795 case 0:
796 strcpy((char *) &(label[1]), "x6-15x4y2+15x2y4-y6");
797 break;
798 case 1:
799 strcpy((char *) &(label[1]), "xz(x4+5y4-10x2y2)");
800 break;
801 case 2:
802 strcpy((char *) &(label[1]), "(x4+y4-6x2y2)(11z2-r2)");
803 break;
804 case 3:
805 strcpy((char *) &(label[1]), "(x3-3xy2)z(11z2-3r2)");
806 break;
807 case 4:
808 strcpy((char *) &(label[1]), "(x2-y2)(33z4-18z2r2+r4)");
809 break;
810 case 5:
811 strcpy((char *) &(label[1]), "xz(33z4-30z2r2+5r4");
812 break;
813 case 6:
814 strcpy((char *) &(label[1]), "231z6-315z4r2+105z2r4-5r6");
815 break;
816 case 7:
817 strcpy((char *) &(label[1]), "yz(33z4-30z2r2+5r4");
818 break;
819 case 8:
820 strcpy((char *) &(label[1]), "xy(33z4-18z2r2+r4)");
821 break;
822 case 9:
823 strcpy((char *) &(label[1]), "z(11z2-3r2)(3x2y-y3)");
824 break;
825 case 10:
826 strcpy((char *) &(label[1]), "xy(x2-y2)(11z2-r2)");
827 break;
828 case 11:
829 strcpy((char *) &(label[1]), "yz(y4+5x4-10x2y2)");
830 break;
831 case 12:
832 strcpy((char *) &(label[1]), "xy(3x4-10x2y2+3y4");
833 break;
834 }
835 break;
836 default:
837 wxLogMessage(_("Undefined shell type not implemented in BasisShell::GetLabel"));
838 break;
839 }
840 }
841 //Read in a general basis set from a GAMESS log file
ParseGAMESSBasisSet(BufferFile * Buffer,long NumAtoms,const mpAtom * Atoms)842 BasisSet * BasisSet::ParseGAMESSBasisSet(BufferFile * Buffer, long NumAtoms, const mpAtom * Atoms) {
843 long NumShells=0, NextAtom=0, ShellStart=0,
844 fShellStart=0, BasisEndPos, nTotFuncs=0, LinePos, EndPos, iatom, items;
845 char LineText[kMaxLineLength+1];
846 Buffer->GetLine(LineText);
847
848 Buffer->LocateKeyWord("CONTRACTION COEFFICIENT", 21);
849 Buffer->SkipnLines(2);
850 long StartPos = Buffer->GetFilePos();
851 //go to the end of the basis and get the total # of shells to dimension mem
852 //output was changed in June, 2000
853 bool ShellsFound = Buffer->LocateKeyWord("TOTAL NUMBER OF BASIS SET SHELLS", 32);
854 if (!ShellsFound) ShellsFound = Buffer->LocateKeyWord("TOTAL NUMBER OF SHELLS", 22);
855 if (ShellsFound) {
856 BasisEndPos = Buffer->GetFilePos();
857 Buffer->GetLine(LineText);
858 LinePos = FindKeyWord(LineText, "=", 1) + 1;
859 sscanf(&(LineText[LinePos]),"%ld",&NumShells);
860 } else return NULL;
861 if (NumShells<=0) return NULL; //invalid number of shells
862 if (NumAtoms <= 0) return NULL;
863 BasisSet * Basis = new BasisSet(NumAtoms, NumShells);
864 if (Basis == NULL) return NULL;//No memory!
865
866 Buffer->SetFilePos(StartPos);
867 Buffer->SkipnLines(2); //advance to the first primitive line
868 //Note the atom label is skipped since it is just a general text string and the
869 //Basis set is always given in the same order as the coordinates
870 //The ishell loop could be a do-while structure, but this ensures
871 //that I don't go off the end of the shell array
872 for (long ishell=0; ishell<NumShells; ishell++) {
873 float fExp, fCoef1, fCoef2, fCoef3, scrap;
874 char ShellText[10];
875 long fShell, scrapl, NumPrims, CoefsperPrim;
876 short ShellType;
877 BasisShell foo;
878
879 StartPos = Buffer->GetFilePos();
880 EndPos = Buffer->FindBlankLine();
881 NumPrims = Buffer->GetNumLines(EndPos-StartPos);
882 Buffer->GetLine(LineText);
883 items = sscanf(LineText, "%ld%5s",&fShell, ShellText);
884 if (items != 2) throw DataError();
885 CoefsperPrim = 1;
886 switch (ShellText[0]) {
887 case 'S':
888 ShellType=SShell;
889 break;
890 case 'P':
891 ShellType=PShell;
892 break;
893 case 'D':
894 ShellType=DShell;
895 break;
896 case 'F':
897 ShellType=FShell;
898 break;
899 case 'G':
900 ShellType=GShell;
901 break;
902 case 'H':
903 ShellType=HShell;
904 break;
905 case 'I':
906 ShellType=IShell;
907 break;
908 case 'L':
909 ShellType=LShell;
910 CoefsperPrim = 2;
911 break;
912 default:
913 wxString msg;
914 msg.Printf(_("Unknown shell type encountered, shell label = %s"), ShellText);
915 wxLogMessage(msg);
916 throw DataError();
917 }
918 Basis->Shells.push_back(foo);
919 Basis->Shells[ishell].NumPrims = (short) NumPrims;
920 Basis->Shells[ishell].ShellType = ShellType;
921 std::vector<float> & Exponent = Basis->Shells[ishell].Exponent;
922 std::vector<float> & NormCoef = Basis->Shells[ishell].NormCoef;
923 std::vector<float> & InputCoef = Basis->Shells[ishell].InputCoef;
924 Exponent.reserve(NumPrims);
925 if (ShellType >= 0) {
926 NormCoef.reserve(NumPrims);
927 InputCoef.reserve(NumPrims);
928 } else {
929 NormCoef.reserve(2*NumPrims);
930 InputCoef.reserve(2*NumPrims);
931 }
932 bool OldStyle = false;
933 if (FindKeyWord(LineText, "(", 1) != -1) OldStyle=true;
934 for (long iprim=0; iprim<NumPrims; iprim++) {
935 InputCoef.push_back(0.0);
936 NormCoef.push_back(0.0);
937 if (ShellType<0) {
938 InputCoef.push_back(0.0);
939 NormCoef.push_back(0.0);
940 }
941 Exponent.push_back(0.0);
942 }
943 for (long iprim=0; iprim<NumPrims; iprim++) {
944 if (ShellType>=0) {
945 if (OldStyle) {
946 items = sscanf(LineText, "%ld%5s%ld%f%f (%f)",&fShell, ShellText,
947 &scrapl, &fExp, &fCoef1, &scrap);
948 if (items != 6) throw DataError();
949 InputCoef[iprim] = scrap;
950 } else {
951 items = sscanf(LineText, "%ld%5s%ld%f%f",&fShell, ShellText,
952 &scrapl, &fExp, &fCoef1);
953 if (items != 5) throw DataError();
954 InputCoef[iprim] = fCoef1;
955 }
956 } else {
957 if (OldStyle) {
958 items = sscanf(LineText, "%ld%5s%ld%f%f (%f)%f (%f)",&fShell, ShellText,
959 &scrapl, &fExp, &fCoef1,&scrap,&fCoef2,&fCoef3);
960 NormCoef[iprim+NumPrims] = fCoef2;
961 InputCoef[iprim+NumPrims] = fCoef3;
962 InputCoef[iprim] = scrap;
963 if (items != 8) throw DataError();
964 } else {
965 items = sscanf(LineText, "%ld%5s%ld%f%f %f",&fShell, ShellText,
966 &scrapl, &fExp, &fCoef1,&fCoef2);
967 NormCoef[iprim+NumPrims] = fCoef2; //This is not normalized, but we take care of that later
968 InputCoef[iprim+NumPrims] = fCoef2;
969 InputCoef[iprim] = fCoef1;
970 if (items != 6) throw DataError();
971 }
972 }
973 Exponent[iprim] = fExp;
974 NormCoef[iprim] = fCoef1;
975 //Coefficients are already normalized when read in from in front of parenthesis
976 //so this code is unneeded but might be useful in the future
977 /* ExpNorm = sqrt(sqrt(fExp));
978 ExpNorm2 = ExpNorm*ExpNorm;
979 ExpNorm *= ExpNorm2;
980 //Normalization constant is: [fExp^(n+3/2) 2^(2n+3/2)/(pi^3/2 n!!)]^1/2
981 //This factor is correct for functions of the form xx, xxx, etc but is off
982 //for the xy, xxy, etc functions which must be multiplied by sqrt3, sqrt5 etc later
983 switch (ShellType) {
984 case -1: //L shells
985 NormCoef[iprim+NumPrims] = fCoef2*PNorm*ExpNorm*ExpNorm2;
986 case 0: //S shell
987 NormCoef[iprim] *= SNorm*ExpNorm;
988 break;
989 case 1: //P shell
990 NormCoef[iprim] *= PNorm*ExpNorm*ExpNorm2;
991 break;
992 case 2: //D shell
993 NormCoef[iprim] *= DNorm*fExp*ExpNorm;
994 break;
995 case 3: //F shell
996 NormCoef[iprim] *= FNorm*fExp*ExpNorm*ExpNorm2;
997 break;
998 case 4: //G shell
999 NormCoef[iprim] *= GNorm*fExp*fExp*ExpNorm;
1000 break;
1001 }*/
1002 Buffer->GetLine(LineText);
1003 }
1004 StartPos = Buffer->GetFilePos();
1005 Buffer->GetLine(LineText);
1006 Buffer->SetFilePos(StartPos);
1007 long iscanerr = sscanf(LineText, "%ld %5s %ld %f (%f) ",&scrapl, ShellText,
1008 &scrapl, &fExp, &scrap);
1009 nTotFuncs += Basis->Shells[ishell].GetNumFuncs(false);
1010 Basis->NumShells++;
1011 if (iscanerr<3) {//Found next atom label so tidy up the current atom
1012 Buffer->SkipnLines(2); //Skip Atom label and blank line
1013 long nGotShells = ishell-ShellStart+1;
1014 long nMappedAtoms = (fShell-fShellStart-nGotShells)/nGotShells + 1;
1015 Basis->NumFuncs += nMappedAtoms*nTotFuncs;
1016 //skip over any MM atoms since they aren't included in the basis
1017 while (Atoms[NextAtom].IsSIMOMMAtom() && (NextAtom < NumAtoms)) NextAtom++;
1018 nTotFuncs = 0; //reset the function counter for the next atom
1019 for (iatom=0; iatom<nMappedAtoms; iatom++) {
1020 Basis->BasisMap[2*NextAtom] = ShellStart;
1021 Basis->BasisMap[2*NextAtom+1] = ishell;
1022 NextAtom++;
1023 }
1024 ShellStart = ishell+1;
1025 fShellStart = fShell;
1026 }
1027 if (fShell >= NumShells) break;
1028 }
1029 return Basis;
1030 } /*ParseGAMESSBasisSet*/
1031
1032 #include <iostream>
ReadMolDenBasisSet(BufferFile * Buffer,long NumAtoms)1033 bool BasisSet::ReadMolDenBasisSet(BufferFile * Buffer, long NumAtoms) {
1034 bool result = true;
1035 // bool use5D=false, use7F=false, use10G=false;
1036 char LineText[kMaxLineLength];
1037 long shellCount = 0;
1038
1039 Buffer->SkipnLines(1);
1040 for (int i=0; i<NumAtoms; i++) {
1041 long shellStart = shellCount;
1042 long atomNum=-1;
1043 Buffer->GetLine(LineText);
1044 sscanf(LineText, "%ld", &atomNum);
1045 if (atomNum != (i+1)) return false; //incorrect order? Need to log a message...
1046 Buffer->GetLine(LineText);
1047 while (!IsBlank(LineText)) {
1048 char token[kMaxLineLength];
1049 int numPrims;
1050 if (sscanf(LineText, "%s %d", token, &numPrims) == 2) {
1051 short ShellType;
1052 if (strlen(token) == 1) {
1053 switch (toupper(token[0])) {
1054 case 'S':
1055 ShellType = SShell;
1056 break;
1057 case 'P':
1058 ShellType = PShell;
1059 break;
1060 case 'D':
1061 ShellType = DShell;
1062 break;
1063 case 'F':
1064 ShellType = FShell;
1065 break;
1066 case 'G':
1067 ShellType = GShell;
1068 break;
1069 case 'H':
1070 ShellType = HShell;
1071 break;
1072 case 'I':
1073 ShellType = IShell;
1074 break;
1075 default:
1076 wxString msg;
1077 msg.Printf(_("Unknown shell type encountered, shell label = %s"), token);
1078 wxLogMessage(msg);
1079 return false; //invalid type
1080 }
1081 } else if (strlen(token) == 2) {
1082 if ((toupper(token[0])=='S')&&(toupper(token[1])=='P'))
1083 ShellType = -1;
1084 else return false;
1085 } else return false; //invalid shell type
1086 if (numPrims <= 0) return false;
1087
1088 BasisShell t;
1089 Shells.push_back(t);
1090 Shells[shellCount].NumPrims = (short) numPrims;
1091 Shells[shellCount].ShellType = ShellType;
1092 std::vector<float> & Exponent = Shells[shellCount].Exponent;
1093 std::vector<float> & NormCoef = Shells[shellCount].NormCoef;
1094 std::vector<float> & InputCoef = Shells[shellCount].InputCoef;
1095 Exponent.reserve(numPrims);
1096 if (ShellType >= 0) {
1097 NormCoef.reserve(numPrims);
1098 InputCoef.reserve(numPrims);
1099 } else {
1100 NormCoef.reserve(2*numPrims);
1101 InputCoef.reserve(2*numPrims);
1102 }
1103 for (long iprim=0; iprim<numPrims; iprim++) {
1104 InputCoef.push_back(0.0);
1105 NormCoef.push_back(0.0);
1106 if (ShellType<0) {
1107 InputCoef.push_back(0.0);
1108 NormCoef.push_back(0.0);
1109 }
1110 Exponent.push_back(0.0);
1111 }
1112 for (long iprim=0; iprim<numPrims; iprim++) {
1113 float fexp, fcoef, fcoefp;
1114 Buffer->GetLine(LineText);
1115 ConvertExponentStyle(LineText);
1116 if (ShellType >= 0) {
1117 if (sscanf(LineText, "%f %f", &fexp, &fcoef) != 2) return false;
1118 } else {
1119 if (sscanf(LineText, "%f %f %f", &fexp, &fcoef, &fcoefp) == 3) {
1120 InputCoef[iprim+numPrims] = fcoefp;
1121 } else return false;
1122 }
1123 Exponent[iprim] = fexp;
1124 InputCoef[iprim] = fcoef;
1125 }
1126 NumFuncs += Shells[shellCount].GetNumFuncs(false);
1127 NumShells ++;
1128 shellCount ++;
1129 } else return false;
1130 Buffer->GetLine(LineText);
1131 }
1132 BasisMap[2*i] = shellStart;
1133 BasisMap[2*i+1] = shellCount - 1;
1134 }
1135 return result;
1136 }
1137