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