1 /* OrbitalsNWChem.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/Constants.h"
26 #include "../Utils/Zlm.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
40 /********************************************************************************/
41 typedef enum
42 {
43 GABEDIT_ORBLOCALTYPE_BOYS=0,
44 GABEDIT_ORBLOCALTYPE_EDMISTON,
45 GABEDIT_ORBLOCALTYPE_PIPEK,
46 GABEDIT_ORBLOCALTYPE_UNKNOWN
47 } GabEditOrbLocalType;
48
49 static gchar* titlesLocalOrb[GABEDIT_ORBLOCALTYPE_PIPEK+1]=
50 {
51 "BOYS ORBITAL LOCALIZATION",
52 "EDMISTON-RUEDENBERG ENERGY LOCALIZATION",
53 "MOLECULAR ORBITALS LOCALIZED BY THE POPULATION METHOD"
54 };
55
56 typedef enum
57 {
58 GABEDIT_ORBTYPE_ALPHA = 0,
59 GABEDIT_ORBTYPE_BETA,
60 GABEDIT_ORBTYPE_RESTRICTED,
61 GABEDIT_ORBTYPE_MCSCF,
62 GABEDIT_ORBTYPE_EIGENVECTORS,
63 GABEDIT_ORBTYPE_BOYS_ALPHA,
64 GABEDIT_ORBTYPE_BOYS_BETA,
65 GABEDIT_ORBTYPE_BOYS,
66 GABEDIT_ORBTYPE_EDMISTON_ALPHA,
67 GABEDIT_ORBTYPE_EDMISTON_BETA,
68 GABEDIT_ORBTYPE_EDMISTON,
69 GABEDIT_ORBTYPE_PIPEK_ALPHA,
70 GABEDIT_ORBTYPE_PIPEK_BETA,
71 GABEDIT_ORBTYPE_PIPEK,
72 } GabEditOrbType;
73 static gchar* titlesOrb[GABEDIT_ORBTYPE_PIPEK+1]=
74 {
75 "Final Alpha Molecular Orbital Analysis",
76 "Final Beta Molecular Orbital Analysis",
77 "Final Molecular Orbital Analysis",
78 "MCSCF OPTIMIZED ORBITALS",
79 "EIGENVECTORS",
80 "***** ALPHA ORBITAL LOCALIZATION *****",
81 "****** BETA ORBITAL LOCALIZATION *****",
82 "THE BOYS LOCALIZED ORBITALS ARE",
83 "***** ALPHA ORBITAL LOCALIZATION *****",
84 "****** BETA ORBITAL LOCALIZATION *****",
85 "EDMISTON-RUEDENBERG ENERGY LOCALIZED ORBITALS",
86 "***** ALPHA ORBITAL LOCALIZATION *****",
87 "****** BETA ORBITAL LOCALIZATION *****",
88 "THE PIPEK-MEZEY POPULATION LOCALIZED ORBITALS ARE"
89 };
90 /********************************************************************************/
91 static gboolean sphericalBasis = FALSE;
92 /********************************************************************************/
read_geomorb_nwchem_file_geom(gchar * FileName)93 static gulong read_geomorb_nwchem_file_geom(gchar *FileName)
94 {
95 gulong lineg = gl_read_nwchem_file_geomi(FileName,-1);
96 if(nCenters != 0 )
97 {
98 DefineType();
99 /* PrintGeomOrb();*/
100 }
101 return lineg;
102 }
103 /********************************************************************************/
DefineNWChemCartBasis()104 static void DefineNWChemCartBasis()
105 {
106 gint i,j,k,n;
107 gint l1,l2,l3;
108 gint L;
109 gint *l[3]={NULL,NULL,NULL};
110 gint m;
111
112 NAOrb = 0;
113 for(i=0;i<nCenters;i++)
114 {
115 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
116 {
117 L=Type[GeomOrb[i].NumType].Ao[j].L;
118 NAOrb += (L+1)*(L+2)/2;
119 }
120 }
121
122 AOrb = g_malloc(NAOrb*sizeof(CGTF));
123 if(SAOrb) g_free(SAOrb);
124 SAOrb = NULL;
125
126 k=-1;
127 for(i=0;i<nCenters;i++)
128 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
129 {
130 L = Type[GeomOrb[i].NumType].Ao[j].L;
131
132 for(m=0;m<3;m++)
133 {
134 if(l[m])
135 g_free(l[m]);
136 l[m] = g_malloc((L+1)*(L+2)/2*sizeof(gint));
137 }
138 switch(L)
139 {
140 case 1 :
141 m=0;
142 l[0][m] = 1;l[1][m] = 0;l[2][m] = 0; /* X */
143 m++;
144 l[0][m] = 0;l[1][m] = 1;l[2][m] = 0; /* Y */
145 m++;
146 l[0][m] = 0;l[1][m] = 0;l[2][m] = 1; /* Z */
147 break;
148 case 2 :
149 m=0; l[0][m] = 2;l[1][m] = 0;l[2][m] = 0; /* XX */
150 m++; l[0][m] = 1;l[1][m] = 1;l[2][m] = 0; /* XY */
151 m++; l[0][m] = 0;l[1][m] = 2;l[2][m] = 0; /* YY */
152 m++; l[0][m] = 1;l[1][m] = 0;l[2][m] = 1; /* XZ */
153 m++; l[0][m] = 0;l[1][m] = 1;l[2][m] = 1; /* YZ */
154 m++; l[0][m] = 0;l[1][m] = 0;l[2][m] = 2; /* ZZ */
155 break;
156 case 3 :
157 m=0; l[0][m] = 3;l[1][m] = 0;l[2][m] = 0; /* XXX */
158 m++; l[0][m] = 2;l[1][m] = 1;l[2][m] = 0; /* XXY */
159 m++; l[0][m] = 1;l[1][m] = 2;l[2][m] = 0; /* YYX */
160 m++; l[0][m] = 0;l[1][m] = 3;l[2][m] = 0; /* YYY */
161 m++; l[0][m] = 2;l[1][m] = 0;l[2][m] = 1; /* XXZ */
162 m++; l[0][m] = 1;l[1][m] = 1;l[2][m] = 1; /* XYZ */
163 m++; l[0][m] = 0;l[1][m] = 2;l[2][m] = 1; /* YYZ */
164 m++; l[0][m] = 1;l[1][m] = 0;l[2][m] = 2; /* ZZX */
165 m++; l[0][m] = 0;l[1][m] = 1;l[2][m] = 2; /* ZZY */
166 m++; l[0][m] = 0;l[1][m] = 0;l[2][m] = 3; /* ZZZ */
167 break;
168 case 4 :
169 /* I think that Qchem can n ot print orb with G cartz */
170 m=0; l[0][m] = 4;l[1][m] = 0;l[2][m] = 0; /* XXXX */
171 m++; l[0][m] = 0;l[1][m] = 4;l[2][m] = 0; /* YYYY */
172 m++; l[0][m] = 0;l[1][m] = 0;l[2][m] = 4; /* ZZZZ */
173 m++; l[0][m] = 3;l[1][m] = 1;l[2][m] = 0; /* XXXY */
174 m++; l[0][m] = 3;l[1][m] = 0;l[2][m] = 1; /* XXXZ */
175 m++; l[0][m] = 1;l[1][m] = 3;l[2][m] = 0; /* YYYX */
176 m++; l[0][m] = 0;l[1][m] = 3;l[2][m] = 1; /* YYYZ */
177 m++; l[0][m] = 1;l[1][m] = 0;l[2][m] = 3; /* ZZZX */
178 m++; l[0][m] = 0;l[1][m] = 1;l[2][m] = 3; /* ZZZY */
179 m++; l[0][m] = 2;l[1][m] = 2;l[2][m] = 0; /* XXYY */
180 m++; l[0][m] = 2;l[1][m] = 0;l[2][m] = 2; /* XXZZ */
181 m++; l[0][m] = 0;l[1][m] = 2;l[2][m] = 2; /* YYZZ */
182 m++; l[0][m] = 2;l[1][m] = 1;l[2][m] = 1; /* XXYZ */
183 m++; l[0][m] = 1;l[1][m] = 2;l[2][m] = 1; /* YYXZ */
184 m++; l[0][m] = 1;l[1][m] = 1;l[2][m] = 2; /* ZZXY */
185 break;
186 default :
187 m=0;
188 for(l3=Type[GeomOrb[i].NumType].Ao[j].L;l3>=0;l3--)
189 for(l2=Type[GeomOrb[i].NumType].Ao[j].L-l3;l2>=0;l2--)
190 {
191 l1 = Type[GeomOrb[i].NumType].Ao[j].L-l2-l3;
192 l[0][m] = l1;
193 l[1][m] = l2;
194 l[2][m] = l3;
195 m++;
196 }
197 }
198 for(m=0;m<(L+1)*(L+2)/2;m++)
199 {
200 l1 = l[0][m];
201 l2 = l[1][m];
202 l3 = l[2][m];
203 k++;
204 AOrb[k].numberOfFunctions=Type[GeomOrb[i].NumType].Ao[j].N;
205 AOrb[k].NumCenter = i;
206 AOrb[k].Gtf =g_malloc(AOrb[k].numberOfFunctions*sizeof(GTF));
207 for(n=0;n<AOrb[k].numberOfFunctions;n++)
208 {
209 AOrb[k].Gtf[n].Ex = Type[GeomOrb[i].NumType].Ao[j].Ex[n];
210 AOrb[k].Gtf[n].Coef = Type[GeomOrb[i].NumType].Ao[j].Coef[n];
211 AOrb[k].Gtf[n].C[0] = GeomOrb[i].C[0];
212 AOrb[k].Gtf[n].C[1] = GeomOrb[i].C[1];
213 AOrb[k].Gtf[n].C[2] = GeomOrb[i].C[2];
214 AOrb[k].Gtf[n].l[0] = l1;
215 AOrb[k].Gtf[n].l[1] = l2;
216 AOrb[k].Gtf[n].l[2] = l3;
217 }
218
219 }
220 }
221
222 NOrb = NAOrb;
223 DefineAtomicNumOrb();
224 /* DefineNorb();*/
225 }
226 /**********************************************/
DefineNWChemSphericalBasis()227 static void DefineNWChemSphericalBasis()
228 {
229 gint i,j,k;
230 gint c;
231 gint kl;
232 gint L,M;
233 CGTF *temp;
234 Zlm Stemp;
235 gint N,Nc,n;
236 gint klbeg;
237 gint klend;
238 gint klinc;
239
240
241 NOrb = 0;
242 for(i=0;i<nCenters;i++)
243 {
244 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
245 {
246 L=Type[GeomOrb[i].NumType].Ao[j].L;
247 NOrb += 2*L+1;
248 }
249 }
250
251 temp = g_malloc(NOrb*sizeof(CGTF));
252
253 k=-1;
254 for(i=0;i<nCenters;i++)
255 for(j=0;j<Type[GeomOrb[i].NumType].Norb;j++)
256 {
257 L =Type[GeomOrb[i].NumType].Ao[j].L;
258
259 if(L==1)
260 {
261 klbeg = L;
262 klend = -L;
263 klinc = -1;
264 }
265 else
266 {
267 klbeg = -L;
268 klend = L;
269 klinc = +1;
270 }
271 for(kl = klbeg;(klbeg == -L && kl<=klend) || (klbeg == L && kl>=klend);kl +=klinc)
272 {
273 M = kl;
274 k++;
275 Stemp = getZlm(L,M);
276
277 temp[k].numberOfFunctions=Stemp.numberOfCoefficients*Type[GeomOrb[i].NumType].Ao[j].N;
278 temp[k].NumCenter=i;
279 temp[k].Gtf =g_malloc(temp[k].numberOfFunctions*sizeof(GTF));
280 Nc=-1;
281 for(N=0;N<Type[GeomOrb[i].NumType].Ao[j].N;N++)
282 for(n=0;n<Stemp.numberOfCoefficients;n++)
283 {
284 Nc++;
285
286 temp[k].Gtf[Nc].Ex = Type[GeomOrb[i].NumType].Ao[j].Ex[N];
287 temp[k].Gtf[Nc].Coef = Type[GeomOrb[i].NumType].Ao[j].Coef[N]*Stemp.lxyz[n].Coef;
288 for(c=0;c<3;c++)
289 {
290 temp[k].Gtf[Nc].C[c] = GeomOrb[i].C[c];
291 temp[k].Gtf[Nc].l[c] = Stemp.lxyz[n].l[c];
292 }
293 }
294 if(L==0) break;
295 if(L==1 && kl == -1)
296 {
297 /* switch -1 and 0 */
298 CGTF dum = temp[k];
299 temp[k] = temp[k-1];
300 temp[k-1] = dum;
301 }
302 }
303 }
304 for(i=0;i<NAOrb;i++)
305 g_free(AOrb[i].Gtf);
306 g_free(AOrb);
307 NAOrb = NOrb;
308 AOrb = temp;
309 if(SAOrb) g_free(SAOrb);
310 SAOrb = NULL;
311 DefineAtomicNumOrb();
312 }
313 /********************************************************************************/
read_basis_from_a_nwchem_output_file(gchar * FileName,gint * nrs,gulong lineg)314 static gchar** read_basis_from_a_nwchem_output_file(gchar *FileName, gint* nrs, gulong lineg)
315 {
316 gchar **strbasis;
317 gchar *t;
318 FILE *fd;
319 gint nrows=0;
320 gboolean OK = FALSE;
321 gulong line = 0;
322
323 if ((!FileName) || (strcmp(FileName,"") == 0))
324 {
325 Message(_("Sorry No file selected\n"),_("Error"),TRUE);
326 return NULL;
327 }
328
329 fd = FOpen(FileName, "rb");
330 if(fd ==NULL)
331 {
332 gchar buffer[BSIZE];
333 sprintf(buffer,_("Sorry, I can not open '%s' file\n"),FileName);
334 Message(buffer,_("Error"),TRUE);
335 return NULL;
336 }
337 t=g_malloc(BSIZE*sizeof(gchar));
338 while(!feof(fd))
339 {
340 line++;
341 if(!fgets(t,BSIZE,fd))break;
342 if(line==lineg) break;
343 }
344 if(line!=lineg)
345 {
346 gchar buffer[BSIZE];
347 sprintf(buffer,_("Sorry, No basis in '%s' file after the last geometry\n"),FileName);
348 Message(buffer,_("Error"),TRUE);
349 if(t) g_free(t);
350 fclose(fd);
351 return NULL;
352 }
353
354 while(!feof(fd))
355 {
356 if(!fgets(t,BSIZE,fd))break;
357 if(strstr(t,"Basis")&&strstr(t,"ao basis"))
358 {
359 if(!fgets(t,BSIZE,fd))break;
360 if(!strstr(t,"---")) continue;
361 OK = TRUE;
362 break;
363 }
364 }
365 if(!OK)
366 {
367 g_free(t);
368 Message(_("Sorry I can read basis from your NWChem output file."),_("Error"),TRUE);
369 return NULL;
370 }
371
372 strbasis=g_malloc(sizeof(gchar*));;
373 while(!feof(fd))
374 {
375 if(!fgets(t,BSIZE,fd))break;
376 if(strstr(t,"Summary of") ) break;
377 nrows++;
378 strbasis = g_realloc(strbasis,nrows*sizeof(gchar*));
379 strbasis[nrows-1] = g_strdup(t);
380 }
381 if(nrows == 0)
382 {
383 g_free(t);
384 g_free(strbasis);
385 Message(_("Sorry I can read basis from this file\n"),_("Error"),TRUE);
386 return NULL;
387 }
388
389 /*
390 Debug("End of read \n");
391 Debug("Atomic basis nrows = %d \n",nrows);
392 for(i=0;i<nrows;i++)
393 {
394 Debug("%s",strbasis[i]);
395 }
396 */
397
398 fclose(fd);
399 g_free(t);
400 *nrs = nrows;
401 return strbasis;
402 }
403 /**********************************************/
get_num_type_from_symbol(gchar * symbol)404 static gint get_num_type_from_symbol(gchar* symbol)
405 {
406 gint k;
407 for(k=0;k<nCenters;k++)
408 {
409 if(strcmp(symbol,GeomOrb[k].Symb)==0)
410 return (gint)GeomOrb[k].NumType;
411 }
412 return -1;
413 }
414 /**********************************************/
addOneBasis(gint i,gint j,gchar * shell,gint ncont,gdouble * ex,gdouble * coef)415 static gboolean addOneBasis(gint i,gint j,gchar *shell,gint ncont, gdouble* ex, gdouble* coef)
416 {
417 gint jj;
418 Type[i].Ao[j].N = ncont;
419 Type[i].Ao[j].Ex=g_malloc(Type[i].Ao[j].N*sizeof(gdouble));
420 Type[i].Ao[j].Coef=g_malloc(Type[i].Ao[j].N*sizeof(gdouble));
421 for(jj=0;jj<Type[i].Ao[j].N;jj++)
422 {
423 Type[i].Ao[j].Ex[jj] = ex[jj];
424 Type[i].Ao[j].Coef[jj] = coef[jj];
425 }
426 switch(shell[0])
427 {
428 /* L =SP with NWChem */
429 case 'l' :
430 case 'L' :
431 case 's' :
432 case 'S' : Type[i].Ao[j].L=0;break;
433 case 'p' :
434 case 'P' : Type[i].Ao[j].L=1;break;
435 case 'd' :
436 case 'D' : Type[i].Ao[j].L=2;break;
437 case 'f' :
438 case 'F' : Type[i].Ao[j].L=3;break;
439 case 'g' :
440 case 'G' : Type[i].Ao[j].L=4;break;
441 case 'h' :
442 case 'H' : Type[i].Ao[j].L=5;break;
443 case 'i' :
444 case 'I' : Type[i].Ao[j].L=6;break;
445 case 'j' :
446 case 'J' : Type[i].Ao[j].L=7;break;
447 case 'k' :
448 case 'K' : Type[i].Ao[j].L=8;break;
449
450 default : return FALSE;
451 }
452 return TRUE;
453 }
454 /**********************************************/
DefineNWChemBasisType(gchar ** strbasis,gint nrows)455 static gboolean DefineNWChemBasisType(gchar** strbasis, gint nrows)
456 {
457 gchar sym[50];
458 gchar shell[10];
459 gchar t[10];
460 gint i;
461 gint j;
462 gint nconts;
463 gint k;
464 gdouble *ex=NULL;
465 gdouble *coef1=NULL;
466 gdouble *coef2=NULL;
467 gchar* temp[10];
468 gint ne;
469 gboolean Ok;
470 gint jj;
471 /* gboolean newAtom = TRUE;*/
472
473 if(Ntype<1) return FALSE;
474 if(nrows<1) return FALSE;
475 ex = g_malloc(nrows*sizeof(gdouble));
476 coef1 = g_malloc(nrows*sizeof(gdouble));
477 coef2 = g_malloc(nrows*sizeof(gdouble));
478 for(i=0;i<10;i++) temp[i] = g_malloc(BSIZE*sizeof(gchar));
479
480 /*
481 for(k=0;k<nCenters;k++)
482 {
483 printf("%s %d\n",GeomOrb[k].Symb,GeomOrb[k].NumType);
484 }
485 */
486
487 Type = g_malloc(Ntype*sizeof(TYPE));
488 for(i=0;i<Ntype;i++)
489 {
490 Type[i].Ao = NULL;
491 Type[i].Norb=0;
492 }
493 for(k=0;k<nCenters;k++)
494 {
495 sprintf(sym,"%s",GeomOrb[k].Symb);
496 i = GeomOrb[k].NumType;
497 /* printf("numType = %d\n",k);*/
498 Type[i].Symb=g_strdup(sym);
499 Type[i].N=GetNelectrons(sym);
500 }
501 /* set number of basis by type */
502 i = -1;
503 sprintf(shell,"S");
504 nconts=-1;
505 Ok = TRUE;
506 nconts = 0;
507 /* newAtom = TRUE;*/
508 for(k=0;k<nrows-1;k++)
509 {
510 if(strstr(strbasis[k],"----")) continue;
511 if(strstr(strbasis[k],"Exponent")) continue;
512 sscanf(strbasis[k],"%s",t);
513 if(!this_is_a_backspace(strbasis[k]) &&!isdigit(t[0])) /* new atom */
514 {
515 /* newAtom = FALSE;*/
516 t[0] = toupper(t[0]);
517 if (2==strlen(t)) t[1]=tolower(t[1]);
518 i=get_num_type_from_symbol(t);
519 if(i<0)
520 {
521 Ok = FALSE;
522 break;
523 }
524 /* printf("symbl=%s\n",t);*/
525 Type[i].Norb=0;
526 continue;
527 }
528 if(this_is_a_backspace(strbasis[k])) continue;/* new shell */
529 nconts=0;
530 while(!this_is_a_backspace(strbasis[k]))
531 {
532 /* printf("%s ",strbasis[k]);*/
533 ne = sscanf(strbasis[k],"%s %s %s %s",temp[0],shell,temp[2],temp[3]);
534 if(ne!=4)
535 {
536 Ok = FALSE;
537 break;
538 }
539 for(j=0;j<ne;j++)
540 {
541 gchar* d=strstr(temp[j],"D");
542 if(d) *d='e';
543 }
544 ex[nconts]=atof(temp[2]);
545 coef1[nconts]=atof(temp[3]);
546 nconts++;
547 k++;
548 }
549 /* printf("%s %d\n",shell,nconts);*/
550 if(nconts != 0)
551 {
552 gint j = Type[i].Norb;
553 uppercase(shell);
554 if(strcmp(shell,"SP")==0) Type[i].Norb+=2;
555 else Type[i].Norb++;
556 if(Type[i].Ao == NULL)
557 Type[i].Ao=g_malloc(Type[i].Norb*sizeof(AO));
558 else
559 Type[i].Ao=g_realloc(Type[i].Ao,Type[i].Norb*sizeof(AO));
560 for(jj=j;jj< Type[i].Norb;jj++)
561 {
562 Type[i].Ao[jj].Ex = NULL;
563 Type[i].Ao[jj].Coef = NULL;
564 }
565
566
567 if(!addOneBasis(i,j,shell,nconts, ex, coef1))
568 {
569 Ok = FALSE;
570 break;
571 }
572 if(strcmp(shell,"SP")==0)
573 {
574 if(!addOneBasis(i,j+1,"P",nconts, ex, coef2))
575 {
576 Ok = FALSE;
577 break;
578 }
579 }
580
581 /*
582 printf("shell =%s ",shell);
583 printf("nconts =%d\n",nconts);
584 */
585 nconts=0;
586 }
587 }
588 if(!Ok)
589 {
590 if(Type)
591 for(i=0;i<Ntype;i++)
592 {
593 if(Type[i].Ao != NULL)
594 {
595 for(j=0;j<Type[i].Norb;j++)
596 {
597 if(Type[i].Ao[j].Ex != NULL) g_free(Type[i].Ao[j].Ex);
598 if(Type[i].Ao[j].Coef != NULL) g_free(Type[i].Ao[j].Coef);
599 }
600 g_free(Type[i].Ao);
601 }
602 }
603
604 if(Type) g_free(Type);
605 if(ex) g_free(ex);
606 if(coef1) g_free(coef1);
607 if(coef2) g_free(coef2);
608 for(j=0;j<10;j++) g_free(temp[j]);
609 return FALSE;
610 }
611 if(ex) g_free(ex);
612 if(coef1) g_free(coef1);
613 if(coef2) g_free(coef2);
614 for(j=0;j<10;j++) g_free(temp[j]);
615
616 return TRUE;
617 }
618 /********************************************************************************/
get_local_orbital_type(gchar * fileName)619 static GabEditOrbLocalType get_local_orbital_type(gchar *fileName)
620 {
621 gchar *t;
622 FILE *file;
623
624 if ((!fileName) || (strcmp(fileName,"") == 0)) return GABEDIT_ORBLOCALTYPE_UNKNOWN;
625
626 t=g_malloc(BSIZE);
627 file = FOpen(fileName, "rb");
628 if(file ==NULL) return GABEDIT_ORBLOCALTYPE_UNKNOWN;
629 while(!feof(file))
630 {
631 GabEditOrbLocalType i;
632 if(!fgets(t,BSIZE,file))break;
633 for(i=GABEDIT_ORBLOCALTYPE_BOYS;i<=GABEDIT_ORBLOCALTYPE_PIPEK;i++)
634 {
635 if(strstr( t,titlesLocalOrb[i]))
636 {
637 g_free(t);
638 return i;
639 }
640 }
641 }
642 g_free(t);
643 return GABEDIT_ORBLOCALTYPE_UNKNOWN;
644 }
645 /********************************************************************************/
read_last_orbitals_in_nwchem_file(gchar * fileName,GabEditOrbType itype,gulong lineg)646 static gboolean read_last_orbitals_in_nwchem_file(gchar *fileName,GabEditOrbType itype,gulong lineg)
647 {
648 gchar *t;
649 gboolean OK;
650 gchar *dum;
651 FILE *fd;
652 gint i;
653 gint numorb;
654 gchar *pdest = NULL;
655 gint j;
656 gdouble **CoefOrbitals;
657 gdouble *EnerOrbitals;
658 gdouble *OccOrbitals;
659 gchar **SymOrbitals;
660 gulong line = 0;
661
662 if ((!fileName) || (strcmp(fileName,"") == 0))
663 {
664 Message(_("Sorry No file selected\n"),_("Error"),TRUE);
665 return FALSE;
666 }
667
668 t=g_malloc(BSIZE*sizeof(gchar));
669 fd = FOpen(fileName, "rb");
670 if(fd ==NULL)
671 {
672 gchar buffer[BSIZE];
673 sprintf(buffer,_("Sorry, I can not open '%s' file\n"),fileName);
674 Message(buffer,_("Error"),TRUE);
675 return FALSE;
676 }
677 while(!feof(fd))
678 {
679 line++;
680 if(!fgets(t,BSIZE,fd))break;
681 if(line==lineg) break;
682 }
683 if(line!=lineg)
684 {
685 gchar buffer[BSIZE];
686 sprintf(buffer,_("Sorry, No orbitals in '%s' file after the last geometry\n"),fileName);
687 Message(buffer,_("Error"),TRUE);
688 if(t) g_free(t);
689 fclose(fd);
690 return FALSE;
691 }
692 dum=g_malloc(BSIZE*sizeof(gchar));
693
694 /* Debug("Norb = %d\n",NOrb);*/
695 CoefOrbitals = CreateTable2(NOrb);
696 EnerOrbitals = g_malloc(NOrb*sizeof(gdouble));
697 OccOrbitals = g_malloc(NOrb*sizeof(gdouble));
698 for(i=0;i<NOrb;i++) OccOrbitals[i] = 0.0;
699 SymOrbitals = g_malloc(NOrb*sizeof(gchar*));
700 numorb =1;
701 do
702 {
703 OK=FALSE;
704 while(!feof(fd))
705 {
706 if(!fgets(t,BSIZE,fd))break;
707 switch(itype)
708 {
709 case GABEDIT_ORBTYPE_ALPHA :
710 pdest = strstr( t, titlesOrb[itype]);
711 break;
712 case GABEDIT_ORBTYPE_BETA :
713 pdest = strstr( t, titlesOrb[itype]);
714 break;
715 case GABEDIT_ORBTYPE_RESTRICTED:
716 pdest = strstr( t, titlesOrb[itype] );
717 break;
718 case GABEDIT_ORBTYPE_MCSCF:
719 pdest = strstr( t,titlesOrb[itype]);
720 break;
721 case GABEDIT_ORBTYPE_EIGENVECTORS:
722 pdest = strstr( t, titlesOrb[itype] );
723 {
724 gchar dump1[50];
725 gchar dump2[50];
726 gint k = sscanf(t,"%s %s",dump1,dump2);
727 if(k!=1 || strcmp(dump1,titlesOrb[itype])!=0) pdest=NULL;
728 }
729 break;
730 case GABEDIT_ORBTYPE_BOYS_ALPHA:
731 pdest = strstr( t,titlesOrb[itype]);
732 break;
733 case GABEDIT_ORBTYPE_BOYS_BETA:
734 pdest = strstr( t,titlesOrb[itype]);
735 break;
736 case GABEDIT_ORBTYPE_BOYS:
737 pdest = strstr( t,titlesOrb[itype]);
738 break;
739 case GABEDIT_ORBTYPE_EDMISTON_ALPHA:
740 pdest = strstr( t,titlesOrb[itype]);
741 break;
742 case GABEDIT_ORBTYPE_EDMISTON_BETA:
743 pdest = strstr( t,titlesOrb[itype]);
744 break;
745 case GABEDIT_ORBTYPE_EDMISTON:
746 pdest = strstr( t,titlesOrb[itype]);
747 break;
748 case GABEDIT_ORBTYPE_PIPEK_ALPHA:
749 pdest = strstr( t,titlesOrb[itype]);
750 break;
751 case GABEDIT_ORBTYPE_PIPEK_BETA:
752 pdest = strstr( t,titlesOrb[itype]);
753 break;
754 case GABEDIT_ORBTYPE_PIPEK:
755 pdest = strstr( t,titlesOrb[itype]);
756 break;
757 }
758 if ( pdest != NULL )
759 {
760 numorb++;
761 OK = TRUE;
762 break;
763 }
764 }
765 if(!OK && (numorb == 1) )
766 {
767 if(
768 itype==GABEDIT_ORBTYPE_BETA ||
769 itype==GABEDIT_ORBTYPE_EIGENVECTORS ||
770 itype==GABEDIT_ORBTYPE_BOYS_BETA ||
771 itype==GABEDIT_ORBTYPE_EDMISTON_BETA ||
772 itype==GABEDIT_ORBTYPE_PIPEK_BETA
773 )
774 {
775 gchar buffer[BSIZE];
776 sprintf(buffer,_("Sorry, I can not read orbitals from '%s' file. Remove \"print none\" from your input file!\n"),fileName);
777 Message(buffer,_("Error"),TRUE);
778 }
779 FreeTable2(CoefOrbitals,NOrb);
780 g_free(EnerOrbitals);
781 g_free(SymOrbitals);
782 fclose(fd);
783 g_free(t);
784 return FALSE;
785 }
786 if(!OK)
787 {
788 /* Debug("End of read \n");*/
789 fclose(fd);
790 g_free(t);
791 g_free(dum);
792
793 switch(itype)
794 {
795 case GABEDIT_ORBTYPE_ALPHA :
796 case GABEDIT_ORBTYPE_BOYS_ALPHA:
797 case GABEDIT_ORBTYPE_EDMISTON_ALPHA:
798 case GABEDIT_ORBTYPE_PIPEK_ALPHA:
799 CoefAlphaOrbitals = CoefOrbitals;
800 EnerAlphaOrbitals = EnerOrbitals;
801 SymAlphaOrbitals = SymOrbitals;
802 OccAlphaOrbitals = OccOrbitals;
803 NAlphaOcc = 0;
804 for(i=0;i<NOrb;i++) if(OccAlphaOrbitals[i]>0) NAlphaOcc++;
805 NAlphaOrb = NOrb;
806 break;
807
808 case GABEDIT_ORBTYPE_BETA :
809 case GABEDIT_ORBTYPE_BOYS_BETA:
810 case GABEDIT_ORBTYPE_EDMISTON_BETA:
811 case GABEDIT_ORBTYPE_PIPEK_BETA:
812 CoefBetaOrbitals = CoefOrbitals;
813 EnerBetaOrbitals = EnerOrbitals;
814 SymBetaOrbitals = SymOrbitals;
815 OccBetaOrbitals = OccOrbitals;
816 NBetaOcc = 0;
817 for(i=0;i<NOrb;i++) if(OccBetaOrbitals[i]>0) NBetaOcc++;
818 NBetaOrb = NOrb;
819 break;
820
821 case GABEDIT_ORBTYPE_RESTRICTED:
822 case GABEDIT_ORBTYPE_MCSCF:
823 case GABEDIT_ORBTYPE_EIGENVECTORS:
824 case GABEDIT_ORBTYPE_BOYS:
825 case GABEDIT_ORBTYPE_EDMISTON:
826 case GABEDIT_ORBTYPE_PIPEK:
827 CoefAlphaOrbitals = CoefOrbitals;
828 EnerAlphaOrbitals = EnerOrbitals;
829 SymAlphaOrbitals = SymOrbitals;
830 OccAlphaOrbitals = OccOrbitals;
831 NAlphaOcc = 0;
832 for(i=0;i<NOrb;i++) if(OccAlphaOrbitals[i]>0) NAlphaOcc++;
833 for(i=0;i<NOrb;i++) OccAlphaOrbitals[i]/=2;
834 NAlphaOrb = NOrb;
835
836 CoefBetaOrbitals = CoefOrbitals;
837 EnerBetaOrbitals = EnerOrbitals;
838 OccBetaOrbitals = OccAlphaOrbitals;
839 SymBetaOrbitals = SymOrbitals;
840 NBetaOcc = NAlphaOcc;
841 NBetaOrb = NOrb;
842 break;
843 }
844 return TRUE;
845 }
846 switch(itype)
847 {
848 case GABEDIT_ORBTYPE_BOYS_ALPHA:
849 case GABEDIT_ORBTYPE_BOYS_BETA:
850 while(!feof(fd))
851 {
852 if(!fgets(t,BSIZE,fd))break;
853 if(strstr(t,titlesOrb[GABEDIT_ORBTYPE_BOYS]))
854 {
855 { char* e = fgets(t,BSIZE,fd);}
856 break;
857 }
858 }
859 break;
860 case GABEDIT_ORBTYPE_EDMISTON_ALPHA:
861 case GABEDIT_ORBTYPE_EDMISTON_BETA:
862 while(!feof(fd))
863 {
864 if(!fgets(t,BSIZE,fd))break;
865 if(strstr(t,titlesOrb[GABEDIT_ORBTYPE_EDMISTON]))
866 {
867 { char* e = fgets(t,BSIZE,fd);}
868 break;
869 }
870 }
871 break;
872 case GABEDIT_ORBTYPE_PIPEK_ALPHA:
873 case GABEDIT_ORBTYPE_PIPEK_BETA:
874 while(!feof(fd))
875 {
876 if(!fgets(t,BSIZE,fd))break;
877 if(strstr(t,titlesOrb[GABEDIT_ORBTYPE_PIPEK]))
878 {
879 { char* e = fgets(t,BSIZE,fd);}
880 break;
881 }
882 }
883 break;
884 case GABEDIT_ORBTYPE_ALPHA :
885 case GABEDIT_ORBTYPE_BETA :
886 case GABEDIT_ORBTYPE_RESTRICTED:
887 { char* e = fgets(t,BSIZE,fd);}
888 { char* e = fgets(t,BSIZE,fd);}
889 break;
890 case GABEDIT_ORBTYPE_MCSCF:
891 case GABEDIT_ORBTYPE_EIGENVECTORS:
892 break;
893 case GABEDIT_ORBTYPE_BOYS:
894 case GABEDIT_ORBTYPE_EDMISTON:
895 case GABEDIT_ORBTYPE_PIPEK:
896 { char* e = fgets(t,BSIZE,fd);}
897 }
898 /*
899 printf("NOrb=%d\n",NOrb);
900 printf("NAOrb=%d\n",NAOrb);
901 */
902 for(j=0;j<NOrb;j++)
903 {
904 EnerOrbitals[j]=0;
905 OccOrbitals[j]=0;
906 SymOrbitals[j] = g_strdup("DELETED");
907 for(i=0;i<NOrb;i++) CoefOrbitals[j][i]=0.0;
908 }
909 while(fgets(t,BSIZE,fd)&&strstr(t,"Vector"))
910 {
911 gchar* tmp = t;
912 /* printf("%s",t);*/
913 tmp = strstr(t,"Vector")+strlen("Vector");
914 for(i=0;i<strlen(t);i++) if(t[i] =='D') t[i] = 'e';
915 i = atoi(tmp)-1;
916 OccOrbitals[i] = atof(strstr(t,"Occ=")+strlen("Occ="));
917 EnerOrbitals[i] = atof(strstr(t,"E=")+strlen("E="));
918 if(strstr(t,"Symmetry=")) sscanf(strstr(t,"Symmetry=")+strlen("Symmetry="),"%s",SymOrbitals[i]);
919 else SymOrbitals[i] = g_strdup("UNK");
920 if(!fgets(t,BSIZE,fd))break;/* MO Center= */
921 if(!fgets(t,BSIZE,fd))break; /* Bfn. Coefficient */
922 if(!fgets(t,BSIZE,fd))break; /* ------------ */
923 while(fgets(t,BSIZE,fd)&&!this_is_a_backspace(t))
924 {
925 gint j;
926 gdouble c;
927 gint k;
928 /* printf("%s",t);*/
929 k = sscanf(t,"%d %lf", &j, &c);
930 if(k==2 && j>=0 && j<=NOrb) CoefOrbitals[i][j-1]=c;
931
932 if(strlen(t)>44)
933 {
934 k = sscanf(t+40,"%d %lf", &j, &c);
935 if(k==2 && j>=0 && j<=NOrb) CoefOrbitals[i][j-1]=c;
936 }
937 }
938 }
939 /* Debug("End ncart\n");*/
940
941 }while(!feof(fd));
942
943 /* Debug("End of read \n");*/
944 fclose(fd);
945 g_free(t);
946 g_free(dum);
947
948 CoefAlphaOrbitals = CoefOrbitals;
949 EnerAlphaOrbitals = EnerOrbitals;
950 return TRUE;
951 }
952 /********************************************************************************/
read_nwchem_orbitals(gchar * FileName)953 void read_nwchem_orbitals(gchar* FileName)
954 {
955 gint typefile;
956 gchar *t = NULL;
957 gint nrs = 0;
958 gchar** strbasis=NULL;
959 gint i;
960 gboolean Ok;
961 GabEditOrbLocalType typeLocal;
962 gint typebasis = -1;
963 /* gint j,jj; */
964 gulong lineg = 0;
965
966 typefile =get_type_file_orb(FileName);
967 if(typefile==GABEDIT_TYPEFILE_UNKNOWN) return;
968
969
970 if(typefile != GABEDIT_TYPEFILE_NWCHEM)
971 {
972 gchar buffer[BSIZE];
973 sprintf(buffer,_("Sorry, I can not read this format from '%s' file\n"),FileName);
974 Message(buffer,_("Error"),TRUE);
975 return ;
976 }
977
978 free_data_all();
979 t = get_name_file(FileName);
980 set_status_label_info(_("File name"),t);
981 g_free(t);
982 set_status_label_info(_("File type"),"NWChem");
983 set_status_label_info(_("Mol. Orb."),"Reading");
984
985 free_orbitals();
986
987 typebasis = get_type_basis_in_nwchem_file(FileName);
988 if(typebasis == -1)
989 {
990 gchar buffer[BSIZE];
991 sprintf(buffer,
992 _(
993 "Sorry, Gabedit does not support mixed spherical and contaminant cartezian basis functions\n\n"
994 )
995 );
996 Message(buffer,_("Error"),TRUE);
997 set_status_label_info(_("File name"),_("Nothing"));
998 set_status_label_info(_("File type"),_("Nothing"));
999 set_status_label_info(_("Mol. Orb."),_("Nothing"));
1000 return;
1001 }
1002
1003 lineg = read_geomorb_nwchem_file_geom(FileName);
1004 if(lineg == 0)
1005 {
1006 free_geometry();
1007 set_status_label_info(_("File name"),_("Nothing"));
1008 set_status_label_info(_("File type"),_("Nothing"));
1009 set_status_label_info(_("Mol. Orb."),_("Nothing"));
1010 return;
1011 }
1012 strbasis=read_basis_from_a_nwchem_output_file(FileName, &nrs,lineg);
1013 /* printf("nrs=%d\n",nrs);*/
1014 if(strbasis==NULL)
1015 {
1016 if(GeomOrb)
1017 {
1018 init_atomic_orbitals();
1019 for(i=0;i<nCenters;i++) GeomOrb[i].Prop = prop_atom_get("H");
1020 free_geometry();
1021 }
1022 set_status_label_info(_("File name"),_("Nothing"));
1023 set_status_label_info(_("File type"),_("Nothing"));
1024 set_status_label_info(_("Mol. Orb."),_("Nothing"));
1025 return;
1026 }
1027 /*
1028 {
1029 gint i;
1030 for(i=0;i<nrs;i++) printf("%s\n",strbasis[i]);
1031 }
1032 */
1033
1034 set_status_label_info(_("Mol. Orb."),"Reading");
1035 InitializeAll();
1036 if(!DefineNWChemBasisType(strbasis,nrs))
1037 {
1038 gchar buffer[BSIZE];
1039 sprintf(buffer,_("Sorry, I can not read basis from '%s' file\n"),FileName);
1040 Message(buffer,_("Error"),TRUE);
1041 set_status_label_info(_("File name"),_("Nothing"));
1042 set_status_label_info(_("File type"),_("Nothing"));
1043 set_status_label_info(_("Mol. Orb."),_("Nothing"));
1044 return;
1045 }
1046 for(i=0;i<Ntype;i++)
1047 if(Type[i].Ao == NULL)
1048 {
1049 gchar buffer[BSIZE];
1050 sprintf(buffer,_("Sorry, I can not read '%s' file, problem with basis set \n"),FileName);
1051 Message(buffer,_("Error"),TRUE);
1052 return;
1053 }
1054
1055 /*
1056 for(i=0;i<Ntype;i++)
1057 for(j=0;j<Type[i].Norb;j++)
1058 {
1059 printf(" L = %d\n",Type[i].Ao[j].L);
1060 for(jj=0;jj<Type[i].Ao[j].N;jj++)
1061 {
1062 printf(" %lf ",Type[i].Ao[j].Ex[jj]);
1063 printf(" %lf ",Type[i].Ao[j].Coef[jj]);
1064 }
1065 printf("\n");
1066 }
1067 */
1068 DefineType();
1069 buildBondsOrb();
1070 RebuildGeom = TRUE;
1071 reset_grid_limits();
1072 init_atomic_orbitals();
1073 set_status_label_info(_("Geometry"),_("Ok"));
1074 glarea_rafresh(GLArea); /* for geometry*/
1075
1076
1077 if(typebasis == 1)
1078 {
1079 DefineNWChemSphericalBasis();
1080 sphericalBasis = TRUE;
1081 }
1082 else
1083 {
1084 DefineNWChemCartBasis();
1085 sphericalBasis = FALSE;
1086 }
1087
1088 /*PrintAllBasis();*/
1089 NormaliseAllBasis();
1090 /*PrintAllBasis();*/
1091 DefineNOccs();
1092
1093 typeLocal = get_local_orbital_type(FileName);
1094 if(typeLocal!=GABEDIT_ORBLOCALTYPE_UNKNOWN)
1095 {
1096 if(typeLocal==GABEDIT_ORBLOCALTYPE_BOYS)
1097 {
1098 Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_BOYS_ALPHA,lineg);
1099 if(Ok) Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_BOYS_BETA,lineg);
1100 else Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_BOYS,lineg);
1101
1102 }
1103 else
1104 if(typeLocal==GABEDIT_ORBLOCALTYPE_EDMISTON)
1105 {
1106 Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_EDMISTON_ALPHA,lineg);
1107 if(Ok) Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_EDMISTON_BETA,lineg);
1108 else Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_EDMISTON,lineg);
1109 }
1110 else
1111 {
1112 Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_PIPEK_ALPHA,lineg);
1113 if(Ok) Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_PIPEK_BETA,lineg);
1114 else Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_PIPEK,lineg);
1115
1116 }
1117
1118 }
1119 else
1120 {
1121 Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_ALPHA,lineg);
1122 if(Ok)
1123 {
1124 Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_BETA,lineg);
1125 }
1126 else
1127 {
1128 if(!Ok) Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_RESTRICTED,lineg);
1129
1130 if(!Ok) Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_MCSCF,lineg);
1131 if(!Ok) Ok = read_last_orbitals_in_nwchem_file(FileName,GABEDIT_ORBTYPE_EIGENVECTORS,lineg);
1132 }
1133 }
1134
1135 if(Ok)
1136 {
1137 /* PrintAllOrb(CoefAlphaOrbitals);*/
1138 set_status_label_info(_("Mol. Orb."),_("Ok"));
1139 glarea_rafresh(GLArea); /* for geometry*/
1140 NumSelOrb = NAlphaOcc-1;
1141 create_list_orbitals();
1142 }
1143 else
1144 {
1145 free_orbitals();
1146 set_status_label_info(_("File name"),_("Nothing"));
1147 set_status_label_info(_("File type"),_("Nothing"));
1148 set_status_label_info(_("Mol. Orb."),_("Nothing"));
1149 }
1150
1151 }
1152