1 /* Basis.c */
2 /**********************************************************************************************************
3 Copyright (c) 2002-2013 Abdul-Rahman Allouche. All rights reserved
4
5 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
6 documentation files (the Gabedit), to deal in the Software without restriction, including without limitation
7 the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
8 and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
9
10 The above copyright notice and this permission notice shall be included in all copies or substantial portions
11 of the Software.
12
13 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
14 TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
15 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
16 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
17 DEALINGS IN THE SOFTWARE.
18 ************************************************************************************************************/
19
20 #include "../../Config.h"
21 #include "GlobalOrb.h"
22 #include "UtilsOrb.h"
23 #include "../Utils/Utils.h"
24 #include "../Utils/UtilsInterface.h"
25 #include "../Utils/Constants.h"
26 #include "../Utils/Zlm.h"
27 #include "../Utils/GTF.h"
28
29 /********************************************************************************/
save_basis_gabedit_format(FILE * file)30 void save_basis_gabedit_format(FILE* file)
31 {
32 gint i;
33 gint j;
34 gint k;
35 gint c;
36 fprintf(file,"[Basis]\n");
37 if(!Type) return;
38 for(c = 0;c<nCenters; c++)
39 {
40 i = GeomOrb[c].NumType;
41 fprintf(file," %d 0\n", c+1);
42 for(j=0;j<Type[i].Norb;j++)
43 {
44 fprintf(file," %c %d 1.0\n",GetSymmetry(Type[i].Ao[j].L),Type[i].Ao[j].N);
45 for(k=0;k<Type[i].Ao[j].N;k++)
46 fprintf(file," %14.6f %14.6f\n",Type[i].Ao[j].Ex[k],Type[i].Ao[j].Coef[k]);
47 }
48 fprintf(file,"\n");
49 }
50 fprintf(file,"\n");
51 }
52 /**********************************************/
ReadCommandLines()53 gint ReadCommandLines()
54 {
55 gint taille=BSIZE;
56 char *t = g_malloc(taille);
57
58 if(forb == NULL)
59 return 0;
60 while(!feof(forb))
61 {
62 { char* e = fgets(t,taille,forb);}
63 if(!strcmp(t," ") || !strcmp(t,"\n") )
64 return 1;
65 if( t[0] != '#' )
66 {
67 Debug("\n\nERROR : theyre firsts lines is not a commads lines\n\n");
68 return 0;
69 }
70 }
71 return 0;
72 }
73 /**********************************************/
OpenDataFile(char * NameFile)74 gint OpenDataFile(char * NameFile)
75 {
76 forb = FOpen(NameFile, "rb");
77 if(forb == NULL)
78 {
79 Debug("\n\nERROR Can not open %s data file\n\n",NameFile);
80 return 0;
81 }
82 return 1;
83 }
84 /**********************************************/
PrintBasis()85 void PrintBasis()
86 {
87 gint i;
88 gint j,k;
89 Debug("\n\n");
90 printLineChar('*',72);
91 for(i=0;i<Ntype;i++)
92 {
93 Debug("Basis for center %s \n",Type[i].Symb);
94 Debug("====================\n");
95 Debug("Norb\t\t= %d\n",Type[i].Norb);
96 Debug("\n");
97 Debug("- --------- --------- \n");
98 Debug("l Ex Coef \n");
99 Debug("- --------- --------- \n");
100 for(j=0;j<Type[i].Norb;j++)
101 { k=0;
102 Debug("%c %14.6f %14.6f\n",GetSymmetry(Type[i].Ao[j].L),Type[i].Ao[j].Ex[k],Type[i].Ao[j].Coef[k]);
103 for(k=1;k<Type[i].Ao[j].N;k++)
104 Debug(" %14.6f %14.6f\n",Type[i].Ao[j].Ex[k],Type[i].Ao[j].Coef[k]);
105
106 }
107 Debug("\n");
108 /*
109 if(Type[i].nps>0)
110 {
111 Debug("Pseudo potential for center %s\n",Type[i].Symb);
112 Debug("==============================\n\n");
113 Debug("-\t- --------- ---------\n");
114 Debug("l\tn Eps Tau\n");
115 Debug("-\t- --------- ---------\n");
116 }
117 else
118 Debug("No Pseudo potential for %s\n",Type[i].Symb);
119 for(j=0;j<Type[i].nps;j++)
120 {
121 Debug("%d\t%d %9.6f %9.6f ",Type[i].ps[j].l,Type[i].ps[j].n,Type[i].ps[j].eps,Type[i].ps[j].tau);
122 Debug("\n");
123
124 }
125 */
126 printLineChar('*',72);
127 Debug("\n");
128 }
129 }
130
131 /**********************************************/
ReadOneBasis(gint i,gint j,char * t,gint * nsym)132 gboolean ReadOneBasis(gint i,gint j,char *t,gint *nsym)
133 {
134 gint k;
135 gint n;
136 gint taille = BSIZE;
137 gchar *sym = g_malloc(10*sizeof(gchar));
138 gint l=1;
139
140 /*Debug("i = %d j= %d\n",i,j);*/
141 /*Debug("t One %s\n",t);*/
142 /*sscanf(t,"%s %d",sym,&Type[i].Ao[j].N);*/
143 sscanf(t,"%s %d",sym,&n);
144 /* Debug("n = %d\n",n);*/
145 Type[i].Ao[j].N=n;
146 /*Debug("N = %d\n",Type[i].Ao[j].N);*/
147 Type[i].Ao[j].Ex=g_malloc(Type[i].Ao[j].N*sizeof(gdouble));
148 Type[i].Ao[j].Coef=g_malloc(Type[i].Ao[j].N*sizeof(gdouble));
149 /*Debug("avant sym ==\n");*/
150 if(strlen(sym)==2)
151 {
152 l=2;
153 Type[i].Ao=g_realloc(Type[i].Ao,(j+2)*sizeof(AO));
154 Type[i].Ao[j+1].N = Type[i].Ao[j].N;
155 Type[i].Ao[j+1].Ex=g_malloc(Type[i].Ao[j].N*sizeof(gdouble));
156 Type[i].Ao[j+1].Coef=g_malloc(Type[i].Ao[j].N*sizeof(gdouble));
157 }
158 *nsym = l;
159 /*Debug("nsym = %d\n",l);*/
160
161
162 for(k=0;k<Type[i].Ao[j].N;k++)
163 {
164 { char* e = fgets(t,taille,forb);}
165 for(n=0;n<(gint)strlen(t);n++)
166 if(t[n]=='D') t[n] = 'e';
167 /*Debug("t de One = %s\n",t);*/
168
169 if(l==1) sscanf(t,"%lf %lf ",&Type[i].Ao[j].Ex[k],&Type[i].Ao[j].Coef[k]);
170 else
171 {
172 sscanf(t,"%lf %lf %lf ",&Type[i].Ao[j].Ex[k],&Type[i].Ao[j].Coef[k],&Type[i].Ao[j+1].Coef[k]);
173 Type[i].Ao[j+1].Ex[k] = Type[i].Ao[j].Ex[k];
174 }
175 }
176 /*Debug("end k\n");*/
177
178 switch(sym[0])
179 {
180 case 's' :
181 case 'S' : Type[i].Ao[j].L=0;break;
182 case 'p' :
183 case 'P' : Type[i].Ao[j].L=1;break;
184 case 'd' :
185 case 'D' : Type[i].Ao[j].L=2;break;
186 case 'f' :
187 case 'F' : Type[i].Ao[j].L=3;break;
188 case 'g' :
189 case 'G' : Type[i].Ao[j].L=4;break;
190 case 'h' :
191 case 'H' : Type[i].Ao[j].L=5;break;
192 case 'i' :
193 case 'I' : Type[i].Ao[j].L=6;break;
194 case 'j' :
195 case 'J' : Type[i].Ao[j].L=7;break;
196 case 'k' :
197 case 'K' : Type[i].Ao[j].L=8;break;
198 case 'l' :
199 case 'L' : Type[i].Ao[j].L=9;break;
200 default : g_free(sym);
201 g_free(Type[i].Ao[j].Ex);
202 g_free(Type[i].Ao[j].Coef);
203 return FALSE;
204 }
205 /*
206 if(l == 2)
207 Debug(" sym[1] =%c \n",sym[1]);
208 */
209 if(l == 2)
210 switch(sym[1])
211 {
212 case 's' :
213 case 'S' : Type[i].Ao[j+1].L=0;break;
214 case 'p' :
215 case 'P' : Type[i].Ao[j+1].L=1;break;
216 case 'd' :
217 case 'D' : Type[i].Ao[j+1].L=2;break;
218 case 'f' :
219 case 'F' : Type[i].Ao[j+1].L=3;break;
220 case 'g' :
221 case 'G' : Type[i].Ao[j+1].L=4;break;
222 case 'h' :
223 case 'H' : Type[i].Ao[j+1].L=5;break;
224 case 'i' :
225 case 'I' : Type[i].Ao[j+1].L=6;break;
226 case 'j' :
227 case 'J' : Type[i].Ao[j+1].L=7;break;
228 case 'k' :
229 case 'K' : Type[i].Ao[j+1].L=8;break;
230 case 'l' :
231 case 'L' : Type[i].Ao[j+1].L=9;break;
232
233 default : g_free(sym);
234 g_free(Type[i].Ao[j+1].Ex);
235 g_free(Type[i].Ao[j+1].Coef);
236 return FALSE;
237 }
238 /*Debug("end readone basis \n");*/
239 g_free(sym);
240 return TRUE;
241 }
242 /**********************************************/
DefineBasisType(gchar * fileName)243 gboolean DefineBasisType(gchar *fileName)
244 {
245 gchar *sym;
246 gchar *t;
247 gchar *pdest;
248 gint taille=BSIZE;
249 gint i;
250 gint j;
251 gboolean ok;
252 gint nsym;
253 long int geompos = 0;
254
255 /* Debug("debut de DefineBasisType\n");*/
256 if ((!fileName) || (strcmp(fileName,"") == 0))
257 {
258 Message("Sorry\n No file selected","Error",TRUE);
259 return FALSE;
260 }
261
262 t=g_malloc(taille*sizeof(gchar));
263 forb = FOpen(fileName, "rb");
264 if(forb == NULL)
265 {
266 gchar buffer[BSIZE];
267 sprintf(buffer,"Sorry, I can not open '%s' file\n",fileName);
268 Message(buffer,"Error",TRUE);
269 return FALSE;
270 }
271 ok = FALSE;
272 while(!feof(forb))
273 {
274 { char* e = fgets(t,taille,forb);}
275 pdest = strstr(t,"asis set"); /* Basis for g98 and basis for g03 */
276 if(pdest != NULL)
277 {
278 ok = TRUE;
279 geompos = ftell(forb);
280 /* break; NO for get the last basis */
281 }
282 }
283 if(!ok)
284 {
285 gchar* t = g_strdup_printf("Sorry,\nI can not read atomic basis and/or molecular orbitals from\n \"%s\"\n file\n"
286 "\nFor display density and/or molecular orbitals from gaussian output file,\n"
287 "the following keywords are required:\n"
288 "#P Gfinput IOP(6/7=3) Pop=full\n"
289 "\n",
290 fileName);
291 Message(t,"Error",TRUE);
292 g_free(t);
293 return FALSE;
294 }
295 fseek(forb, geompos, SEEK_SET);
296
297 sym=g_malloc(10*sizeof(gchar));
298
299 t=(gchar*)g_malloc(taille*sizeof(gchar));
300
301 if(forb !=NULL)
302 {
303 /* Debug("Ntype = %d\n",Ntype);*/
304 Type = g_malloc(Ntype*sizeof(TYPE));
305 for(i=0;i<Ntype;i++)
306 Type[i].Ao = NULL;
307 while(!feof(forb))
308 {
309 { char* e = fgets(t,taille,forb);}
310 if(!strcmp(t,"\n") || !strcmp(t," ") )
311 {
312 g_free(sym);
313 g_free(t);
314 for(i=0;i<Ntype;i++)
315 if(Type[i].Ao == NULL)
316 {
317 gchar buffer[BSIZE];
318 sprintf(buffer,"Sorry, I can not read '%s' file, problem with basis set\n",fileName);
319 Message(buffer,"Error",TRUE);
320 return FALSE;
321 }
322 return TRUE;
323 }
324
325 if( strcmp(t,"\n")!= 0 && strcmp(t," ") && t[1] !='*')
326 {
327 sscanf(t,"%s",sym);
328 sym[0]=toupper(sym[0]);
329 if(strlen(sym)>1)
330 sym[1]=tolower(sym[1]);
331 /*i=GetNumType(sym);*/
332 i=atoi(sym)-1;
333 if(i>-1)
334 {
335 sym = g_strdup( GeomOrb[i].Symb);
336 i = GeomOrb[i].NumType;
337 Type[i].Symb=g_strdup(sym);
338 Type[i].N=GetNelectrons(sym);
339 j=-1;
340 while(!feof(forb))
341 {
342 { char* e = fgets(t,taille,forb);}
343 if(!strcmp(t,"\n") || t[1]=='*') break;
344 j++;
345 Type[i].Norb=j+1;
346 if(j == 0) Type[i].Ao=g_malloc(sizeof(AO));
347 else Type[i].Ao=g_realloc(Type[i].Ao,(j+1)*sizeof(AO));
348
349 ok= ReadOneBasis(i,j,t,&nsym);
350 if(nsym==2)
351 {
352 Type[i].Norb=j+2;
353 j++;
354 }
355 if(!ok)
356 {
357 j = j - nsym;
358 if(j==0) Type[i].Ao=g_malloc(sizeof(AO));
359 else Type[i].Ao=g_realloc(Type[i].Ao,(j+1)*sizeof(AO));
360 Type[i].Norb=j+1;
361 }
362 }
363 }
364 else
365 {
366 if(!strcmp(t,"\n") || !strcmp(t," ") )
367 {
368 g_free(sym);
369 g_free(t);
370 for(i=0;i<Ntype;i++)
371 if(Type[i].Ao == NULL)
372 {
373 gchar buffer[BSIZE];
374 sprintf(buffer,"Sorry, I can not read '%s' file, problem with basis set\n",fileName);
375 Message(buffer,"Error",TRUE);
376 return FALSE;
377 }
378 return TRUE;
379 }
380 break;
381 }
382 }
383 }
384 }
385 for(i=0;i<Ntype;i++)
386 if(Type[i].Ao == NULL)
387 {
388 gchar buffer[BSIZE];
389 sprintf(buffer,"Sorry, I can not read '%s' file, problem with basis set\n",fileName);
390 Message(buffer,"Error",TRUE);
391 return FALSE;
392 }
393 return TRUE;
394 }
395 /**********************************************/
resortAtoms(gint * numAtoms)396 static void resortAtoms(gint* numAtoms)
397 {
398 TypeGeomOrb* newGeom = NULL;
399 gint i;
400 /* printf("begin resortAtoms\n ");*/
401 if(!numAtoms)return;
402 if(nCenters<1)return;
403 /*
404 printf("Sorting ");
405 for(i=0;i<nCenters;i++) printf("%d ",numAtoms[i]);
406 printf("\n");
407 */
408 for(i=0;i<nCenters;i++) if(numAtoms[i] == -1) return;
409 newGeom = g_malloc(nCenters*sizeof(TypeGeomOrb));
410 for(i=0;i<nCenters;i++) newGeom[i] = GeomOrb[numAtoms[i]];
411 for(i=0;i<nCenters;i++) GeomOrb[i] = newGeom[i];
412 g_free(newGeom);
413 }
414 /**********************************************/
getNumberOfBasisCenters(gchar * fileName,gchar * title)415 static gint getNumberOfBasisCenters(gchar *fileName, gchar* title)
416 {
417 gchar t[BSIZE];
418 gint i;
419 gboolean ok;
420 gint nAtoms = 0;
421
422 if ((!fileName) || (strcmp(fileName,"") == 0)) return nAtoms;
423
424 forb = FOpen(fileName, "rb");
425 if(forb == NULL) return nAtoms;
426 ok = FALSE;
427 while(!feof(forb))
428 {
429 { char* e = fgets(t,BSIZE,forb);}
430 if(strstr(t,title) != NULL)
431 {
432 ok = TRUE;
433 break;
434 }
435 }
436 if(!ok)
437 {
438 fclose(forb);
439 return nAtoms;
440 }
441
442 nAtoms = 0;
443 while(!feof(forb))
444 {
445 { char* e = fgets(t,BSIZE,forb);}
446 if(this_is_a_backspace(t) || strstr(t,"[")) break;
447
448 /* printf("tt = %s\n",t);*/
449 i=atoi(t);
450 if(i>0)
451 {
452 nAtoms++;
453 while(!feof(forb))
454 {
455 { char* e = fgets(t,BSIZE,forb);}
456 if(this_is_a_backspace(t) || strstr(t,"[")) break;
457 }
458 }
459 else
460 {
461 if(this_is_a_backspace(t) || strstr(t,"[")) break;
462 }
463 if(strstr(t,"[")) break;
464 }
465 fclose(forb);
466 return nAtoms;
467 }
468 /**********************************************/
DefineGabeditMoldenBasisType(gchar * fileName,gchar * title)469 gboolean DefineGabeditMoldenBasisType(gchar *fileName,gchar* title)
470 {
471 gchar *sym;
472 gchar *t;
473 gchar *pdest;
474 gint taille=BSIZE;
475 gint i;
476 gint j;
477 gboolean ok;
478 gint nsym;
479 gint* numAtoms = NULL;
480 gint nAtoms = 0;
481 gint nC = getNumberOfBasisCenters(fileName, title);
482
483 if ((!fileName) || (strcmp(fileName,"") == 0))
484 {
485 Message("Sorry\n No file selected","Error",TRUE);
486 return FALSE;
487 }
488
489 t=g_malloc(taille*sizeof(gchar));
490 forb = FOpen(fileName, "rb");
491 if(forb == NULL)
492 {
493 gchar buffer[BSIZE];
494 sprintf(buffer,"Sorry, I can not open '%s' file\n",fileName);
495 Message(buffer,"Error",TRUE);
496 return FALSE;
497 }
498 ok = FALSE;
499 while(!feof(forb))
500 {
501 { char* e = fgets(t,taille,forb);}
502 pdest = strstr(t,title);
503 if(pdest != NULL)
504 {
505 ok = TRUE;
506 break;
507 }
508 }
509 if(!ok)
510 {
511 gchar buffer[BSIZE];
512 sprintf(buffer,"Sorry\nI can not read basis from '%s' file\n",fileName);
513 Message(buffer,"Error",TRUE);
514 return FALSE;
515 }
516
517 sym=g_malloc(10*sizeof(char));
518 /* printf("nC = %d\n",nC);*/
519 /* basis available for all centers */
520 if(nC==nCenters)
521 {
522 Ntype = nCenters;
523 for(i=0;i<nCenters;i++) GeomOrb[i].NumType= i;
524 }
525
526 if(forb !=NULL)
527 {
528 /* Debug("Ntype = %d\n",Ntype);*/
529 numAtoms = g_malloc(nCenters*sizeof(gint));
530 for(i=0;i<nCenters;i++) numAtoms[i] = -1;
531 nAtoms = 0;
532 Type = g_malloc(Ntype*sizeof(TYPE));
533 for(i=0;i<Ntype;i++)
534 {
535 Type[i].Ao = NULL;
536 Type[i].Norb=0;
537 }
538 sprintf(t," ");
539 while(!feof(forb))
540 {
541 if(strstr(t,"[")) break;
542 { char* e = fgets(t,taille,forb);}
543 /* Debug("tav = %s\n",t);*/
544 if(this_is_a_backspace(t) || strstr(t,"["))
545 {
546 g_free(sym);
547 g_free(t);
548 for(i=0;i<Ntype;i++)
549 if(Type[i].Ao == NULL)
550 {
551 gchar buffer[BSIZE];
552 sprintf(buffer,"Sorry, I can not read '%s' file, problem with basis set (1)\n",fileName);
553 Message(buffer,"Error",TRUE);
554 printf("AO pour i = %d\n", i);
555 if(numAtoms) g_free(numAtoms);
556 return FALSE;
557 }
558 resortAtoms(numAtoms);
559 if(numAtoms) g_free(numAtoms);
560 return TRUE;
561 }
562
563 if(!this_is_a_backspace(t))
564 {
565 /*Debug("tap = %s\n",t);*/
566 i=atoi(t)-1;
567 if(i>-1 && i<nCenters) numAtoms[nAtoms] = i;
568 nAtoms++;
569 /*Debug("i1 = %d \n",i);*/
570 if(i>-1)
571 {
572 sym = g_strdup( GeomOrb[i].Symb);
573 i = GeomOrb[i].NumType;
574 /* printf("numType = %d\n",i);*/
575 Type[i].Symb=g_strdup(sym);
576 Type[i].N=GetNelectrons(sym);
577 j=-1;
578 while(!feof(forb))
579 {
580 { char* e = fgets(t,taille,forb);}
581 /* Debug("t = %s\n",t);*/
582 if(this_is_a_backspace(t) || strstr(t,"["))
583 {
584 /*Debug("This is a backspace\n");*/
585 break;
586 }
587 /*
588 else
589 Debug("This is not a backspace\n");
590 */
591 j++;
592 Type[i].Norb=j+1;
593 /*
594 Debug("debut Alloc %d %d \n",i,j);
595 Debug("point %d \n",Type[i].Ao);
596 */
597 if(j == 0) Type[i].Ao=g_malloc(sizeof(AO));
598 else Type[i].Ao=g_realloc(Type[i].Ao,(j+1)*sizeof(AO));
599
600 /*
601 Debug("debut ReadOne i=%d j = %d \n",i,j);
602 Debug("debut t = %s \n",t);
603 */
604 ok= ReadOneBasis(i,j,t,&nsym);
605 /*Debug("nsym apres = %d\n",nsym);*/
606 if(nsym==2)
607 {
608 Type[i].Norb = j+2;
609 j++;
610 }
611 if(!ok)
612 {
613 j = j - nsym;
614 if(j==0) Type[i].Ao=g_malloc(sizeof(AO));
615 else Type[i].Ao=g_realloc(Type[i].Ao,(j+1)*sizeof(AO));
616 Type[i].Norb=j+1;
617 }
618
619 }
620 }
621 else
622 {
623 /*Debug("else = %s\n",t);*/
624 if(this_is_a_backspace(t) || strstr(t,"["))
625 {
626 g_free(sym);
627 g_free(t);
628 for(i=0;i<Ntype;i++)
629 if(Type[i].Ao == NULL)
630 {
631 gchar buffer[BSIZE];
632 sprintf(buffer,"Sorry, I can not read '%s' file, problem with basis set (2)\n",fileName);
633 Message(buffer,"Error",TRUE);
634 if(numAtoms) g_free(numAtoms);
635 return FALSE;
636 }
637 resortAtoms(numAtoms);
638 if(numAtoms) g_free(numAtoms);
639 return TRUE;
640 }
641 break;
642 }
643 }
644 }
645 }
646 for(i=0;i<Ntype;i++)
647 if(Type[i].Ao == NULL)
648 {
649 gchar buffer[BSIZE];
650 sprintf(buffer,"Sorry, I can not read '%s' file, problem with basis set (3)\n",fileName);
651 Message(buffer,"Error",TRUE);
652 if(numAtoms) g_free(numAtoms);
653 return FALSE;
654 }
655 resortAtoms(numAtoms);
656 if(numAtoms) g_free(numAtoms);
657 return TRUE;
658 }
659 /**********************************************/
DefineMoldenBasisType(gchar * fileName)660 gboolean DefineMoldenBasisType(gchar *fileName)
661 {
662 /* Debug("begin of DefineMoldenBasisType\n");*/
663 return DefineGabeditMoldenBasisType(fileName,"[GTO]");
664 }
665 /**********************************************/
DefineGabeditBasisType(gchar * fileName)666 gboolean DefineGabeditBasisType(gchar *fileName)
667 {
668 /* Debug("begin of DefineGabeditBasisType\n");*/
669 return DefineGabeditMoldenBasisType(fileName,"[Basis]");
670 }
671 /**********************************************/
PrintAllBasis()672 void PrintAllBasis()
673 {
674 gint k,j,n;
675 gint l;
676 char *XYZ[]={"x","y","z"};
677
678
679 for(k=0;k<NAOrb;k++)
680 {
681
682 for(n=0;n<AOrb[k].numberOfFunctions;n++)
683 {
684
685 l=0;
686 for(j=0;j<3;j++)
687 l += AOrb[k].Gtf[n].l[j];
688
689 Debug("%c",GetSymmetry(l));
690 for(j=0;j<3;j++)
691 {
692 switch(AOrb[k].Gtf[n].l[j])
693 {
694 case 0:break;
695 case 1:Debug("%s",XYZ[j]);break;
696 default :Debug("%s%d",XYZ[j],AOrb[k].Gtf[n].l[j]);
697 }
698 }
699 Debug("\t%9.6f %9.6f \tC=%9.6f %9.6f %9.6f",AOrb[k].Gtf[n].Ex,AOrb[k].Gtf[n].Coef, AOrb[k].Gtf[n].C[0],AOrb[k].Gtf[n].C[1],AOrb[k].Gtf[n].C[2]);
700 Debug("\n");
701 }
702 Debug("\n");
703 }
704 }
705 /**********************************************/
NormaliseAllBasis()706 void NormaliseAllBasis()
707 {
708 gint k,n;
709
710 for(k=0;k<NAOrb;k++)
711 for(n=0;n<AOrb[k].numberOfFunctions;n++)
712 normaliseRadialGTF(&AOrb[k].Gtf[n]);
713
714 for(k=0;k<NAOrb;k++)
715 normaliseCGTF(&AOrb[k]);
716 }
717 /**********************************************/
NormaliseAllNoRadBasis()718 void NormaliseAllNoRadBasis()
719 {
720 gint k;
721
722 for(k=0;k<NAOrb;k++)
723 normaliseCGTF(&AOrb[k]);
724 }
725 /**********************************************/
726 /**********************************************/
DefineAtomicNumOrb()727 void DefineAtomicNumOrb()
728 {
729 gint i;
730 gint j;
731 for(i=0;i<nCenters;i++)
732 {
733 /* Debug("i= %d \n",i);*/
734 GeomOrb[i].NAOrb = 0;
735 GeomOrb[i].NAlphaOrb = 0;
736 GeomOrb[i].NBetaOrb = 0;
737 GeomOrb[i].NumOrb = NULL;
738 GeomOrb[i].CoefAlphaOrbitals = NULL;
739 GeomOrb[i].OccAlphaOrbitals = NULL;
740 GeomOrb[i].EnerAlphaOrbitals = NULL;
741 GeomOrb[i].SymAlphaOrbitals = NULL;
742 GeomOrb[i].CoefBetaOrbitals = NULL;
743 GeomOrb[i].EnerBetaOrbitals = NULL;
744 GeomOrb[i].OccBetaOrbitals = NULL;
745 GeomOrb[i].SymBetaOrbitals = NULL;
746 }
747 /* Debug("End Geom init \n");*/
748
749 for(j=0;j<NAOrb;j++)
750 {
751 if(AOrb) i = AOrb[j].NumCenter;
752 else if(SAOrb) i = SAOrb[j].NumCenter;
753 GeomOrb[i].NAOrb++;
754 if(!GeomOrb[i].NumOrb)
755 GeomOrb[i].NumOrb = g_malloc(GeomOrb[i].NAOrb*sizeof(gint));
756 else
757 GeomOrb[i].NumOrb = g_realloc(GeomOrb[i].NumOrb,GeomOrb[i].NAOrb*sizeof(gint));
758 GeomOrb[i].NumOrb[GeomOrb[i].NAOrb-1] = j;
759 }
760 }
761 /**********************************************/
DefineCartBasis()762 void DefineCartBasis()
763 {
764 gint i,j,k,n;
765 gint l1,l2,l3;
766 gint L;
767 gint *l[3]={NULL,NULL,NULL};
768 gint m;
769
770 NAOrb = 0;
771 for(i=0;i<nCenters;i++)
772 {
773 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
774 {
775 L=Type[GeomOrb[i].NumType].Ao[j].L;
776 NAOrb += (L+1)*(L+2)/2;
777 }
778 }
779
780 AOrb = g_malloc(NAOrb*sizeof(CGTF));
781 if(SAOrb) g_free(SAOrb);
782 SAOrb = NULL;
783
784 k=-1;
785 for(i=0;i<nCenters;i++)
786 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
787 {
788 L = Type[GeomOrb[i].NumType].Ao[j].L;
789
790 for(m=0;m<3;m++)
791 {
792 if(l[m])
793 g_free(l[m]);
794 l[m] = g_malloc((L+1)*(L+2)/2*sizeof(gint));
795 }
796 switch(L)
797 {
798 case 1 :
799 m=0;
800 l[0][m] = 1;l[1][m] = 0;l[2][m] = 0; /* X */
801 m++;
802 l[0][m] = 0;l[1][m] = 1;l[2][m] = 0; /* Y */
803 m++;
804 l[0][m] = 0;l[1][m] = 0;l[2][m] = 1; /* Z */
805 break;
806 case 2 :
807 m=0;
808 l[0][m] = 2;l[1][m] = 0;l[2][m] = 0; /* XX */
809 m++;
810 l[0][m] = 0;l[1][m] = 2;l[2][m] = 0; /* YY */
811 m++;
812 l[0][m] = 0;l[1][m] = 0;l[2][m] = 2; /* ZZ */
813 m++;
814 l[0][m] = 1;l[1][m] = 1;l[2][m] = 0; /* XY */
815 m++;
816 l[0][m] = 1;l[1][m] = 0;l[2][m] = 1; /* XZ */
817 m++;
818 l[0][m] = 0;l[1][m] = 1;l[2][m] = 1; /* YZ */
819 break;
820 case 3 :
821 m=0;
822 l[0][m] = 3;l[1][m] = 0;l[2][m] = 0; /* XXX */
823 m++;
824 l[0][m] = 0;l[1][m] = 3;l[2][m] = 0; /* YYY */
825 m++;
826 l[0][m] = 0;l[1][m] = 0;l[2][m] = 3; /* ZZZ */
827 m++;
828 l[0][m] = 1;l[1][m] = 2;l[2][m] = 0; /* XYY */
829 m++;
830 l[0][m] = 2;l[1][m] = 1;l[2][m] = 0; /* XXY */
831 m++;
832 l[0][m] = 2;l[1][m] = 0;l[2][m] = 1; /* XXZ */
833 m++;
834 l[0][m] = 1;l[1][m] = 0;l[2][m] = 2; /* XZZ */
835 m++;
836 l[0][m] = 0;l[1][m] = 1;l[2][m] = 2; /* YZZ */
837 m++;
838 l[0][m] = 0;l[1][m] = 2;l[2][m] = 1; /* YYZ */
839 m++;
840 l[0][m] = 1;l[1][m] = 1;l[2][m] = 1; /* XYZ */
841 break;
842 default :
843 m=0;
844 for(l3=Type[GeomOrb[i].NumType].Ao[j].L;l3>=0;l3--)
845 for(l2=Type[GeomOrb[i].NumType].Ao[j].L-l3;l2>=0;l2--)
846 {
847 l1 = Type[GeomOrb[i].NumType].Ao[j].L-l2-l3;
848 l[0][m] = l1;
849 l[1][m] = l2;
850 l[2][m] = l3;
851 m++;
852 }
853 }
854 for(m=0;m<(L+1)*(L+2)/2;m++)
855 {
856 l1 = l[0][m];
857 l2 = l[1][m];
858 l3 = l[2][m];
859 k++;
860 AOrb[k].numberOfFunctions=Type[GeomOrb[i].NumType].Ao[j].N;
861 AOrb[k].NumCenter = i;
862 AOrb[k].Gtf =g_malloc(AOrb[k].numberOfFunctions*sizeof(GTF));
863 for(n=0;n<AOrb[k].numberOfFunctions;n++)
864 {
865 AOrb[k].Gtf[n].Ex = Type[GeomOrb[i].NumType].Ao[j].Ex[n];
866 AOrb[k].Gtf[n].Coef = Type[GeomOrb[i].NumType].Ao[j].Coef[n];
867 AOrb[k].Gtf[n].C[0] = GeomOrb[i].C[0];
868 AOrb[k].Gtf[n].C[1] = GeomOrb[i].C[1];
869 AOrb[k].Gtf[n].C[2] = GeomOrb[i].C[2];
870 AOrb[k].Gtf[n].l[0] = l1;
871 AOrb[k].Gtf[n].l[1] = l2;
872 AOrb[k].Gtf[n].l[2] = l3;
873 }
874
875 }
876 }
877
878 NOrb = NAOrb;
879 DefineAtomicNumOrb();
880 /* DefineNorb();*/
881 }
882 /**********************************************/
DefineSphericalBasis()883 void DefineSphericalBasis()
884 {
885 gint i,j,k;
886 gint c;
887 gint kl;
888 gint L,M;
889 CGTF *temp;
890 Zlm Stemp;
891 gint N,Nc,n;
892 gint inc;
893 gint klbeg;
894 gint klend;
895 gint klinc;
896
897
898 NOrb = 0;
899 for(i=0;i<nCenters;i++)
900 {
901 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
902 {
903 L=Type[GeomOrb[i].NumType].Ao[j].L;
904 NOrb += 2*L+1;
905 }
906 }
907
908 temp = g_malloc(NOrb*sizeof(CGTF));
909
910 k=-1;
911 for(i=0;i<nCenters;i++)
912 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
913 {
914 L =Type[GeomOrb[i].NumType].Ao[j].L;
915 /*Debug("L=%d \n",L);*/
916
917 /*Debug("L =%d \n",L);*/
918 if(L==1)
919 {
920 klbeg = L;
921 klend = 0;
922 klinc = -1;
923 }
924 else
925 {
926 klbeg = 0;
927 klend = L;
928 klinc = +1;
929 }
930 for(kl = klbeg;(klbeg == 0 && kl<=klend) || (klbeg == L && kl>=klend);kl +=klinc)
931 {
932 if(kl!=0)
933 inc = 2*kl;
934 else
935 inc = 1;
936 for(M=kl;M>=-kl;M -=inc)
937 {
938 /*Debug("L =%d kl=%d M=%d \n",L,kl,M);*/
939 k++;
940 Stemp = getZlm(L,M);
941
942 temp[k].numberOfFunctions=Stemp.numberOfCoefficients*Type[GeomOrb[i].NumType].Ao[j].N;
943 temp[k].NumCenter=i;
944 /*Debug("M=%d N=%d\n",M,temp[k].N);*/
945 temp[k].Gtf =g_malloc(temp[k].numberOfFunctions*sizeof(GTF));
946 Nc=-1;
947 for(N=0;N<Type[GeomOrb[i].NumType].Ao[j].N;N++)
948 for(n=0;n<Stemp.numberOfCoefficients;n++)
949 {
950 Nc++;
951
952 temp[k].Gtf[Nc].Ex = Type[GeomOrb[i].NumType].Ao[j].Ex[N];
953 temp[k].Gtf[Nc].Coef = Type[GeomOrb[i].NumType].Ao[j].Coef[N]*Stemp.lxyz[n].Coef;
954 for(c=0;c<3;c++)
955 {
956 temp[k].Gtf[Nc].C[c] = GeomOrb[i].C[c];
957 temp[k].Gtf[Nc].l[c] = Stemp.lxyz[n].l[c];
958 }
959 }
960 if(L==0)
961 break;
962 }
963 if(L==0)
964 break;
965 }
966 }
967 for(i=0;i<NAOrb;i++)
968 g_free(AOrb[i].Gtf);
969 g_free(AOrb);
970 NAOrb = NOrb;
971 AOrb = temp;
972 if(SAOrb) g_free(SAOrb);
973 SAOrb = NULL;
974 DefineAtomicNumOrb();
975 /* DefineNorb();*/
976 }
977 /***************************************************************/
getlTable(gint L,gint * nCoefs,gdouble ** coefs,gint ** l[3])978 static void getlTable(gint L, gint* nCoefs, gdouble** coefs, gint** l[3])
979 {
980 gint m;
981 gint l1,l2,l3;
982 gint n;
983 /* Spehrical D,F, ...*/
984 if(L<-1)
985 {
986 gint klbeg=0, klend=-L, klinc=1;
987 gint inc;
988 gint kl;
989 gint M;
990 gint c;
991 Zlm Stemp;
992 gint m = -1;
993 L = abs(L);
994 if(L==1)
995 {
996 klbeg = L;
997 klend = 0;
998 klinc = -1;
999 }
1000 else
1001 {
1002 klbeg = 0;
1003 klend = L;
1004 klinc = +1;
1005 }
1006 for(kl = klbeg;(klbeg == 0 && kl<=klend) || (klbeg == L && kl>=klend);kl +=klinc)
1007 {
1008 if(kl!=0) inc = 2*kl;
1009 else inc = 1;
1010 for(M=kl;M>=-kl;M -=inc)
1011 {
1012 if(L==1) m = M+abs(L);
1013 else m++;
1014 Stemp = getZlm(L,M);
1015 nCoefs[m] = Stemp.numberOfCoefficients;
1016 for(n=0;n<Stemp.numberOfCoefficients;n++)
1017 {
1018 coefs[m][n] = Stemp.lxyz[n].Coef;
1019 for(c=0;c<3;c++) l[c][m][n] = Stemp.lxyz[n].l[c];
1020 }
1021 if(L==0) break;
1022 }
1023 if(L==0) break;
1024 }
1025 return;
1026 }
1027 /* Cartesian S,P,D,F,..*/
1028 L = abs(L); /* for P */
1029 for(m=0;m<(L+1)*(L+2)/2;m++)
1030 {
1031 nCoefs[m] = 1;
1032 coefs[m][0] = 1.0;
1033 }
1034 switch(L)
1035 {
1036 case 0 :
1037 m=0;
1038 l[0][m][0] = 0;l[1][m][0] = 0;l[2][m][0] = 0;
1039 break;
1040 case 1 :
1041 m=0;
1042 l[0][m][0] = 1;l[1][m][0] = 0;l[2][m][0] = 0; /* X */
1043 m++;
1044 l[0][m][0] = 0;l[1][m][0] = 1;l[2][m][0] = 0; /* Y */
1045 m++;
1046 l[0][m][0] = 0;l[1][m][0] = 0;l[2][m][0] = 1; /* Z */
1047 break;
1048 case 2 :
1049 m=0;
1050 l[0][m][0] = 2;l[1][m][0] = 0;l[2][m][0] = 0; /* XX */
1051 m++;
1052 l[0][m][0] = 0;l[1][m][0] = 2;l[2][m][0] = 0; /* YY */
1053 m++;
1054 l[0][m][0] = 0;l[1][m][0] = 0;l[2][m][0] = 2; /* ZZ */
1055 m++;
1056 l[0][m][0] = 1;l[1][m][0] = 1;l[2][m][0] = 0; /* XY */
1057 m++;
1058 l[0][m][0] = 1;l[1][m][0] = 0;l[2][m][0] = 1; /* XZ */
1059 m++;
1060 l[0][m][0] = 0;l[1][m][0] = 1;l[2][m][0] = 1; /* YZ */
1061 break;
1062 case 3 :
1063 m=0;
1064 l[0][m][0] = 3;l[1][m][0] = 0;l[2][m][0] = 0; /* XXX */
1065 m++;
1066 l[0][m][0] = 0;l[1][m][0] = 3;l[2][m][0] = 0; /* YYY */
1067 m++;
1068 l[0][m][0] = 0;l[1][m][0] = 0;l[2][m][0] = 3; /* ZZZ */
1069 m++;
1070 l[0][m][0] = 1;l[1][m][0] = 2;l[2][m][0] = 0; /* XYY */
1071 m++;
1072 l[0][m][0] = 2;l[1][m][0] = 1;l[2][m][0] = 0; /* XXY */
1073 m++;
1074 l[0][m][0] = 2;l[1][m][0] = 0;l[2][m][0] = 1; /* XXZ */
1075 m++;
1076 l[0][m][0] = 1;l[1][m][0] = 0;l[2][m][0] = 2; /* XZZ */
1077 m++;
1078 l[0][m][0] = 0;l[1][m][0] = 1;l[2][m][0] = 2; /* YZZ */
1079 m++;
1080 l[0][m][0] = 0;l[1][m][0] = 2;l[2][m][0] = 1; /* YYZ */
1081 m++;
1082 l[0][m][0] = 1;l[1][m][0] = 1;l[2][m][0] = 1; /* XYZ */
1083 break;
1084 default :
1085 m=0;
1086 for(l3=abs(L);l3>=0;l3--)
1087 for(l2=abs(L)-l3;l2>=0;l2--)
1088 {
1089 l1 = abs(L)-l2-l3;
1090 l[0][m][0] = l1;
1091 l[1][m][0] = l2;
1092 l[2][m][0] = l3;
1093 m++;
1094 }
1095 }
1096 }
1097 /********************************************************************************/
readBasisFromGaussianFChk(gchar * fileName)1098 gboolean readBasisFromGaussianFChk(gchar *fileName)
1099 {
1100 FILE *file;
1101 gint n;
1102 gint nS;
1103 gint c;
1104 gint nShells = 0;
1105 gint nPrimitives = 0;
1106 gint lMax = 0;
1107 gint contMax = 0;
1108 gint* shellTypes = NULL;
1109 gint* nPrimitivesByShell = NULL;
1110 gint* numAtoms = NULL;
1111 gdouble* primitiveExponents = NULL;
1112 gdouble* contractionsCoefs = NULL;
1113 gdouble* contractionsCoefsSP = NULL;
1114 gdouble* coordinatesForShells = NULL;
1115 gint** l[3] = {NULL,NULL,NULL};
1116 gdouble** coefs = NULL;
1117 gint* nCoefs = NULL;
1118 CGTF *temp = NULL;
1119 gint kOrb, kPrimitive;
1120 gint m;
1121 gint nBasis;
1122 gboolean sp = FALSE;
1123 gint llMax = 0;
1124
1125 file = FOpen(fileName, "rb");
1126 if(file ==NULL)
1127 {
1128 Message(_("Sorry\nI can not open this file"),_("Error"),TRUE);
1129 return FALSE;
1130 }
1131
1132 nBasis = get_one_int_from_fchk_gaussian_file(file,"Number of basis functions ");
1133 if(nBasis<1)
1134 {
1135 Message(_("Sorry\nI can not read the number of basis functions"),_("Error"),TRUE);
1136 fclose(file);
1137 return FALSE;
1138 }
1139
1140 nShells = get_one_int_from_fchk_gaussian_file(file,"Number of contracted shells ");
1141 if(nShells<1)
1142 {
1143 Message(_("Sorry\nI can not the number of contracted shells"),_("Error"),TRUE);
1144 fclose(file);
1145 return FALSE;
1146 }
1147 nPrimitives = get_one_int_from_fchk_gaussian_file(file,"Number of primitive shells ");
1148 if(nPrimitives<1)
1149 {
1150 Message(_("Sorry\nI can not the number of primitive shells"),_("Error"),TRUE);
1151 fclose(file);
1152 return FALSE;
1153 }
1154 rewind(file);
1155 lMax = get_one_int_from_fchk_gaussian_file(file,"Highest angular momentum ");
1156 if(lMax<0)
1157 {
1158 Message(_("Sorry\nI can not the value of the highest angular momentum"),_("Error"),TRUE);
1159 fclose(file);
1160 return FALSE;
1161 }
1162 rewind(file);
1163 contMax = get_one_int_from_fchk_gaussian_file(file,"Largest degree of contraction ");
1164 if(contMax<1)
1165 {
1166 Message(_("Sorry\nI can not the value of the largest degree of contraction"),_("Error"),TRUE);
1167 fclose(file);
1168 return FALSE;
1169 }
1170 rewind(file);
1171 shellTypes = get_array_int_from_fchk_gaussian_file(file, "Shell types ", &n);
1172 if(!shellTypes || n!=nShells)
1173 {
1174 Message(_("Sorry\nI can not read the shell types"),_("Error"),TRUE);
1175 if(shellTypes) g_free(shellTypes);
1176 fclose(file);
1177 return FALSE;
1178 }
1179 for(nS = 0;nS<nShells;nS++) if( shellTypes[nS]==-1) { sp = TRUE; break;}
1180
1181 NOrb = 0;
1182 for(nS=0;nS<nShells;nS++)
1183 {
1184 if(shellTypes[nS]<-1) NOrb += 2*abs(shellTypes[nS])+1; /* Sperical D, F, G, ...*/
1185 else if(shellTypes[nS]==-1) NOrb += 4; /* This a SP.*/
1186 else NOrb += (shellTypes[nS]+1)*(shellTypes[nS]+2)/2; /* Cartezian S,P,D,F,G,..*/
1187 }
1188 if(NOrb != nBasis)
1189 {
1190 Message(_("Sorry\nThe number of basis function in fch file is not equal to that computed by Gabedit!"),_("Error"),TRUE);
1191 if(shellTypes) g_free(shellTypes);
1192 fclose(file);
1193 return FALSE;
1194 }
1195 nPrimitivesByShell = get_array_int_from_fchk_gaussian_file(file, "Number of primitives per shell ", &n);
1196 if(!nPrimitivesByShell || n!=nShells)
1197 {
1198 Message(_("Sorry\nI can not read the number of primitives per shell"),_("Error"),TRUE);
1199 if(nPrimitivesByShell) g_free(nPrimitivesByShell);
1200 fclose(file);
1201 return FALSE;
1202 }
1203 numAtoms = get_array_int_from_fchk_gaussian_file(file, "Shell to atom map ", &n);
1204 if(!numAtoms || n!=nShells)
1205 {
1206 Message(_("Sorry\nI can not read the atoms number for shell"),_("Error"),TRUE);
1207 if(numAtoms) g_free(numAtoms);
1208 fclose(file);
1209 return FALSE;
1210 }
1211 primitiveExponents = get_array_real_from_fchk_gaussian_file(file, "Primitive exponents ", &n);
1212 if(!primitiveExponents || n != nPrimitives)
1213 {
1214 Message(_("Sorry\nI can not read the primitive exponents "),_("Error"),TRUE);
1215 if(primitiveExponents) g_free(primitiveExponents);
1216 fclose(file);
1217 return FALSE;
1218 }
1219 contractionsCoefs = get_array_real_from_fchk_gaussian_file(file, "Contraction coefficients ", &n);
1220 if(!contractionsCoefs || n != nPrimitives)
1221 {
1222 Message(_("Sorry\nI can not read the contraction coefficients "),_("Error"),TRUE);
1223 if(contractionsCoefs) g_free(contractionsCoefs);
1224 fclose(file);
1225 return FALSE;
1226 }
1227 if(sp)
1228 {
1229 contractionsCoefsSP = get_array_real_from_fchk_gaussian_file(file, "P(S=P) Contraction coefficients ", &n);
1230 if(!contractionsCoefsSP || n != nPrimitives)
1231 {
1232 Message(_("Sorry\nI can not read the P(S=P) contraction coefficients "),_("Error"),TRUE);
1233 if(contractionsCoefsSP) g_free(contractionsCoefsSP);
1234 fclose(file);
1235 return FALSE;
1236 }
1237 }
1238 coordinatesForShells = get_array_real_from_fchk_gaussian_file(file, "Coordinates of each shell ", &n);
1239 if(!contractionsCoefs || n != nShells*3)
1240 {
1241 Message(_("Sorry\nI can not read the coordinates of each shell "),_("Error"),TRUE);
1242 if(coordinatesForShells) g_free(coordinatesForShells);
1243 fclose(file);
1244 return FALSE;
1245 }
1246 fclose(file);
1247 /* printf("close file\n");*/
1248
1249 llMax = (lMax+1)*(lMax+2)/2;
1250 nCoefs = g_malloc(llMax*sizeof(gint));
1251 coefs = g_malloc(llMax*sizeof(gdouble*));
1252 for(m=0;m<llMax;m++) coefs[m] = g_malloc(llMax*sizeof(gdouble));
1253
1254 for(c=0;c<3;c++)
1255 {
1256 l[c] = g_malloc(llMax*sizeof(gint*));
1257 for(m=0;m<llMax;m++) l[c][m] = g_malloc(llMax*sizeof(gint));
1258 }
1259
1260 temp = g_malloc(NOrb*sizeof(CGTF));
1261 kOrb = 0;
1262 kPrimitive = 0;
1263 for(nS = 0;nS<nShells; nS++)
1264 {
1265 gint nM = 0;
1266 /* printf("begin primitive nS = %d\n",nS);*/
1267 if(shellTypes[nS]<-1) nM = 2*abs(shellTypes[nS])+1; /* Sperical D, F, G, ...*/
1268 else if(shellTypes[nS]==-1) nM = 1; /* This a SP. Make S befor */
1269 else nM = (shellTypes[nS]+1)*(shellTypes[nS]+2)/2;
1270
1271 /* printf("nM = %d\n",nM);*/
1272 if(shellTypes[nS]==-1) getlTable(0, nCoefs, coefs, l); /* This a SP. Make S befor */
1273 else getlTable(shellTypes[nS], nCoefs, coefs, l);
1274 /* printf("end getlTable\n");*/
1275 for(m=0;m<nM;m++)
1276 {
1277 gint ip,j,n;
1278 temp[kOrb].numberOfFunctions=nCoefs[m]*nPrimitivesByShell[nS];
1279 temp[kOrb].NumCenter=numAtoms[nS]-1;
1280 /* printf("m = %d nCoef = %d nPrim = %d\n",m,nCoefs[m],nPrimitivesByShell[nS]);*/
1281 /*Debug("M=%d N=%d\n",M,temp[k].N);*/
1282 temp[kOrb].Gtf =g_malloc(temp[kOrb].numberOfFunctions*sizeof(GTF));
1283 j = -1;
1284 for(ip=0;ip<nPrimitivesByShell[nS];ip++)
1285 for(n=0;n<nCoefs[m];n++)
1286 {
1287 j++;
1288 temp[kOrb].Gtf[j].Ex = primitiveExponents[kPrimitive+ip];
1289 temp[kOrb].Gtf[j].Coef = contractionsCoefs[kPrimitive+ip]*coefs[m][n];
1290 for(c=0;c<3;c++)
1291 {
1292 temp[kOrb].Gtf[j].C[c] = coordinatesForShells[c+nS*3];
1293 temp[kOrb].Gtf[j].l[c] = l[c][m][n];
1294 }
1295 }
1296 kOrb++;
1297 }
1298 if(shellTypes[nS]==-1) /* This a SP. Now make P*/
1299 {
1300 getlTable(-1, nCoefs, coefs, l);
1301 nM = 3;
1302 for(m=0;m<nM;m++)
1303 {
1304 gint ip,j,n;
1305 /* printf("P : m = %d nCoef = %d nPrim = %d\n",m,nCoefs[m],nPrimitivesByShell[nS]);*/
1306 temp[kOrb].numberOfFunctions=nCoefs[m]*nPrimitivesByShell[nS];
1307 temp[kOrb].NumCenter=numAtoms[nS]-1;
1308 temp[kOrb].Gtf =g_malloc(temp[kOrb].numberOfFunctions*sizeof(GTF));
1309 j = -1;
1310 for(ip=0;ip<nPrimitivesByShell[nS];ip++)
1311 for(n=0;n<nCoefs[m];n++)
1312 {
1313 j++;
1314 temp[kOrb].Gtf[j].Ex = primitiveExponents[kPrimitive+ip];
1315 temp[kOrb].Gtf[j].Coef = contractionsCoefsSP[kPrimitive+ip]*coefs[m][n];
1316 for(c=0;c<3;c++)
1317 {
1318 temp[kOrb].Gtf[j].C[c] = coordinatesForShells[c+nS*3];
1319 temp[kOrb].Gtf[j].l[c] = l[c][m][n];
1320 }
1321 }
1322 kOrb++;
1323 }
1324 }
1325 /* printf("end primitive nS = %d\n",nS);*/
1326 kPrimitive += nPrimitivesByShell[nS];
1327 }
1328 if(nCoefs) g_free(nCoefs);
1329 if(coefs)
1330 {
1331 for(m=0;m<llMax;m++) if(coefs[m]) g_free(coefs[m]);
1332 g_free(coefs);
1333 }
1334 for(c=0;c<3;c++)
1335 if(l[c])
1336 {
1337 for(m=0;m<llMax;m++) if(l[c][m]) g_free(l[c][m]);
1338 g_free(l[c]);
1339 }
1340 if(shellTypes) g_free(shellTypes);
1341 if(nPrimitivesByShell) g_free(nPrimitivesByShell);
1342 if(numAtoms) g_free(numAtoms);
1343 if(primitiveExponents) g_free(primitiveExponents);
1344 if(contractionsCoefs) g_free(contractionsCoefs);
1345 if(coordinatesForShells) g_free(coordinatesForShells);
1346 for(kOrb=0;kOrb<NAOrb;kOrb++) g_free(AOrb[kOrb].Gtf);
1347 if(AOrb) g_free(AOrb);
1348 NAOrb = NOrb;
1349 AOrb = temp;
1350 if(SAOrb) g_free(SAOrb);
1351 SAOrb = NULL;
1352 DefineAtomicNumOrb();
1353 /* printf("End all read Basis\n");*/
1354 return TRUE;
1355 }
1356