1 /* UtilsCIF.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 <math.h>
22 #include <stdlib.h>
23 #include <stdio.h>
24 #include <ctype.h>
25 #include "../Common/Global.h"
26 #include "../Utils/Constants.h"
27 #include "../Utils/Utils.h"
28 #include "../Crystallography/Crystallo.h"
29
30 /*************************************************************************************/
g_list_free_all(GList * list,GDestroyNotify free_func)31 static void g_list_free_all (GList * list, GDestroyNotify free_func)
32 {
33 g_list_foreach (list, (GFunc) free_func, NULL);
34 g_list_free (list);
35 }
36 /********************************************************************************/
isItNumber(gchar c)37 static gboolean isItNumber(gchar c)
38 {
39 gchar numb[] = {'0','1','2','3','4','5','6','7','8','9'};
40 gint nums = sizeof(numb)/sizeof(gchar);
41 gint i;
42 for(i=0;i<nums;i++) if(c==numb[i]) return TRUE;
43 return FALSE;
44 }
45 /********************************************************************************/
getStrNum(gchar * str,gint n)46 static gchar* getStrNum(gchar* str, gint n)
47 {
48 gint i;
49 gint len;
50 gchar bl=' ';
51 gchar a='\'';
52 gchar endl='\n';
53 gchar c;
54 gint j;
55 gchar* s = NULL;
56 gint ibegin, iend;
57 if(!str) return s;
58 len = strlen(str);
59 //fprintf(stderr,"str=%s len=%d\n",str,len);
60 ibegin=0;
61 iend=0;
62 c=bl;
63 j=0;
64 i=0;
65 while(str[i]==bl)i++;
66 for(;i<len;i++)
67 {
68 if(str[i]==a) c=a;
69 if(str[i]==c || str[i]==endl || i==len-1)
70 {
71 j++;
72 iend=i;
73 //fprintf(stderr,"ib=%d ie=%d j =%d n=%d\n",ibegin,iend,j,n);
74 if(j==n) break;
75 while(str[i]==bl)i++;
76 ibegin=i;
77 //fprintf(stderr,"new ib=%d\n",ibegin);
78 if(str[i]==a) c=a;
79 if(ibegin==len-1)
80 {
81 j++;
82 iend=ibegin;
83 }
84 }
85 }
86 //fprintf(stderr,"ibegin=%d iend=%d j =%d n=%d\n",ibegin,iend,j,n);
87 if(j!=n) return s;
88 s = g_malloc(sizeof(gchar*)*(iend-ibegin+2));
89 for(i=ibegin; i<=iend;i++) s[i-ibegin] = str[i];
90 s[iend-ibegin+1] = '\0';
91 len = strlen(s);
92 for(i=0;i<len;i++) if(s[i]==a) for(j=i;j<len;j++) s[j] = s[j+1];
93 //fprintf(stderr,"s=%s\n",s);
94 return s;
95 }
96 /********************************************************************************/
beginWithUnderScore(gchar * str)97 static gboolean beginWithUnderScore(gchar* str)
98 {
99 gint i;
100 gint len;
101 if(!str) return FALSE;
102 len = strlen(str);
103 for(i=0;i<len;i++)
104 {
105 if(str[i]==' ' ) continue;
106 if(str[i]=='_' ) return TRUE;
107 else return FALSE;
108 }
109 return FALSE;
110 }
111 /********************************************************************************/
read_onval(FILE * file,gchar * label,gdouble * pVal)112 static gboolean read_onval(FILE* file, gchar* label, gdouble* pVal)
113 {
114 gchar t[BSIZE];
115 gboolean ok = FALSE;
116 rewind(file);
117 //fprintf(stderr,"label=%s\n",label);
118 while(!feof(file))
119 {
120 if(!fgets(t,BSIZE,file)) break;
121 if(beginWithUnderScore(t)) lowercase(t);
122 if(strstr(t,label)) {
123 gchar* val = getStrNum(t, 2);
124 if(val){
125 *pVal = atof(val);
126 ok = TRUE;
127 g_free(val);
128 }
129 break;
130 }
131 }
132 return ok;
133 }
134 /********************************************************************************/
read_cell(Crystal * crystal,FILE * file)135 static gboolean read_cell(Crystal* crystal, FILE* file)
136 {
137 if(! (read_onval(file, "_cell_length_a", &crystal->a) || read_onval(file, "_cell.length_a", &crystal->a))) return FALSE;
138 if(! (read_onval(file, "_cell_length_b", &crystal->b) || read_onval(file, "_cell.length_b", &crystal->b))) return FALSE;
139 if(! (read_onval(file, "_cell_length_c", &crystal->c) || read_onval(file, "_cell.length_c", &crystal->c))) return FALSE;
140 if(! (read_onval(file, "_cell_angle_alpha", &crystal->alpha) || read_onval(file, "_cell.angle_alpha", &crystal->alpha))) return FALSE;
141 if(! (read_onval(file, "_cell_angle_beta", &crystal->beta) || read_onval(file, "_cell.angle_beta", &crystal->beta))) return FALSE;
142 if(! (read_onval(file, "_cell_angle_gamma", &crystal->gamma) || read_onval(file, "_cell.angle_gamma", &crystal->gamma))) return FALSE;
143 /*
144 fprintf(stderr,"a=%f\n",crystal->a);
145 fprintf(stderr,"b=%f\n",crystal->b);
146 fprintf(stderr,"c=%f\n",crystal->c);
147 fprintf(stderr,"alpha=%f\n",crystal->alpha);
148 fprintf(stderr,"beta=%f\n",crystal->beta);
149 fprintf(stderr,"gamma=%f\n",crystal->gamma);
150 */
151 return TRUE;
152 }
153 /****************************************************************************************/
read_loop_block(FILE * file,gchar * tag)154 static GList* read_loop_block(FILE* file, gchar* tag)
155 {
156 gchar t[BSIZE];
157 gboolean ok=FALSE;
158 gboolean beginData=FALSE;
159 GList * list = NULL;
160 if(!goToStr(file,"loop_") ) return NULL;
161 while(!feof(file))
162 {
163 if(!fgets(t,BSIZE,file)) break;
164 if(beginWithUnderScore(t)) lowercase(t);
165 if(strlen(t)<1 && ok) break;
166 if(t[0]=='#' && ok) break;
167 if(t[0]=='\n' && ok) break;
168 if(beginWithUnderScore(t) && ok && beginData) break;
169 if(strstr(t,"loop_") && ok) break;
170 if(strstr(t,"loop_"))
171 {
172 g_list_free_all(list, g_free);
173 list = NULL;
174 beginData=FALSE;
175 }
176 else
177 {
178 gint i;
179 gint len = strlen(t);
180 gchar* data = NULL;
181 for(i=0;i<len;i++) if(t[i]!=' ' && t[i]!='\n' && t[i]!='\r') break;
182 if(i==len && ok) break;
183 for(i=len-1;i>0;i--) if(t[i]=='\n') {t[i]='\0'; break;}
184 data = g_strdup(t);
185 list=g_list_append(list, (gpointer) data);
186 if(strstr(t,tag)) ok = TRUE;
187 if(!beginWithUnderScore(t)) beginData = TRUE;
188 }
189 }
190 if(!ok)
191 {
192 g_list_free_all(list, g_free);
193 list = NULL;
194 }
195 return list;
196 }
197 /****************************************************************************************/
print_block(GList * list)198 static void print_block(GList* list)
199 {
200 gchar t[BSIZE];
201 gboolean ok=FALSE;
202 GList * l = NULL;
203
204 for(l = g_list_first(list); l != NULL; l = l->next)
205 {
206 gchar* line = (gchar*)l->data;
207 fprintf(stderr,"%s\n",line);
208 }
209 }
210 /****************************************************************************************/
buildAtomsListFromBlockGeom(Crystal * crystal,GList * listGeom,gchar * typex,gchar * typey,gchar * typez)211 static gchar* buildAtomsListFromBlockGeom(Crystal* crystal, GList* listGeom, gchar* typex, gchar* typey, gchar* typez)
212 {
213 gint numCol[4];
214 GList * l = NULL;
215 gint i=0;
216 gint numOcc = -1;
217 gboolean ok = FALSE;
218 gint nOccDiffOne = 0;
219 for(i=0;i<4;i++) numCol[i] = -1;
220 i = 0;
221 for(l = g_list_first(listGeom); l != NULL; l = l->next)
222 {
223 gchar* line = (gchar*)l->data;
224 if(beginWithUnderScore(line))
225 {
226 i++;
227 if(strstr(line, typex)) numCol[1] = i;
228 if(strstr(line, typey)) numCol[2] = i;
229 if(strstr(line, typez)) numCol[3] = i;
230 if(strstr(line, "_atom_site_label") && numCol[0]<0) numCol[0] = i;
231 if(strstr(line, "_atom_site.label") && numCol[0]<0) numCol[0] = i;
232 if(strstr(line, "_atom_site.label_atom")) numCol[0] = i;
233 if(strstr(line, "_atom_site_type_symbol")) numCol[0] = i;
234 if(strstr(line, "_atom_site.type_symbol")) numCol[0] = i;
235 if(strstr(line, "_atom_site_occupancy")) numOcc = i;
236 }
237 }
238 //fprintf(stderr,"numCol0 = %d\n",numCol[0]);
239
240 ok = TRUE;
241 for(i=0;i<4;i++) if(numCol[i]<0) { ok = FALSE; break; }
242 if(!ok)
243 {
244 gchar* error = g_strdup_printf("Error in Geom block : no %s or no _atom_site_label\n",typex);
245 fprintf(stderr,"%s\n",error);
246 return error;
247 }
248 if(crystal->atoms) g_list_free_all (crystal->atoms, crystalloFreeAtom);
249 crystal->atoms = NULL;
250
251 ok = TRUE;
252 for(l = g_list_first(listGeom); l != NULL; l = l->next)
253 {
254 gint len;
255 gint k;
256 gchar* line = (gchar*)l->data;
257 CrystalloAtom* cifAtom = NULL;
258 gchar* val = NULL;
259 if(beginWithUnderScore(line)) continue;
260 //fprintf(stderr,"=====> %s\n",line);
261 cifAtom = g_malloc(sizeof(CrystalloAtom));
262 /* symbol */
263 val = getStrNum(line, numCol[0]);
264 //fprintf(stderr,"=====> val=%s\n",val);
265 if(!val) { ok = FALSE; break;}
266 len = strlen(val);
267 for(i=0;i<len;i++)
268 {
269 if(isItNumber(val[i])) {val[i]='\0';break;}
270 if(val[i]=='\0') break;
271 }
272 for(i=len-1;i>=1;i--)
273 if(val[i]==' ') {val[i]='\0';break;}
274
275 len = strlen(val);
276
277 k=-1;
278 for(i=0;i<=len;i++)
279 {
280 if(val[i]==' ') continue;
281 //fprintf(stderr,"=====> val2=%s\n",val);
282 k++;
283 cifAtom->symbol[k]=val[i];
284 }
285
286 crystalloInitAtom(cifAtom, cifAtom->symbol);
287 g_free(val);
288 for(i=1;i<=3;i++)
289 {
290 val = getStrNum(line, numCol[i]);
291 //fprintf(stderr,"=====> val=%s numCol=%d\n",val, numCol[i]);
292 if(!val) { ok = FALSE; break;}
293 cifAtom->C[i-1]= atof(val);
294 g_free(val);
295 }
296 if(!ok) break;
297 if(numOcc>-1)
298 {
299 val = getStrNum(line, numOcc);
300 //fprintf(stderr,"=====> val=%s\n",val);
301 if(val && fabs(atof(val)-1)>1e-6) nOccDiffOne++;
302 }
303 crystal->atoms=g_list_append(crystal->atoms, (gpointer) cifAtom);
304 }
305 if(!ok)
306 {
307 gchar* error = g_strdup_printf("Error during read x, y, z and symbol");
308 fprintf(stderr,"%s\n",error);
309 return error;
310 }
311 if(nOccDiffOne>0)
312 {
313 gchar* error = g_strdup_printf("Error : gabedit cannot yet read CIF file with fractional occupancy");
314 fprintf(stderr,"%s\n",error);
315 return error;
316 }
317 return NULL;
318 }
319 /********************************************************************************/
buildSymOperatorsFromBlockSym(Crystal * crystal,GList * listSymOp)320 static gchar* buildSymOperatorsFromBlockSym(Crystal* crystal, GList* listSymOp)
321 {
322
323 gint numCol = -1;
324 GList * l = NULL;
325 gint i=0;
326 gboolean ok = FALSE;
327 i = 0;
328 for(l = g_list_first(listSymOp); l != NULL; l = l->next)
329 {
330 gchar* line = (gchar*)l->data;
331 if(beginWithUnderScore(line))
332 {
333 i++;
334 if(strstr(line, "_symmetry_equiv_pos_as_xyz")) { numCol = i; break;}
335 }
336 }
337
338 ok = TRUE;
339 if(numCol<0)
340 {
341 gchar* error = g_strdup_printf("Error in sym block : no _symmetry_equiv_pos_as_xyz");
342 fprintf(stderr,"%s\n",error);
343 ok = FALSE;
344 return error;
345 }
346 if(crystal->operators) g_list_free_all (crystal->operators, crystalloFreeSymOp);
347 crystal->operators = NULL;
348
349 ok = TRUE;
350 for(l = g_list_first(listSymOp); l != NULL; l = l->next)
351 {
352
353 gint len;
354 gint i,j;
355 gchar* line = (gchar*)l->data;
356 CrystalloSymOp* cifSymOp = NULL;
357 gchar* val = NULL;
358 if(beginWithUnderScore(line)) continue;
359 //fprintf(stderr,"=====> %s\n",line);
360 cifSymOp = g_malloc(sizeof(CrystalloSymOp));
361 val = getStrNum(line, numCol);
362 //fprintf(stderr,"=====> val=%s\n",val);
363 if(!val) { ok = FALSE; break;}
364 len = strlen(val);
365 for(i=0;i<len;i++) if(val[i]==' ')
366 {
367 for(j=i;j<len-1;j++) val[j] = val[j+1];
368 val[len-1] = '\0';
369 }
370 len = strlen(val);
371 for(i=0;i<len;i++) if(val[i]==',') val[i]=' ';
372 cifSymOp->S[0] = getStrNum(val, 1);
373 cifSymOp->S[1] = getStrNum(val, 2);
374 cifSymOp->S[2] = getStrNum(val, 3);
375 /*
376 fprintf(stderr,"=====> val=%s\n",val);
377 fprintf(stderr,"=====> S0=%s\n",cifSymOp->S[0]);
378 fprintf(stderr,"=====> S1=%s\n",cifSymOp->S[1]);
379 fprintf(stderr,"=====> S2=%s\n",cifSymOp->S[2]);
380 */
381 g_free(val);
382 if(!cifSymOp->S[0] || !cifSymOp->S[1]|| !cifSymOp->S[2]){ g_free(cifSymOp); ok = FALSE; break;}
383 for(i=0;i<3;i++) cifSymOp->w[i] = 0;
384 for(i=0;i<3;i++) cifSymOp->W[i][i] = 0;
385 for(i=0;i<3;i++) for(j=i+1;j<3;j++) cifSymOp->W[i][j] = cifSymOp->W[j][j] = 0;
386 crystalloBuildWwFromStr(cifSymOp);
387 crystal->operators=g_list_append(crystal->operators, (gpointer) cifSymOp);
388 }
389 if(!ok)
390 {
391 gchar* error = g_strdup_printf("Error : problem with symmetry operators");
392 fprintf(stderr,"%s\n",error);
393 return error;
394 }
395 return NULL;
396 }
397 /********************************************************************************/
read_geom_fract_from_cif_file(Crystal * crystal,FILE * file)398 static gchar* read_geom_fract_from_cif_file(Crystal* crystal, FILE* file)
399 {
400 gchar* error = NULL;
401
402 GList* listStrGeom = read_loop_block(file, "_atom_site_fract_x");
403 if(!listStrGeom)
404 {
405 error = g_strdup_printf("Error in geometry block : _atom_site_fract_x");
406 return error;
407 }
408 //print_block(listStrGeom);
409 error = buildAtomsListFromBlockGeom(crystal, listStrGeom,"_atom_site_fract_x","_atom_site_fract_y","_atom_site_fract_z");
410 g_list_free_all(listStrGeom, g_free);
411 if(!error)
412 {
413 crystalloAddTvectorsToGeom(crystal);
414 //crystalloPrintAtoms(crystal->atoms);
415 }
416 return error;
417 }
418 /********************************************************************************/
readOneTransfMatrixElemen(FILE * file,gdouble * e,gint i,gint j)419 static gchar* readOneTransfMatrixElemen(FILE* file, gdouble* e, gint i, gint j)
420 {
421 gchar buffer[BSIZE];
422 sprintf(buffer,"_atom_sites_fract_tran_matrix_%d%d",i+1,j+1);
423 if(read_onval(file, buffer, e)) return NULL;
424 sprintf(buffer,"_atom_sites.fract_transf_matrix[%d][%d]",i+1,j+1);
425 if(read_onval(file, buffer, e)) return NULL;
426 return g_strdup_printf("Error I cannot read %s",buffer);
427 }
428 /********************************************************************************/
readOneTransfVectorElemen(FILE * file,gdouble * e,gint i)429 static gchar* readOneTransfVectorElemen(FILE* file, gdouble* e, gint i)
430 {
431 gchar buffer[BSIZE];
432 sprintf(buffer,"_atom_sites_fract_tran_vector_%d",i+1);
433 if(read_onval(file, buffer, e)) return NULL;
434 sprintf(buffer,"_atom_sites.fract_transf_vector[%d]",i+1);
435 if(read_onval(file, buffer, e)) return NULL;
436 return g_strdup_printf("Error I cannot read %s",buffer);
437 }
438 /********************************************************************************/
readTransfMatrix(FILE * file,gdouble W[][3],gdouble w[])439 static gchar* readTransfMatrix(FILE* file, gdouble W[][3], gdouble w[])
440 {
441 gchar* error = NULL;
442 gint i,j;
443 for(i=0;i<3;i++)
444 for(j=0;j<3;j++)
445 {
446 error = readOneTransfMatrixElemen(file,&W[i][j],i,j);
447 if(error) return error;
448 }
449 /* if vector not defined, I set it to 0 */
450 for(i=0;i<3;i++)
451 {
452 error = readOneTransfVectorElemen(file,&w[i],i);
453 if(error) { g_free(error), error = NULL;}
454 }
455 return NULL;
456 }
457 /********************************************************************************/
read_geom_cartn_from_cif_file(Crystal * crystal,FILE * file)458 static gchar* read_geom_cartn_from_cif_file(Crystal* crystal, FILE* file)
459 {
460 gchar* error = NULL;
461 gdouble W[3][3];
462 gdouble w[] = {0,0,0};
463
464 GList* listStrGeom = read_loop_block(file, "_atom_site.cartn_x");
465 if(!listStrGeom)
466 {
467 error = g_strdup_printf("Error in geom block : no _atom_site.cartn_x");
468 return error;
469 }
470 //print_block(listStrGeom);
471 error = readTransfMatrix(file,W,w);
472 if(error) return error;
473
474 error = buildAtomsListFromBlockGeom(crystal, listStrGeom,"_atom_site.cartn_x","_atom_site.cartn_y","_atom_site.cartn_z");
475 if(error) return error;
476
477 crystalloCartnToFractWw(crystal->atoms,W,w);
478
479 g_list_free_all(listStrGeom, g_free);
480 if(!error)
481 {
482 crystalloAddTvectorsToGeom(crystal);
483 //crystalloPrintAtoms(crystal->atoms);
484 }
485 return error;
486 }
487 /****************************************************************************************/
read_and_apply_symop_from_cif_file(Crystal * crystal,FILE * file)488 static gchar* read_and_apply_symop_from_cif_file(Crystal* crystal, FILE* file)
489 {
490 gchar* error = NULL;
491 GList* listStrSymOp = read_loop_block(file, "_symmetry_equiv_pos_as_xyz");
492 if(!listStrSymOp)
493 {
494 error = g_strdup_printf("Error in sym block : no _symmetry_equiv_pos_as_xyz");
495 return error;
496 }
497 //print_block(listStrSymOp);
498 error = buildSymOperatorsFromBlockSym(crystal, listStrSymOp);
499 g_list_free_all(listStrSymOp, g_free);
500 return error;
501 }
502 /****************************************************************************************/
read_geometry_from_cif_file(Crystal * crystal,gchar * fileName,gboolean applySymOp)503 gchar* read_geometry_from_cif_file(Crystal* crystal, gchar* fileName, gboolean applySymOp)
504 {
505 FILE* file;
506 gchar* error = NULL;
507 initCrystal(crystal);
508 if(!fileName)
509 {
510 error = g_strdup_printf("I cannot open %s file", fileName);
511 return error;
512 }
513 file = FOpen(fileName, "rb");
514 if(!read_cell(crystal, file))
515 {
516
517 error = g_strdup_printf(" Error : I cannot read Cell info\n");
518 fprintf(stderr,"%s\n",error);
519 return error;
520 }
521 error = read_geom_fract_from_cif_file(crystal, file);
522 if(error && strstr(error,"occupancy")) return error;
523 if(error)
524 {
525 gchar* error2 = read_geom_cartn_from_cif_file(crystal, file);
526 if(error2)
527 {
528 gchar* e = g_strdup_printf("%s\n%s",error,error2);
529 g_free(error);
530 g_free(error2);
531 error = e;
532 return error;
533 }
534 else {
535 g_free(error);
536 error = NULL;
537 }
538 }
539 if(applySymOp)
540 {
541 error = read_and_apply_symop_from_cif_file(crystal, file);
542 if(error)
543 {
544 fprintf(stderr,"%s\n",error);
545 g_free(error);
546 error= NULL;
547 }
548 }
549 if(!error)
550 {
551 //crystalloPrintSymOp(crystal->operators);
552 if(applySymOp) crystalloApplySymOperators(&crystal->atoms,crystal->operators);
553 crystalloSetAtomsInBox(crystal->atoms);
554 crystalloRemoveAtomsWithSmallDistance(&crystal->atoms);
555 crystalloFractToCartn(crystal->atoms);
556 //crystalloPrintAtoms(crystal->atoms);
557 }
558 //fprintf(stderr,"End reading\n");
559 return error;
560 }
561 /****************************************************************************************/
562