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