1 /* OrbitalsNBO.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 "../Utils/AtomsProp.h"
23 #include "../Utils/UtilsInterface.h"
24 #include "../Utils/Utils.h"
25 #include "../Utils/Zlm.h"
26 #include "../Utils/Constants.h"
27 #include "../Geometry/GeomGlobal.h"
28 #include "GeomDraw.h"
29 #include "GLArea.h"
30 #include "UtilsOrb.h"
31 #include "Basis.h"
32 #include "GeomOrbXYZ.h"
33 #include "AtomicOrbitals.h"
34 #include "StatusOrb.h"
35 #include "Basis.h"
36 #include "Orbitals.h"
37 #include "GeomOrbXYZ.h"
38 #include "BondsOrb.h"
39 #include "LabelsGL.h"
40
41 /********************************************************************************/
42 typedef enum
43 {
44 GABEDIT_ORBLOCALTYPE_BOYS=0,
45 GABEDIT_ORBLOCALTYPE_EDMISTON,
46 GABEDIT_ORBLOCALTYPE_PIPEK,
47 GABEDIT_ORBLOCALTYPE_UNKNOWN
48 } GabEditOrbLocalType;
49
50 typedef enum
51 {
52 GABEDIT_ORBTYPE_ALPHA = 0,
53 GABEDIT_ORBTYPE_BETA,
54 GABEDIT_ORBTYPE_MOLECULAR,
55 GABEDIT_ORBTYPE_MCSCF,
56 GABEDIT_ORBTYPE_EIGENVECTORS,
57 GABEDIT_ORBTYPE_BOYS_ALPHA,
58 GABEDIT_ORBTYPE_BOYS_BETA,
59 GABEDIT_ORBTYPE_BOYS,
60 GABEDIT_ORBTYPE_EDMISTON_ALPHA,
61 GABEDIT_ORBTYPE_EDMISTON_BETA,
62 GABEDIT_ORBTYPE_EDMISTON,
63 GABEDIT_ORBTYPE_PIPEK_ALPHA,
64 GABEDIT_ORBTYPE_PIPEK_BETA,
65 GABEDIT_ORBTYPE_PIPEK,
66 } GabEditOrbType;
67 /********************************************************************************/
68 static gboolean sphericalBasis = TRUE;
69 /********************************************************************************/
get_charges_from_nbo_output_file(FILE * file,gint nAtoms)70 static void get_charges_from_nbo_output_file(FILE* file,gint nAtoms)
71 {
72 /* charges not available in a nbo file. I set it to 0.0 */
73 gint i;
74 for(i=0;i<nAtoms;i++)
75 GeomOrb[i].partialCharge = 0.0;
76 }
77 /********************************************************************************/
goToLine(FILE * file,char * nextString)78 static gboolean goToLine(FILE* file,char* nextString)
79 {
80 static char t[BSIZE];
81 while(!feof(file))
82 {
83 if(!fgets(t,BSIZE,file))break;
84 if (strstr(t,nextString)) return TRUE;
85 }
86 return FALSE;
87 }
88 /********************************************************************************/
setTitle(FILE * file)89 static void setTitle(FILE* file)
90 {
91 gchar t[BSIZE];
92
93 if(!file)
94 {
95 set_label_title("",0,0);
96 return;
97 }
98 rewind(file);
99 if(fgets(t,BSIZE,file))
100 if(fgets(t,BSIZE,file))
101 set_label_title(t,0,0);
102 }
103 /********************************************************************************/
read_geomorb_nbo_file_geom(gchar * fileName)104 static gboolean read_geomorb_nbo_file_geom(gchar *fileName)
105 {
106 gchar *tmp = NULL;
107 FILE *file;
108 gint i;
109 gint k;
110 gint j=0;
111 gint uni=1;
112 gint z = 0;
113 static gchar t[BSIZE];
114 gint nAtoms, nShell, nExp;
115
116 if ((!fileName) || (strcmp(fileName,"") == 0))
117 {
118 Message(_("Sorry\n No file selected"),_("Error"),TRUE);
119 return FALSE;
120 }
121
122 file = FOpen(fileName, "rb");
123 if(file == NULL)
124 {
125 Message(_("Sorry\nI cannot open this file"),_("Error"),TRUE);
126 return FALSE;
127 }
128 if(!goToLine(file,"--------")) return FALSE;
129 if(!fgets(t,BSIZE,file)) return FALSE;
130 sscanf(t,"%d %d %d",&nAtoms,&nShell,&nExp);
131 if(!goToLine(file,"--------")) return FALSE;
132 if(nAtoms<1) return FALSE;
133 init_dipole();
134 GeomOrb=g_malloc(nAtoms*sizeof(TypeGeomOrb));
135 uni = 1;
136
137 tmp = get_name_file(fileName);
138 set_status_label_info(_("File name"),tmp);
139 g_free(tmp);
140 set_status_label_info(_("File type"),"NBO");
141 set_status_label_info(_("Geometry"),_("Reading"));
142
143 j = 0;
144 for(k=0;k<nAtoms;k++)
145 {
146 if(!fgets(t,BSIZE,file)) break;
147 sscanf(t,"%d %lf %lf %lf",&z, &GeomOrb[j].C[0], &GeomOrb[j].C[1], &GeomOrb[j].C[2]);
148 if(z<=0) GeomOrb[j].Symb=g_strdup("X");
149 else GeomOrb[j].Symb=get_symbol_using_z(z);
150 if(uni==1) for(i=0;i<3;i++) GeomOrb[j].C[i] *= ANG_TO_BOHR;
151 GeomOrb[j].Prop = prop_atom_get(GeomOrb[j].Symb);
152 GeomOrb[j].nuclearCharge = get_atomic_number_from_symbol(GeomOrb[j].Symb);
153 GeomOrb[j].variable = TRUE;
154 j++;
155
156 }
157 nCenters = 0;
158 if(k==nAtoms) nCenters = nAtoms;
159 if(nCenters !=0) get_charges_from_nbo_output_file(file,nCenters);
160 fclose(file);
161 if(nCenters == 0 )
162 {
163 g_free(GeomOrb);
164 sprintf(t,_("Sorry, I can not read this format from '%s' file\n"),fileName);
165 Message(t,_("Error"),TRUE);
166 set_status_label_info(_("File name"),_("Nothing"));
167 set_status_label_info(_("File type"),_("Nothing"));
168 set_status_label_info(_("Mol. Orb."),_("Nothing"));
169 RebuildGeom = TRUE;
170 return FALSE;
171 }
172 else
173 {
174 /* DefineType();*/
175 gint i;
176 if(Type)
177 for(i=0;i<Ntype;i++)
178 {
179 if(Type[i].Ao != NULL)
180 {
181 for(j=0;j<Type[i].Norb;j++)
182 {
183 if(Type[i].Ao[j].Ex != NULL) g_free(Type[i].Ao[j].Ex);
184 if(Type[i].Ao[j].Coef != NULL) g_free(Type[i].Ao[j].Coef);
185 }
186 g_free(Type[i].Ao);
187 }
188 }
189 if(Type) g_free(Type);
190 Ntype =nCenters;
191 Type = g_malloc(Ntype*sizeof(TYPE));
192 for(i=0;i<nCenters;i++)
193 {
194 GeomOrb[i].NumType = i;
195 Type[i].Symb=g_strdup(GeomOrb[i].Symb);
196 Type[i].N=GetNelectrons(GeomOrb[i].Symb);
197 Type[i].Norb = 0;
198 Type[i].Ao = NULL;
199 }
200 buildBondsOrb();
201 RebuildGeom = TRUE;
202 return TRUE;
203 }
204 return TRUE;
205 }
206 /**********************************************/
DefineNBOSphericalBasis()207 static void DefineNBOSphericalBasis()
208 {
209 gint i,j,k;
210 gint c;
211 gint kl;
212 gint L,M;
213 CGTF *temp;
214 Zlm Stemp;
215 gint N,Nc,n;
216 gint inc;
217 gint klbeg;
218 gint klend;
219 gint klinc;
220
221
222 NOrb = 0;
223 for(i=0;i<nCenters;i++)
224 {
225 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
226 {
227 L=Type[GeomOrb[i].NumType].Ao[j].L;
228 NOrb += 2*L+1;
229 }
230 }
231
232 temp = g_malloc(NOrb*sizeof(CGTF));
233
234 k=-1;
235 for(i=0;i<nCenters;i++)
236 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
237 {
238 L =Type[GeomOrb[i].NumType].Ao[j].L;
239 /* Debug("L=%d \n",L);*/
240
241 /*Debug("L =%d \n",L);*/
242 if(L==1)
243 {
244 klbeg = 0;
245 klend = L;
246 klinc = +1;
247 }
248 else
249 {
250 klbeg = 0;
251 klend = L;
252 klinc = +1;
253 }
254 for(kl = klbeg;(klbeg == 0 && kl<=klend) || (klbeg == L && kl>=klend);kl +=klinc)
255 {
256 if(kl!=0) inc = 2*kl;
257 else inc = 1;
258 for(M=kl;M>=-kl;M -=inc)
259 {
260 /*Debug("L =%d kl=%d M=%d \n",L,kl,M);*/
261 k++;
262 Stemp = getZlm(L,M);
263
264 temp[k].numberOfFunctions=Stemp.numberOfCoefficients*Type[GeomOrb[i].NumType].Ao[j].N;
265 temp[k].NumCenter=i;
266 /* Debug("M=%d N=%d\n",M,temp[k].N);*/
267 temp[k].Gtf =g_malloc(temp[k].numberOfFunctions*sizeof(GTF));
268 Nc=-1;
269 for(N=0;N<Type[GeomOrb[i].NumType].Ao[j].N;N++)
270 for(n=0;n<Stemp.numberOfCoefficients;n++)
271 {
272 Nc++;
273 temp[k].Gtf[Nc].Ex = Type[GeomOrb[i].NumType].Ao[j].Ex[N];
274 temp[k].Gtf[Nc].Coef = Type[GeomOrb[i].NumType].Ao[j].Coef[N]*Stemp.lxyz[n].Coef;
275 for(c=0;c<3;c++)
276 {
277 temp[k].Gtf[Nc].C[c] = GeomOrb[i].C[c];
278 temp[k].Gtf[Nc].l[c] = Stemp.lxyz[n].l[c];
279 }
280 }
281 if(L==0) break;
282 }
283 if(L==0) break;
284 }
285 }
286 for(i=0;i<NAOrb;i++) g_free(AOrb[i].Gtf);
287 g_free(AOrb);
288 NAOrb = NOrb;
289 AOrb = temp;
290 if(SAOrb) g_free(SAOrb);
291 SAOrb = NULL;
292 DefineAtomicNumOrb();
293 }
294
295 /********************************************************************************/
DefineNBOCartBasis()296 static void DefineNBOCartBasis()
297 {
298 gint i,j,k,n;
299 gint l1,l2,l3;
300 gint L;
301 gint *l[3]={NULL,NULL,NULL};
302 gint m;
303
304 NAOrb = 0;
305 for(i=0;i<nCenters;i++)
306 {
307 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
308 {
309 L=Type[GeomOrb[i].NumType].Ao[j].L;
310 NAOrb += (L+1)*(L+2)/2;
311 }
312 }
313
314 AOrb = g_malloc(NAOrb*sizeof(CGTF));
315 if(SAOrb) g_free(SAOrb);
316 SAOrb = NULL;
317
318 k=-1;
319 for(i=0;i<nCenters;i++)
320 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
321 {
322 L = Type[GeomOrb[i].NumType].Ao[j].L;
323
324 for(m=0;m<3;m++)
325 {
326 if(l[m])
327 g_free(l[m]);
328 l[m] = g_malloc((L+1)*(L+2)/2*sizeof(gint));
329 }
330 switch(L)
331 {
332 case 1 :
333 m=0;
334 l[0][m] = 1;l[1][m] = 0;l[2][m] = 0; /* X */
335 m++;
336 l[0][m] = 0;l[1][m] = 1;l[2][m] = 0; /* Y */
337 m++;
338 l[0][m] = 0;l[1][m] = 0;l[2][m] = 1; /* Z */
339 break;
340 case 2 :
341 m=0;
342 l[0][m] = 2;l[1][m] = 0;l[2][m] = 0; /* XX */
343 m++;
344 l[0][m] = 0;l[1][m] = 2;l[2][m] = 0; /* YY */
345 m++;
346 l[0][m] = 0;l[1][m] = 0;l[2][m] = 2; /* ZZ */
347 m++;
348 l[0][m] = 1;l[1][m] = 1;l[2][m] = 0; /* XY */
349 m++;
350 l[0][m] = 1;l[1][m] = 0;l[2][m] = 1; /* XZ */
351 m++;
352 l[0][m] = 0;l[1][m] = 1;l[2][m] = 1; /* YZ */
353 break;
354 case 3 :
355 m=0;
356 l[0][m] = 3;l[1][m] = 0;l[2][m] = 0; /* XXX */
357 m++;
358 l[0][m] = 0;l[1][m] = 3;l[2][m] = 0; /* YYY */
359 m++;
360 l[0][m] = 0;l[1][m] = 0;l[2][m] = 3; /* ZZZ */
361 m++;
362 l[0][m] = 2;l[1][m] = 1;l[2][m] = 0; /* XXY */
363 m++;
364 l[0][m] = 2;l[1][m] = 0;l[2][m] = 1; /* XXZ */
365 m++;
366 l[0][m] = 1;l[1][m] = 2;l[2][m] = 0; /* YYX */
367 m++;
368 l[0][m] = 0;l[1][m] = 2;l[2][m] = 1; /* YYZ */
369 m++;
370 l[0][m] = 1;l[1][m] = 0;l[2][m] = 2; /* ZZX */
371 m++;
372 l[0][m] = 0;l[1][m] = 1;l[2][m] = 2; /* ZZY */
373 m++;
374 l[0][m] = 1;l[1][m] = 1;l[2][m] = 1; /* XYZ */
375 break;
376 case 4 :
377 m=0; l[0][m] = 4;l[1][m] = 0;l[2][m] = 0; /* XXXX */
378 m++; l[0][m] = 0;l[1][m] = 4;l[2][m] = 0; /* YYYY */
379 m++; l[0][m] = 0;l[1][m] = 0;l[2][m] = 4; /* ZZZZ */
380 m++; l[0][m] = 3;l[1][m] = 1;l[2][m] = 0; /* XXXY */
381 m++; l[0][m] = 3;l[1][m] = 0;l[2][m] = 1; /* XXXZ */
382 m++; l[0][m] = 1;l[1][m] = 3;l[2][m] = 0; /* YYYX */
383 m++; l[0][m] = 0;l[1][m] = 3;l[2][m] = 1; /* YYYZ */
384 m++; l[0][m] = 1;l[1][m] = 0;l[2][m] = 3; /* ZZZX */
385 m++; l[0][m] = 0;l[1][m] = 1;l[2][m] = 3; /* ZZZY */
386 m++; l[0][m] = 2;l[1][m] = 2;l[2][m] = 0; /* XXYY */
387 m++; l[0][m] = 2;l[1][m] = 0;l[2][m] = 2; /* XXZZ */
388 m++; l[0][m] = 0;l[1][m] = 2;l[2][m] = 2; /* YYZZ */
389 m++; l[0][m] = 2;l[1][m] = 1;l[2][m] = 1; /* XXYZ */
390 m++; l[0][m] = 1;l[1][m] = 2;l[2][m] = 1; /* YYXZ */
391 m++; l[0][m] = 1;l[1][m] = 1;l[2][m] = 2; /* ZZXY */
392 break;
393 default :
394 m=0;
395 for(l3=Type[GeomOrb[i].NumType].Ao[j].L;l3>=0;l3--)
396 for(l2=Type[GeomOrb[i].NumType].Ao[j].L-l3;l2>=0;l2--)
397 {
398 l1 = Type[GeomOrb[i].NumType].Ao[j].L-l2-l3;
399 l[0][m] = l1;
400 l[1][m] = l2;
401 l[2][m] = l3;
402 m++;
403 }
404 }
405 for(m=0;m<(L+1)*(L+2)/2;m++)
406 {
407 l1 = l[0][m];
408 l2 = l[1][m];
409 l3 = l[2][m];
410 k++;
411 AOrb[k].numberOfFunctions=Type[GeomOrb[i].NumType].Ao[j].N;
412 AOrb[k].NumCenter = i;
413 AOrb[k].Gtf =g_malloc(AOrb[k].numberOfFunctions*sizeof(GTF));
414 for(n=0;n<AOrb[k].numberOfFunctions;n++)
415 {
416 AOrb[k].Gtf[n].Ex = Type[GeomOrb[i].NumType].Ao[j].Ex[n];
417 AOrb[k].Gtf[n].Coef = Type[GeomOrb[i].NumType].Ao[j].Coef[n];
418 AOrb[k].Gtf[n].C[0] = GeomOrb[i].C[0];
419 AOrb[k].Gtf[n].C[1] = GeomOrb[i].C[1];
420 AOrb[k].Gtf[n].C[2] = GeomOrb[i].C[2];
421 AOrb[k].Gtf[n].l[0] = l1;
422 AOrb[k].Gtf[n].l[1] = l2;
423 AOrb[k].Gtf[n].l[2] = l3;
424 }
425
426 }
427 }
428
429 NOrb = NAOrb;
430 DefineAtomicNumOrb();
431 /* DefineNorb();*/
432 }
433 /********************************************************************************/
read_basis_from_a_nbo_output_file(gchar * fileName)434 static gboolean read_basis_from_a_nbo_output_file(gchar *fileName)
435 {
436 static gchar t[BSIZE];
437 FILE *file;
438 gint nAtoms = 0, nShell = 0, nExp = 0;
439 gdouble* expo = NULL;
440 gdouble** coefs = NULL;
441 gint* numCenters = NULL;
442 gint* nOrbs = NULL;
443 gint** numTypes = NULL;
444 gint* iPointers = NULL;
445 gint* nGauss = NULL;
446 gint i;
447 gint j;
448 gint k;
449 gint is;
450 gint ie;
451 gint jj;
452 gint lmax = 0;
453
454
455 if ((!fileName) || (strcmp(fileName,"") == 0))
456 {
457 Message(_("Sorry No file selected\n"),_("Error"),TRUE);
458 return FALSE;
459 }
460
461 file = FOpen(fileName, "rb");
462 if(file == NULL)
463 {
464 gchar buffer[BSIZE];
465 sprintf(buffer,_("Sorry, I can not open '%s' file\n"),fileName);
466 Message(buffer,_("Error"),TRUE);
467 return FALSE;
468 }
469 if(!goToLine(file,"--------")) return FALSE;
470 if(!fgets(t,BSIZE,file)) return FALSE;
471 sscanf(t,"%d %d %d",&nAtoms,&nShell,&nExp);
472 if(nAtoms!=nCenters || nAtoms <= 0 || nShell <= 0 || nExp <= 0) return FALSE;
473 if(!goToLine(file,"--------")) return FALSE;
474 if(!goToLine(file,"--------")) return FALSE;
475
476 numCenters = g_malloc(nShell*sizeof(gint));
477 nOrbs = g_malloc(nShell*sizeof(gint));
478 numTypes = g_malloc(nShell*sizeof(gint*));
479 for(j=0;j<nShell;j++) numTypes[j] = g_malloc(nShell*sizeof(gint));
480 iPointers = g_malloc(nShell*sizeof(gint));
481 nGauss = g_malloc(nShell*sizeof(gint));
482
483 lmax = 0;
484 for(is = 0; is<nShell; is++)
485 {
486 if(1!=fscanf(file,"%d",&numCenters[is])) break;
487 numCenters[is]--;
488 if(1!=fscanf(file,"%d",&nOrbs[is])) break;
489 if(1!=fscanf(file,"%d",&iPointers[is])) break;
490 if(1!=fscanf(file,"%d",&nGauss[is])) break;
491 for(k=0;k<nOrbs[is];k++)
492 {
493 if(1!=fscanf(file,"%d",&numTypes[is][k]))break;
494 gint l = numTypes[is][k]/100;
495 if(lmax<l) lmax = l;
496 }
497 if(k!=nOrbs[is]) break;
498 }
499 if(is!=nShell)
500 {
501 if(numCenters) g_free(numCenters);
502 if(nOrbs) g_free(nOrbs);
503 if(numTypes) for(j=0;j<nShell;j++) if(numTypes[j]) g_free(numTypes[j]);
504 if(numTypes) g_free(numTypes);
505 if(iPointers) g_free(iPointers);
506 if(nGauss) g_free(nGauss);
507 return FALSE;
508 }
509 if(!goToLine(file,"--------"))
510 {
511 if(numCenters) g_free(numCenters);
512 if(nOrbs) g_free(nOrbs);
513 if(numTypes) for(j=0;j<nShell;j++) if(numTypes[j]) g_free(numTypes[j]);
514 if(numTypes) g_free(numTypes);
515 if(iPointers) g_free(iPointers);
516 if(nGauss) g_free(nGauss);
517 return FALSE;
518 }
519
520 expo = g_malloc(nExp*sizeof(gdouble));
521 coefs = g_malloc((lmax+1)*sizeof(gdouble*));
522 for(i=0;i<=lmax;i++)
523 {
524 coefs[i] = g_malloc(nExp*sizeof(gdouble));
525 for(j=0;j<nExp;j++) coefs[i][j] = 0.0;
526 }
527 for(ie = 0; ie<nExp; ie++)
528 {
529 {int ii = fscanf(file,"%lf",&expo[ie]);}
530 }
531 if(!fgets(t,BSIZE,file))
532 {
533 if(numCenters) g_free(numCenters);
534 if(nOrbs) g_free(nOrbs);
535 if(numTypes) for(j=0;j<nShell;j++) if(numTypes[j]) g_free(numTypes[j]);
536 if(numTypes) g_free(numTypes);
537 if(iPointers) g_free(iPointers);
538 if(nGauss) g_free(nGauss);
539 if(expo) g_free(expo);
540 if(coefs)
541 {
542 for(i=0;i<=lmax;i++) if(coefs[i]) g_free(coefs[i]);
543 g_free(coefs);
544 }
545 return FALSE;
546 }
547
548 for(i=0;i<=lmax;i++)
549 {
550 for(ie = 0; ie<nExp; ie++)
551 if(1!=fscanf(file,"%lf",&coefs[i][ie])) break;
552 if(!fgets(t,BSIZE,file)) break; /* f orb is not always available */
553 }
554 for(i=0;i<nCenters;i++)
555 {
556 Type[i].Norb = 0;
557 Type[i].Ao = NULL;
558 }
559 for(is = 0; is<nShell; is++)
560 {
561 i = numCenters[is];
562 Type[i].Norb++;
563
564 for(k=1;k<nOrbs[is];k++)
565 if(numTypes[is][k]/100 != numTypes[is][k-1]/100) Type[i].Norb++;
566 }
567
568 for(i=0;i<nCenters;i++)
569 {
570 Type[i].Ao=g_malloc(Type[i].Norb*sizeof(AO));
571 for(j=0;j< Type[i].Norb;j++)
572 {
573 Type[i].Ao[j].Ex = NULL;
574 Type[i].Ao[j].Coef = NULL;
575 }
576 j = 0;
577 for(is = 0; is<nShell; is++)
578 {
579 gint ncont = nGauss[is];
580 if(i != numCenters[is]) continue;
581 for(k=0;k<nOrbs[is];k++)
582 {
583 gint l = numTypes[is][k]/100;
584 if(k>0 && numTypes[is][k]/100 == numTypes[is][k-1]/100) continue;
585 Type[i].Ao[j].L = l;
586 Type[i].Ao[j].N = ncont;
587 Type[i].Ao[j].Ex=g_malloc(Type[i].Ao[j].N*sizeof(gdouble));
588 Type[i].Ao[j].Coef=g_malloc(Type[i].Ao[j].N*sizeof(gdouble));
589 for(jj=0;jj<Type[i].Ao[j].N;jj++)
590 {
591 gint jjP = iPointers[is] - 1 + jj;
592 Type[i].Ao[j].Ex[jj] = expo[jjP];
593 Type[i].Ao[j].Coef[jj] = coefs[l][jjP];
594 }
595 j++;
596 }
597 }
598
599 }
600 /* free tables */
601 if(numCenters) g_free(numCenters);
602 if(nOrbs) g_free(nOrbs);
603 if(expo) g_free(expo);
604 if(coefs)
605 {
606 for(i=0;i<=lmax;i++) if(coefs[i]) g_free(coefs[i]);
607 g_free(coefs);
608 }
609
610 if(numTypes)
611 for(j=0;j<nShell;j++) if(numTypes[j]) g_free(numTypes[j]);
612 if(numTypes) g_free(numTypes);
613 if(iPointers) g_free(iPointers);
614 if(nGauss) g_free(nGauss);
615
616 fclose(file);
617 return TRUE;
618 }
619 /********************************************************************************/
620 /* typeobr = ALPHA or BETA */
read_orbitals_in_nbo_file_alpha_or_beta(gchar * fileName,gint nAlpha,gint nBeta,gchar * typeOrb)621 static gboolean read_orbitals_in_nbo_file_alpha_or_beta(gchar *fileName,gint nAlpha, gint nBeta, gchar* typeOrb)
622 {
623 FILE *file;
624 gint i;
625 gint k;
626 gint j;
627 gdouble **CoefOrbitals;
628 gdouble *EnerOrbitals;
629 gchar **SymOrbitals;
630 gboolean om = FALSE;
631 gdouble dum;
632
633 if ((!fileName) || (strcmp(fileName,"") == 0))
634 {
635 Message(_("Sorry No file selected\n"),_("Error"),TRUE);
636 return FALSE;
637 }
638
639 file = FOpen(fileName, "rb");
640 if(file ==NULL)
641 {
642 gchar buffer[BSIZE];
643 sprintf(buffer,_("Sorry, I can not open '%s' file\n"),fileName);
644 Message(buffer,_("Error"),TRUE);
645 return FALSE;
646 }
647 {
648 gchar t[BSIZE];
649 { char* e = fgets(t,BSIZE,file);}
650 { char* e = fgets(t,BSIZE,file);}
651 if(strstr(t," MOs ")) om = TRUE;
652 }
653
654 if(!goToLine(file,"--------")) return FALSE;
655 if(!goToLine(file,typeOrb)) return FALSE;
656
657
658 CoefOrbitals = CreateTable2(NOrb);
659 EnerOrbitals = g_malloc(NOrb*sizeof(gdouble));
660 SymOrbitals = g_malloc(NOrb*sizeof(gchar*));
661
662 for(j=0;j<NOrb;j++)
663 {
664 EnerOrbitals[j]=0.0;
665 for(i=0;i<NOrb;i++)
666 {
667 if(1!=fscanf(file,"%lf",&dum))break;
668 if(om) CoefOrbitals[i][j] = dum;
669 else CoefOrbitals[j][i] = dum;
670 }
671 if(i!=NOrb)
672 {
673 for(k=j;k<NOrb;k++) SymOrbitals[k] = g_strdup("DELETED");
674 break;
675 }
676 SymOrbitals[j] = g_strdup("UNK");
677 }
678 /*
679 printf("norb read=%d\n",j);
680 */
681 if(strstr(typeOrb,"ALPHA"))
682 {
683 CoefAlphaOrbitals = CoefOrbitals;
684 EnerAlphaOrbitals = EnerOrbitals;
685 SymAlphaOrbitals = SymOrbitals;
686 OccAlphaOrbitals = g_malloc(NOrb*sizeof(gdouble));
687 for(i=0;i<nAlpha;i++) OccAlphaOrbitals[i] = 1.0;
688 for(i=nAlpha;i<NOrb;i++) OccAlphaOrbitals[i] = 0.0;
689 }
690 else{
691 CoefBetaOrbitals = CoefOrbitals;
692 EnerBetaOrbitals = EnerOrbitals;
693 OccBetaOrbitals = OccAlphaOrbitals;
694 SymBetaOrbitals = SymOrbitals;
695 NAlphaOcc = nAlpha;
696 NBetaOcc = nBeta;
697 NAlphaOrb = NOrb;
698 NBetaOrb = NOrb;
699 }
700 setTitle(file);
701 fclose(file);
702
703 return TRUE;
704 }
705 /********************************************************************************/
read_orbitals_in_nbo_file_all(gchar * fileName,gint nAlpha,gint nBeta)706 static gboolean read_orbitals_in_nbo_file_all(gchar *fileName,gint nAlpha, gint nBeta)
707 {
708 FILE *file;
709 gint i;
710 gint k;
711 gint j;
712 gdouble **CoefOrbitals;
713 gdouble *EnerOrbitals;
714 gchar **SymOrbitals;
715 gboolean om = FALSE;
716 gdouble dum;
717
718 if ((!fileName) || (strcmp(fileName,"") == 0))
719 {
720 Message(_("Sorry No file selected\n"),_("Error"),TRUE);
721 return FALSE;
722 }
723
724 file = FOpen(fileName, "rb");
725 if(file ==NULL)
726 {
727 gchar buffer[BSIZE];
728 sprintf(buffer,_("Sorry, I can not open '%s' file\n"),fileName);
729 Message(buffer,_("Error"),TRUE);
730 return FALSE;
731 }
732 {
733 gchar t[BSIZE];
734 { char* e = fgets(t,BSIZE,file);}
735 { char* e = fgets(t,BSIZE,file);}
736 if(strstr(t," MOs ")) om = TRUE;
737 }
738 if(!goToLine(file,"--------")) return FALSE;
739
740
741 CoefOrbitals = CreateTable2(NOrb);
742 EnerOrbitals = g_malloc(NOrb*sizeof(gdouble));
743 SymOrbitals = g_malloc(NOrb*sizeof(gchar*));
744
745 for(j=0;j<NOrb;j++)
746 {
747 EnerOrbitals[j]=0.0;
748 for(i=0;i<NOrb;i++)
749 {
750 if(1!=fscanf(file,"%lf",&dum))break;
751 if(om) CoefOrbitals[i][j] = dum;
752 else CoefOrbitals[j][i] = dum;
753 }
754 if(i!=NOrb)
755 {
756 for(k=j;k<NOrb;k++) SymOrbitals[k] = g_strdup("DELETED");
757 break;
758 }
759 SymOrbitals[j] = g_strdup("UNK");
760 }
761 /*
762 printf("norb read=%d\n",j);
763 */
764 CoefAlphaOrbitals = CoefOrbitals;
765 EnerAlphaOrbitals = EnerOrbitals;
766 SymAlphaOrbitals = SymOrbitals;
767 OccAlphaOrbitals = g_malloc(NOrb*sizeof(gdouble));
768 nAlpha = nBeta;
769 for(i=0;i<nAlpha;i++) OccAlphaOrbitals[i] = 1.0;
770 for(i=nAlpha;i<NOrb;i++) OccAlphaOrbitals[i] = 0.0;
771
772 CoefBetaOrbitals = CoefOrbitals;
773 EnerBetaOrbitals = EnerOrbitals;
774 OccBetaOrbitals = OccAlphaOrbitals;
775 SymBetaOrbitals = SymOrbitals;
776 NAlphaOcc = nAlpha;
777 NBetaOcc = nBeta;
778 NAlphaOrb = NOrb;
779 NBetaOrb = NOrb;
780 setTitle(file);
781 fclose(file);
782
783 return TRUE;
784 }
785 /********************************************************************************/
read_orbitals_in_nbo_file(gchar * fileName,gint nAlpha,gint nBeta)786 static gboolean read_orbitals_in_nbo_file(gchar *fileName,gint nAlpha, gint nBeta)
787 {
788 if(read_orbitals_in_nbo_file_alpha_or_beta(fileName, nAlpha, nBeta, "ALPHA"))
789 {
790 return read_orbitals_in_nbo_file_alpha_or_beta(fileName, nAlpha, nBeta, "BETA");
791 }
792 else
793 return read_orbitals_in_nbo_file_all(fileName, nAlpha, nBeta);
794
795 return TRUE;
796 }
797 /********************************************************************************/
read_nbo_orbitals(gchar * fileName)798 void read_nbo_orbitals(gchar* fileName)
799 {
800 gint typefile;
801 /* gint typebasis=1;*/ /* NBO print OM in cartezian presentation even ISPHER=0 or 1 or -1 */
802 gchar *t = NULL;
803 gint i;
804 gint nAlpha;
805 gint nBeta;
806 gint typebasis = 0;
807 gboolean Ok;
808 gchar* fileName31 = NULL;
809
810
811 typefile =get_type_file_orb(fileName);
812 if(typefile==GABEDIT_TYPEFILE_UNKNOWN) return;
813
814
815 if(typefile != GABEDIT_TYPEFILE_NBO)
816 {
817 gchar buffer[BSIZE];
818 sprintf(buffer,_("Sorry, I can not read this format from '%s' file\n"),fileName);
819 Message(buffer,_("Error"),TRUE);
820 return ;
821 }
822 fileName31 = g_strdup_printf("%s.31",get_suffix_name_file(fileName));
823
824 free_data_all();
825 t = get_name_file(fileName31);
826 set_status_label_info(_("File name"),t);
827 g_free(t);
828 set_status_label_info(_("File type"),"NBO");
829 set_status_label_info(_("Mol. Orb."),_("Reading"));
830
831 free_orbitals();
832
833 /* typebasis =get_type_basis_in_nbo_file(fileName);
834 if(typebasis == -1)
835 {
836 gchar buffer[BSIZE];
837 sprintf(buffer,
838 "Sorry, Gabedit does not support spherical basis with contaminant cartezian function\n\n"
839 "Use ISPHER=-1 or ISPHER=1 in CONTROL block"
840 );
841 Message(buffer,_("Error"),TRUE);
842 set_status_label_info(_("File name"),_("Nothing"));
843 set_status_label_info(_("File type"),_("Nothing"));
844 set_status_label_info(_("Mol. Orb."),_("Nothing"));
845 return;
846 }
847 */
848
849 if(!read_geomorb_nbo_file_geom(fileName31))
850 {
851 free_geometry();
852 set_status_label_info(_("File name"),_("Nothing"));
853 set_status_label_info(_("File type"),_("Nothing"));
854 set_status_label_info(_("Mol. Orb."),_("Nothing"));
855 return;
856 }
857 set_status_label_info(_("Geometry"),_("Ok"));
858 glarea_rafresh(GLArea); /* for geometry*/
859 init_atomic_orbitals();
860 if(!read_basis_from_a_nbo_output_file(fileName31))
861 {
862 gchar buffer[BSIZE];
863 sprintf(buffer,_("Sorry, I cannot read basis from '%s' file\n"),fileName31);
864 Message(buffer,_("Error"),TRUE);
865 if(GeomOrb)
866 {
867 init_atomic_orbitals();
868 for(i=0;i<nCenters;i++) GeomOrb[i].Prop = prop_atom_get("H");
869 free_geometry();
870 }
871 set_status_label_info(_("File name"),_("Nothing"));
872 set_status_label_info(_("File type"),_("Nothing"));
873 set_status_label_info(_("Mol. Orb."),_("Nothing"));
874 return;
875 }
876
877 g_free(fileName31);
878
879 t = get_name_file(fileName);
880 set_status_label_info(_("File name"),t);
881 g_free(t);
882 set_status_label_info(_("Mol. Orb."),_("Reading"));
883 InitializeAll();
884 buildBondsOrb();
885 RebuildGeom = TRUE;
886 reset_grid_limits();
887 init_atomic_orbitals();
888
889
890 if(typebasis == 0)
891 {
892 DefineNBOSphericalBasis();
893 sphericalBasis = TRUE;
894 }
895 else
896 {
897 DefineNBOCartBasis();
898 sphericalBasis = FALSE;
899 }
900
901 /*
902 printf("Not normalized basis\n");
903 PrintAllBasis();
904 */
905 /* NormaliseAllBasis();*/
906 NormaliseAllNoRadBasis();
907 /*
908 printf("Normalized basis\n");
909 PrintAllBasis();
910 */
911 DefineNOccs();
912
913 nAlpha = NOrb;
914 nBeta = NOrb;
915 nAlpha = NAlphaOcc;
916 nBeta = NBetaOcc;
917
918 /*
919 printf("Number of ALPHA occ = %d\n",nAlpha);
920 printf("Number of BETA occ = %d\n",nBeta);
921 printf("NOrb = %d\n",NOrb);
922 */
923 Ok = read_orbitals_in_nbo_file(fileName,nAlpha,nBeta);
924
925 if(Ok)
926 {
927 /*PrintAllOrb(CoefAlphaOrbitals);*/
928 set_status_label_info(_("Mol. Orb."),_("Ok"));
929 glarea_rafresh(GLArea); /* for geometry*/
930 NumSelOrb = NAlphaOcc-1;
931 create_list_orbitals();
932 }
933 else
934 {
935 free_orbitals();
936 set_status_label_info(_("File name"),_("Nothing"));
937 set_status_label_info(_("File type"),_("Nothing"));
938 set_status_label_info(_("Mol. Orb."),_("Nothing"));
939 }
940
941 }
942