1
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 "../Utils/Transformation.h"
29 #include "../Crystallography/Crystallo.h"
30 #include "../Crystallography/GabeditSPG.h"
31
32 static void setTv(GList* atoms, gdouble Tv[][3]);
33 /*************************************************************************************/
g_list_free_all(GList * list,GDestroyNotify free_func)34 static void g_list_free_all (GList * list, GDestroyNotify free_func)
35 {
36 g_list_foreach (list, (GFunc) free_func, NULL);
37 g_list_free (list);
38 }
39 /********************************************************************************/
detMatrixInt3D(gint mat[][3])40 static gint detMatrixInt3D(gint mat[][3])
41 {
42 gint d = 0;
43 gint i,j;
44 for(i=0;i<3;i++)
45 {
46 d+= mat[i][0]*(mat[(i+1)%3][1]*mat[(i+2)%3][2]-mat[(i+2)%3][1]*mat[(i+1)%3][2]);
47 }
48 return d;
49 }
50 /********************************************************************************/
transposeMatrix3D(gdouble mat[][3])51 static void transposeMatrix3D(gdouble mat[][3])
52 {
53 gdouble t;
54 gint i,j;
55 for(i=0;i<3;i++) for(j=i+1;j<3;j++)
56 {
57 t = mat[i][j];
58 mat[i][j] = mat[j][i];
59 mat[j][i] = t;
60 }
61 }
62 /********************************************************************************/
CInverse3(gdouble invmat[][3],gdouble mat[][3])63 static void CInverse3(gdouble invmat[][3], gdouble mat[][3])
64 {
65 gdouble t4,t6,t8,t10,t12,t14,t17;
66
67 t4 = mat[0][0]*mat[1][1];
68 t6 = mat[0][0]*mat[1][2];
69 t8 = mat[0][1]*mat[1][0];
70 t10 = mat[0][2]*mat[1][0];
71 t12 = mat[0][1]*mat[2][0];
72 t14 = mat[0][2]*mat[2][0];
73 t17 = 1/(t4*mat[2][2]-t6*mat[2][1]-t8*mat[2][2]+t10*mat[2][1]+t12*mat[1][2]-t14*mat
74 [1][1]);
75 invmat[0][0] = (mat[1][1]*mat[2][2]-mat[1][2]*mat[2][1])*t17;
76 invmat[0][1] = -(mat[0][1]*mat[2][2]-mat[0][2]*mat[2][1])*t17;
77 invmat[0][2] = -(-mat[0][1]*mat[1][2]+mat[0][2]*mat[1][1])*t17;
78 invmat[1][0] = -(mat[1][0]*mat[2][2]-mat[1][2]*mat[2][0])*t17;
79 invmat[1][1] = (mat[0][0]*mat[2][2]-t14)*t17;
80 invmat[1][2] = -(t6-t10)*t17;
81 invmat[2][0] = -(-mat[1][0]*mat[2][1]+mat[1][1]*mat[2][0])*t17;
82 invmat[2][1] = -(mat[0][0]*mat[2][1]-t12)*t17;
83 invmat[2][2] = (t4-t8)*t17;
84
85 }
86 /********************************************************************************/
prodMatrix3D(gdouble prod[][3],gdouble A[][3],gdouble B[][3])87 static void prodMatrix3D(gdouble prod[][3], gdouble A[][3], gdouble B[][3])
88 {
89 gint i,j,k;
90 for(i=0;i<3;i++)
91 for(j=0;j<3;j++)
92 {
93 prod[i][j] = 0;
94 for(k=0;k<3;k++) prod[i][j] += A[i][k]*B[k][j];
95 }
96 }
97 /********************************************************************************/
crystalloRotate(GList * atoms,gdouble T[][3],gboolean invers)98 gboolean crystalloRotate(GList* atoms, gdouble T[][3], gboolean invers)
99 {
100 gdouble invA[3][3];
101 GList *l = NULL;
102 gdouble C[] = {0,0,0};
103 gint i,j;
104
105 if(!atoms) return FALSE;
106
107 if(invers) CInverse3(invA, T);
108 else for(i=0;i<3;i++) for(j=0;j<3;j++) invA[i][j] = T[i][j];
109
110 /*
111 {
112 gint i,j;
113 fprintf(stderr,"%s :\n","Inv matrix ");
114 for(i=0;i<3;i++)
115 {
116 for(j=0;j<3;j++) fprintf(stderr,"%f ",invA[i][j]);
117 fprintf(stderr,"\n");
118 }
119 }
120 */
121
122 for(l = g_list_first(atoms); l != NULL; l = l->next)
123 {
124 CrystalloAtom* a = (CrystalloAtom*)l->data;
125 // On All if(!(strstr(a->symbol,"Tv") || strstr(a->symbol,"Tv")))
126 {
127 for(j=0;j<3;j++) C[j] = a->C[j];
128
129 for(j=0;j<3;j++) a->C[j] = 0;
130
131 for(j=0;j<3;j++)
132 for(i=0;i<3;i++) a->C[j] += invA[j][i]*C[i];
133 }
134 }
135 return TRUE;
136 }
137 /*************************************************************************************/
crystalloFreeSymOp(gpointer data)138 void crystalloFreeSymOp(gpointer data)
139 {
140 CrystalloSymOp* symOp = (CrystalloSymOp*)data;
141 gint c;
142 if(!symOp) return;
143 for(c=0;c<3;c++) if(symOp->S[c]) g_free(symOp->S[c]);
144 g_free(symOp);
145 }
146 /********************************************************************************/
crystalloFreeAtom(gpointer data)147 void crystalloFreeAtom(gpointer data)
148 {
149 CrystalloAtom* crystalloAtom= (CrystalloAtom*) data;
150 if(!crystalloAtom) return;
151 g_free(crystalloAtom);
152 }
153 /********************************************************************************/
initCrystal(Crystal * crystal)154 void initCrystal(Crystal* crystal)
155 {
156 if(!crystal) return;
157 crystal->atoms = NULL;
158 crystal->operators = NULL;
159 crystal->alpha = 0;
160 crystal->beta = 0;
161 crystal->gamma = 0;
162 crystal->a = 0;
163 crystal->b = 0;
164 crystal->c = 0;
165 }
166 /********************************************************************************/
freeCrystal(Crystal * crystal)167 void freeCrystal(Crystal* crystal)
168 {
169 if(!crystal) return;
170 if(crystal->atoms) g_list_free_all(crystal->atoms, crystalloFreeAtom);
171 if(crystal->operators) g_list_free_all(crystal->operators, crystalloFreeSymOp);
172 /* freeCrystallo(crystal); */
173 }
174 /****************************************************************************************/
crystalloGetDistance2(CrystalloAtom * a1,CrystalloAtom * a2)175 gdouble crystalloGetDistance2(CrystalloAtom* a1, CrystalloAtom* a2)
176 {
177 return (a1->C[0]-a2->C[0])*(a1->C[0]-a2->C[0])
178 + (a1->C[1]-a2->C[1])*(a1->C[1]-a2->C[1])
179 + (a1->C[2]-a2->C[2])*(a1->C[2]-a2->C[2]);
180 }
181 /****************************************************************************************/
crystalloSmallDistance(CrystalloAtom * a1,CrystalloAtom * a2)182 gboolean crystalloSmallDistance(CrystalloAtom* a1, CrystalloAtom* a2)
183 {
184 static gdouble precision = 1e-4;
185 if(crystalloGetDistance2(a1,a2)<precision) return TRUE;
186 return FALSE;
187 }
188 /****************************************************************************************/
crystalloPrintAtoms(GList * atoms)189 void crystalloPrintAtoms(GList* atoms)
190 {
191 GList * l = NULL;
192
193 for(l = g_list_first(atoms); l != NULL; l = l->next)
194 {
195 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
196 if(!crystalloAtom) break;
197 fprintf(stderr,"%s %f %f %f\n",crystalloAtom->symbol, crystalloAtom->C[0],crystalloAtom->C[1], crystalloAtom->C[2]);
198 }
199 }
200 /****************************************************************************************/
crystalloPrintSymOp(GList * operators)201 void crystalloPrintSymOp(GList* operators)
202 {
203 GList * l = NULL;
204
205 for(l = g_list_first(operators); l != NULL; l = l->next)
206 {
207 CrystalloSymOp* crystalloSymOp = (CrystalloSymOp*)l->data;
208 if(crystalloSymOp)
209 {
210 gint i,j;
211 fprintf(stderr,"%s %s %s\n", crystalloSymOp->S[0], crystalloSymOp->S[1],crystalloSymOp->S[2]);
212 for(i=0;i<3;i++)
213 {
214 for(j=0;j<3;j++) fprintf(stderr,"%f ", crystalloSymOp->W[i][j]);
215 fprintf(stderr,"%f\n", crystalloSymOp->w[i]);
216 }
217 }
218 else break;
219 }
220 }
221 /********************************************************************************/
getSignFromStr(gchar * str,gchar xyz)222 static gint getSignFromStr(gchar* str, gchar xyz)
223 {
224 gchar tmp[10];
225 sprintf(tmp,"-%c",xyz);
226 if(strstr(str,tmp)) return -1;
227 sprintf(tmp,"+%c",xyz);
228 if(strstr(str,tmp)) return 1;
229 sprintf(tmp,"%c",xyz);
230 if(strstr(str,tmp)) return 1;
231 return 0;
232 }
233 /********************************************************************************/
isItNumber(gchar c)234 static gboolean isItNumber(gchar c)
235 {
236 gchar numb[] = {'0','1','2','3','4','5','6','7','8','9'};
237 gint nums = sizeof(numb)/sizeof(gchar);
238 gint i;
239 for(i=0;i<nums;i++) if(c==numb[i]) return TRUE;
240 return FALSE;
241 }
242 /********************************************************************************/
getValwFromStr(gchar * str)243 static gdouble getValwFromStr(gchar* str)
244 {
245 gchar tmp[10];
246 gchar* t;
247 gdouble d;
248 gchar* n;
249 t = strstr(str,"/");
250 if(!t) return 0.0;
251 if(strlen(str)==strlen(t)) return 0.0;
252 d = atof(t+1);
253 //fprintf(stderr,"%f\n",d);
254 if(d<=0) return 0.0;
255 for(n=t-1;n>=str;n--)
256 {
257 if(!isItNumber(n[0]))
258 {
259 if(n[0]=='-') return atof(n)/d;
260 return atof(n+1)/d;
261 }
262 }
263 if(isItNumber(str[0]))return atof(str)/d;
264 return 0;
265 }
266 /********************************************************************************/
crystalloBuildWwFromStr(CrystalloSymOp * crystalloSymOp)267 void crystalloBuildWwFromStr(CrystalloSymOp* crystalloSymOp)
268 {
269 gchar xyz[] = {'x','y','z'};
270 gint i,j;
271 for(i=0;i<3;i++) for(j=0;j<3;j++) crystalloSymOp->W[i][j]=getSignFromStr(crystalloSymOp->S[i], xyz[j]);
272 for(i=0;i<3;i++) crystalloSymOp->w[i] = getValwFromStr(crystalloSymOp->S[i]);
273
274 }
275 /********************************************************************************/
crystalloSetAtomsInBox(GList * atoms)276 gboolean crystalloSetAtomsInBox(GList* atoms)
277 {
278 gdouble precision = 1e-2;
279 GList *l = NULL;
280 if(!atoms) return FALSE;
281
282 for(l = g_list_first(atoms); l != NULL; l = l->next)
283 {
284 gint i;
285 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
286 if(strstr(crystalloAtom->symbol,"Tv")) continue;
287 for(i=0;i<3;i++) while(crystalloAtom->C[i]<0) crystalloAtom->C[i]+=1;
288 for(i=0;i<3;i++) while(crystalloAtom->C[i]>1) crystalloAtom->C[i]-=1;
289 for(i=0;i<3;i++)
290 if(fabs(crystalloAtom->C[i]-1)<precision) crystalloAtom->C[i] = 0;
291 }
292 return TRUE;
293 }
294 /********************************************************************************/
crystalloSetCartnAtomsInBox(GList * atoms)295 gboolean crystalloSetCartnAtomsInBox(GList* atoms)
296 {
297 gboolean ok = crystalloCartnToFract(atoms);
298 if(ok) ok=crystalloSetAtomsInBox(atoms);
299 if(ok) ok=crystalloFractToCartn(atoms);
300 return ok;
301 }
302 /********************************************************************************/
setAtomTypes(CrystalloAtom * crystalAtom,gchar * type)303 static void setAtomTypes(CrystalloAtom* crystalAtom, gchar* type)
304 {
305 sprintf(crystalAtom->mmType,"%s",type);
306 sprintf(crystalAtom->pdbType,"%s",type);
307 sprintf(crystalAtom->residueName,"%s",type);
308 }
309 /********************************************************************************/
copyAtomTypes(CrystalloAtom * a1,CrystalloAtom * a2)310 static void copyAtomTypes(CrystalloAtom* a1, CrystalloAtom* a2)
311 {
312 sprintf(a1->mmType,"%s",a2->mmType);
313 sprintf(a1->pdbType,"%s",a2->pdbType);
314 sprintf(a1->residueName,"%s",a2->residueName);
315 }
316 /********************************************************************************/
crystalloInitAtom(CrystalloAtom * a,gchar * symbol)317 void crystalloInitAtom(CrystalloAtom* a, gchar* symbol)
318 {
319 sprintf(a->mmType,"%s",symbol);
320 setAtomTypes(a,symbol);
321 a->residueNumber=1;
322 a->charge=0.0;
323 }
324 /********************************************************************************/
copyAtom(CrystalloAtom * a1,CrystalloAtom * a2)325 static void copyAtom(CrystalloAtom* a1, CrystalloAtom* a2)
326 {
327 gint i;
328 sprintf(a1->symbol,"%s",a2->symbol);
329 copyAtomTypes(a1,a2);
330 a1->residueNumber = a2->residueNumber;
331 a1->charge = a2->charge;
332 for(i=0;i<3;i++) a1->C[i] = a2->C[i];
333 }
334 /********************************************************************************/
crystalloAddTvectorsToGeom(Crystal * crystal)335 gboolean crystalloAddTvectorsToGeom(Crystal* crystal)
336 {
337 gdouble calpha, cbeta, cgamma, sgamma;
338 gdouble cx,cy;
339 gdouble conv=M_PI/180.0;
340
341 CrystalloAtom* crystalAtom = NULL;
342 if(!crystal) return FALSE;
343 if(!crystal->atoms) return FALSE;
344
345 cbeta=cos(crystal->beta*conv);
346 calpha=cos(crystal->alpha*conv);
347 cgamma=cos(crystal->gamma*conv);
348
349 sgamma=sin(crystal->gamma*conv);
350
351 //fprintf(stderr,"gamma=%f\n",crystal->gamma);
352 //fprintf(stderr,"sgamma=%f\n",sgamma);
353
354 /* a1 = a ex*/
355 crystalAtom = g_malloc(sizeof(CrystalloAtom));
356 sprintf(crystalAtom->symbol,"Tv");
357 setAtomTypes(crystalAtom, crystalAtom->symbol);
358 crystalAtom->residueNumber=1;
359 crystalAtom->charge=0.0;
360
361 crystalAtom->C[0]= crystal->a;
362 crystalAtom->C[1]= 0;
363 crystalAtom->C[2]= 0;
364 crystal->atoms=g_list_append(crystal->atoms, (gpointer) crystalAtom);
365
366 /* a2 = b cos(gamma) ex + b sin(gamma) ey */
367 crystalAtom = g_malloc(sizeof(CrystalloAtom));
368 sprintf(crystalAtom->symbol,"Tv");
369 setAtomTypes(crystalAtom, crystalAtom->symbol);
370 crystalAtom->residueNumber=1;
371 crystalAtom->charge=0.0;
372
373 crystalAtom->C[0]= crystal->b*cgamma;
374 crystalAtom->C[1]= crystal->b*sgamma;
375 crystalAtom->C[2]= 0;
376 crystal->atoms=g_list_append(crystal->atoms, (gpointer) crystalAtom);
377
378 /* a3 = cx ex + cy ey + cz ez */
379 /* cx = c cos(beta) */
380 /* cy = c (cos alpha - cos beta cos gamma) /sin gamma */
381 /* cz = sqrt(c^2 - cx^2 - cy^2) */
382
383 crystalAtom = g_malloc(sizeof(CrystalloAtom));
384 sprintf(crystalAtom->symbol,"Tv");
385 setAtomTypes(crystalAtom, crystalAtom->symbol);
386 crystalAtom->residueNumber=1;
387 crystalAtom->charge=0.0;
388
389 cx = crystal->c*cbeta;
390 cy = crystal->c*(calpha-cbeta*cgamma)/sgamma;
391
392 crystalAtom->C[0]= cx;
393 crystalAtom->C[1]= cy;
394 crystalAtom->C[2]= sqrt(crystal->c*crystal->c-cx*cx-cy*cy);
395 crystal->atoms=g_list_append(crystal->atoms, (gpointer) crystalAtom);
396 return TRUE;
397 }
398 /********************************************************************************/
crystalloPrintNumberOfAtoms(GList * atoms)399 void crystalloPrintNumberOfAtoms(GList* atoms)
400 {
401 gint na=0;
402 GList *la = NULL;
403 for(la = g_list_first(atoms); la != NULL; la = la->next) na++;
404 fprintf(stderr," Number of atoms = %d\n", na);
405 }
406 /********************************************************************************/
crystalloRemoveAtomsWithSmallDistance(GList ** patoms)407 gboolean crystalloRemoveAtomsWithSmallDistance(GList** patoms)
408 {
409 GList *li = NULL;
410 GList *lj = NULL;
411 GList *atoms = *patoms;
412
413 if(!atoms) return FALSE;
414
415 for(li = g_list_first(atoms); li != NULL; li = li->next)
416 {
417 CrystalloAtom* ci = (CrystalloAtom*)li->data;
418 lj = li->next;
419 while (lj != NULL)
420 {
421 GList *next = lj->next;
422 CrystalloAtom* cj = (CrystalloAtom*)lj->data;
423 if(crystalloSmallDistance(ci,cj))
424 {
425 atoms = g_list_delete_link(atoms,lj);
426 crystalloFreeAtom(cj);
427 }
428 lj = next;
429 }
430 }
431 fprintf(stderr," After remove atoms with small distance\n");
432 crystalloPrintNumberOfAtoms(atoms);
433 *patoms = atoms;
434 return TRUE;
435 }
436 /********************************************************************************/
crystalloApplySymOperators(GList ** patoms,GList * operators)437 gboolean crystalloApplySymOperators(GList** patoms, GList* operators)
438 {
439 GList *la = NULL;
440 GList *lo = NULL;
441 GList *atoms = *patoms;
442 GList *endAtom = NULL;
443 gboolean ok = FALSE;
444 gint nOp=0;
445 if(!atoms) return FALSE;
446 if(!operators) return FALSE;
447 for(la = g_list_first(atoms); la != NULL; la = la->next) endAtom = la;
448
449 fprintf(stderr," Befor apply symmetry operators\n");
450 crystalloPrintNumberOfAtoms(atoms);
451
452 for(lo = g_list_first(operators); lo != NULL; lo = lo->next)
453 {
454 CrystalloSymOp* crystalloSymOp = (CrystalloSymOp*)lo->data;
455 for(la = g_list_first(atoms); la != NULL; la = la->next)
456 {
457 gint i;
458 gboolean small = FALSE;
459 CrystalloAtom* crystalloAtom = (CrystalloAtom*)la->data;
460 CrystalloAtom* newCrystalloAtom = NULL;
461 if(!strstr(crystalloAtom->symbol,"Tv"))
462 {
463 newCrystalloAtom = g_malloc(sizeof(CrystalloAtom));
464 copyAtom(newCrystalloAtom, crystalloAtom);
465 for(i=0;i<3;i++)
466 {
467 gint j;
468 newCrystalloAtom->C[i] = 0;
469 for(j=0;j<3;j++) newCrystalloAtom->C[i]+= crystalloSymOp->W[i][j]*crystalloAtom->C[j];
470 newCrystalloAtom->C[i] += crystalloSymOp->w[i];
471 }
472 atoms=g_list_append(atoms, (gpointer) newCrystalloAtom);
473 }
474
475 if(la==endAtom) break;
476 }
477 nOp++;
478 }
479 fprintf(stderr," After apply of %d symmetry operators\n",nOp);
480 crystalloPrintNumberOfAtoms(atoms);
481 *patoms = atoms;
482 return TRUE;
483 }
484 /****************************************************************************************/
crystalloCartnToFractWw(GList * atoms,gdouble W[][3],gdouble w[])485 gboolean crystalloCartnToFractWw(GList* atoms, gdouble W[][3], gdouble w[])
486 {
487 GList *l = NULL;
488 if(!atoms) return FALSE;
489
490 for(l = g_list_first(atoms); l != NULL; l = l->next)
491 {
492 gint i,j;
493 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
494 gdouble C[3];
495 if(strstr(crystalloAtom->symbol,"Tv")) continue;
496 if(strstr(crystalloAtom->symbol,"TV")) continue;
497 for(i=0;i<3;i++) C[i] = 0;
498 for(i=0;i<3;i++) for(j=0;j<3;j++) C[i] += W[i][j]*crystalloAtom->C[j];
499 for(i=0;i<3;i++) C[i] += w[i];
500 for(i=0;i<3;i++) crystalloAtom->C[i] = C[i];
501 }
502 return TRUE;
503 }
504 /********************************************************************************/
crystalloCartnToFract(GList * atoms)505 gboolean crystalloCartnToFract(GList* atoms)
506 {
507 GList *l = NULL;
508 gint nTv = 0;
509 gdouble W[3][3];
510 gdouble Tv[3][3];
511 gdouble w[] = {0,0,0};
512
513 if(!atoms) return FALSE;
514
515 for(l = g_list_first(atoms); l != NULL; l = l->next)
516 {
517 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
518 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV"))
519 {
520 gint i;
521 for(i=0;i<3;i++) Tv[i][nTv] = crystalloAtom->C[i];
522 nTv++;
523 }
524 }
525 if(nTv<3) return FALSE;
526 CInverse3(W,Tv);
527 return crystalloCartnToFractWw(atoms, W, w);
528 }
529 /********************************************************************************/
crystalloFractToCartn(GList * atoms)530 gboolean crystalloFractToCartn(GList* atoms)
531 {
532 GList *l = NULL;
533 gint nTv = 0;
534 gdouble Tv[3][3];
535
536 if(!atoms) return FALSE;
537
538 for(l = g_list_first(atoms); l != NULL; l = l->next)
539 {
540 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
541 if(strstr(crystalloAtom->symbol,"Tv"))
542 {
543 gint i;
544 for(i=0;i<3;i++) Tv[nTv][i] = crystalloAtom->C[i];
545 nTv++;
546 }
547 }
548 if(nTv<3) return FALSE;
549
550 for(l = g_list_first(atoms); l != NULL; l = l->next)
551 {
552 gint i,j;
553 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
554 gdouble C[3];
555 if(strstr(crystalloAtom->symbol,"Tv")) continue;
556 for(i=0;i<3;i++) C[i] = 0;
557 for(i=0;i<3;i++) for(j=0;j<3;j++) C[i] += Tv[j][i]*crystalloAtom->C[j];
558 for(i=0;i<3;i++) crystalloAtom->C[i] = C[i];
559 }
560 return TRUE;
561 }
562 /****************************************************************************************/
crystalloAddReplica(GList ** patoms,gint direction,gint nStep,gboolean scaleTv)563 gboolean crystalloAddReplica(GList** patoms, gint direction, gint nStep, gboolean scaleTv)
564 {
565 GList *l = NULL;
566 GList *lend = NULL;
567 GList *lTv = NULL;
568 gint nTv = 0;
569 gdouble Tv[3][3];
570 GList* atoms = *patoms;
571
572 if(!atoms) return FALSE;
573
574 for(l = g_list_first(atoms); l != NULL; l = l->next)
575 {
576 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
577 //fprintf(stderr,"%s\n",crystalloAtom->symbol);
578 if(strstr(crystalloAtom->symbol,"Tv"))
579 {
580 gint i;
581 for(i=0;i<3;i++) Tv[nTv][i] = crystalloAtom->C[i];
582 if(nTv==direction) lTv = l;
583 if(nTv<=2) nTv++;
584 }
585 lend=l;
586 }
587 if(nTv<direction+1) return FALSE;
588 /*
589 fprintf(stderr,"nTv=%d\n",nTv);
590 fprintf(stderr,"end direction=%d\n",direction);
591 fprintf(stderr,"# atoms=%d\n", crystalloNumberOfAtoms(atoms));
592 */
593
594 for(l = g_list_first(atoms); l != NULL; l = l->next)
595 {
596 gint iBegin=(nStep>0)?1:nStep;
597 gint iEnd = (nStep>0)?nStep-1:-1;
598 gint is;
599 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
600 if(!strstr(crystalloAtom->symbol,"Tv"))
601 for(is=iBegin;is<=iEnd;is++)
602 {
603 gint i;
604 CrystalloAtom* newCrystalloAtom = g_malloc(sizeof(CrystalloAtom));
605 copyAtom(newCrystalloAtom, crystalloAtom);
606 for(i=0;i<3;i++) newCrystalloAtom->C[i] += is*Tv[direction][i];
607 atoms=g_list_append(atoms, (gpointer) newCrystalloAtom);
608 }
609 if(l==lend) break;
610 }
611 if(lTv && scaleTv)
612 {
613 gint i;
614 CrystalloAtom* crystalloAtom = (CrystalloAtom*)lTv->data;
615 for(i=0;i<3;i++) crystalloAtom->C[i] = (nStep)*Tv[direction][i];
616 }
617 *patoms = atoms;
618 /*
619 fprintf(stderr,"nTv=%d\n",nTv);
620 fprintf(stderr,"end direction=%d\n",direction);
621 fprintf(stderr,"# atoms=%d\n", crystalloNumberOfAtoms(atoms));
622 */
623 return TRUE;
624 }
625 /************************************************************************************************************/
buildSuperCellSimple(GList ** patoms,gint nReplicas1,gint nReplicas2,gint nReplicas3)626 gboolean buildSuperCellSimple(GList** patoms, gint nReplicas1, gint nReplicas2, gint nReplicas3)
627 {
628 gboolean ok = TRUE;
629 if(ok && nReplicas1>1) ok = crystalloAddReplica(patoms, 0, nReplicas1,TRUE);
630 if(ok && nReplicas2>1) ok = crystalloAddReplica(patoms, 1, nReplicas2,TRUE);
631 if(ok && nReplicas3>1) ok = crystalloAddReplica(patoms, 2, nReplicas3,TRUE);
632 return ok;
633 }
634 /************************************************************************************************************/
buildSuperCell(GList ** patoms,gint P[][3],gdouble p[])635 gboolean buildSuperCell(GList** patoms, gint P[][3], gdouble p[])
636 {
637 gboolean ok = TRUE;
638 gdouble newTv[3][3];
639 gdouble Tv[3][3];
640 gint nTv = crystalloGetTv(*patoms, Tv);
641 gint i,j,c;
642 gint det=detMatrixInt3D(P);
643 if(det<=0)
644 {
645
646 fprintf(stderr,"The determinant of rotation matrix must be > 0\n");
647 fprintf(stderr,"The determinant of your rotation matrix = %d\n",det);
648 return FALSE;
649 }
650
651 for(i=0;i<nTv;i++)
652 {
653 for(j=0;j<nTv;j++)
654 {
655 if(ok) ok = crystalloAddReplica(patoms, j, P[j][i],FALSE);
656 }
657 }
658 if(ok)
659 for(i=0;i<nTv;i++)
660 {
661 for(c=0;c<3;c++)
662 {
663 newTv[i][c] = 0.0;
664 for(j=0;j<nTv;j++) newTv[i][c] += P[j][i]*Tv[j][c];
665 }
666 }
667 if(ok) setTv(*patoms, newTv);
668 if(ok)
669 {
670 gdouble pCart[3];
671 for(c=0;c<3;c++) pCart[c] = 0;
672 for(c=0;c<3;c++) for(j=0;j<nTv;j++) pCart[c] += p[j]*Tv[j][c];
673 /* The origin is shifted by p = p1 a + p2 b + p3 c
674 The atoms shifted by -p
675 */
676 for(c=0;c<3;c++) pCart[c] = -pCart[c];
677
678 ok = crystalloTranslate(*patoms, pCart, TRUE);
679 }
680 if(ok) crystalloSetCartnAtomsInBox(*patoms);
681 return ok;
682 }
683 /****************************************************************************************/
crystalloNumberOfTv(GList * atoms)684 gint crystalloNumberOfTv(GList* atoms)
685 {
686 GList *l = NULL;
687 gint nTv = 0;
688
689 if(!atoms) return 0;
690
691 for(l = g_list_first(atoms); l != NULL; l = l->next)
692 {
693 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
694 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"Tv")) nTv++;
695 }
696 return nTv;
697 }
698 /****************************************************************************************/
crystalloBuildTablesSymbolsXYZ(GList * atoms,gchar ** atomSymbols[],gdouble * positions[])699 gint crystalloBuildTablesSymbolsXYZ(GList* atoms, gchar** atomSymbols[], gdouble* positions[])
700 {
701 gchar** symbols = NULL;
702 gdouble* X = NULL;
703 gdouble* Y = NULL;
704 gdouble* Z = NULL;
705 GList* l;
706 gint nAtoms = 0;
707 gint i;
708
709 if(!atoms) return nAtoms;
710
711 for(l = g_list_first(atoms); l != NULL; l = l->next) nAtoms++;
712 if(nAtoms<1) return nAtoms;
713
714 symbols = g_malloc(nAtoms*sizeof(gchar*));
715 for(i=0;i<nAtoms;i++) symbols[i] = NULL;
716 X = g_malloc(nAtoms*sizeof(gdouble));
717 Y = g_malloc(nAtoms*sizeof(gdouble));
718 Z = g_malloc(nAtoms*sizeof(gdouble));
719
720 i=0;
721 for(l = g_list_first(atoms); l != NULL; l = l->next)
722 {
723 CrystalloAtom* crystalloAtom;
724 crystalloAtom = (CrystalloAtom*)l->data;
725 symbols[i] = g_strdup(crystalloAtom->symbol);
726 //fprintf(stderr,"symb=%s=\n",symbols[i]);
727 X[i] = crystalloAtom->C[0];
728 Y[i] = crystalloAtom->C[1];
729 Z[i] = crystalloAtom->C[2];
730 i++;
731 }
732 atomSymbols[0] = symbols;
733 positions[0] = X;
734 positions[1] = Y;
735 positions[2] = Z;
736 return nAtoms;
737 }
738 /****************************************************************************************/
crystalloNumberOfAtoms(GList * atoms)739 gint crystalloNumberOfAtoms(GList* atoms)
740 {
741 GList *l = NULL;
742 gint nAtoms = 0;
743
744 if(!atoms) return nAtoms;
745
746 for(l = g_list_first(atoms); l != NULL; l = l->next) nAtoms++;
747
748 return nAtoms;
749 }
750 /****************************************************************************************/
dot(gdouble v1[],gdouble v2[])751 static gdouble dot(gdouble v1[], gdouble v2[])
752 {
753 return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
754 }
angle(gdouble v1[],gdouble v2[])755 static gdouble angle(gdouble v1[], gdouble v2[])
756 {
757 static gdouble prec = 1e-12;
758 gdouble cosa = dot(v1,v2);
759 gdouble d1 = dot(v1,v1);
760 gdouble d2 = dot(v2,v2);
761 if(fabs(d1*d2)<1e-12) return 0.0;
762 cosa /= sqrt(d1*d2);
763 if (cosa>1.0-prec) cosa = 1.0-prec;
764 if (cosa<-1.0+prec) cosa = -1.0+prec;
765 /* fprintf(stderr,"d1=%f d2=%f cosa=%f\n",d1,d2,cosa);*/
766 return acos(cosa)*180.0/M_PI;
767 }
cross(gdouble v1[],gdouble v2[],gdouble cross[])768 static void cross(gdouble v1[], gdouble v2[], gdouble cross[])
769 {
770
771 cross[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]);
772 cross[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]);
773 cross[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]);
774 }
normalize(gdouble v[])775 static void normalize(gdouble v[])
776 {
777 gdouble n = dot(v,v);
778 if(n>0)
779 {
780 n=sqrt(n);
781 v[0] /= n;
782 v[1] /= n;
783 v[2] /= n;
784 }
785 }
786 /****************************************************************************************/
computeTs(gdouble TV[][3],gdouble Ts[][3],gint nTv)787 static void computeTs(gdouble TV[][3], gdouble Ts[][3], gint nTv)
788 {
789 gint i,j;
790 for(i=0;i<3;i++) for(j=0;j<3;j++) Ts[i][j] = ((i==j)?1:0.0);
791 if(nTv>2)
792 {
793 gdouble vol=1;
794 cross(TV[1], TV[2], Ts[0]); // b^c
795 cross(TV[2], TV[0], Ts[1]); // c^a
796 cross(TV[0], TV[1], Ts[2]); // a^b
797 vol=dot(TV[2],Ts[2]); // c.(a^b)
798 for(i=0;i<3;i++) for(j=0;j<3;j++) Ts[i][j] /= vol;
799 }
800 }
801 /********************************************************************************/
crystalloGetVolume(GList * atoms)802 gdouble crystalloGetVolume(GList* atoms)
803 {
804 gdouble Tv1[3] = {0,0,0};
805 gdouble Tv2[3] = {0,0,0};
806 gdouble Tv3[3] = {0,0,0};
807 gdouble V[3];
808 GList *l = NULL;
809 gint nTv = 0;
810
811 if(!atoms) return 0;
812
813 for(l = g_list_first(atoms); l != NULL; l = l->next)
814 {
815 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
816 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"Tv"))
817 {
818 gint j;
819 if(nTv==0) for(j=0;j<3;j++) Tv1[j]=crystalloAtom->C[j];
820 if(nTv==1) for(j=0;j<3;j++) Tv2[j]=crystalloAtom->C[j];
821 if(nTv==2) for(j=0;j<3;j++) Tv3[j]=crystalloAtom->C[j];
822 nTv++;
823 }
824 }
825 cross(Tv1,Tv2,V);
826 return fabs(dot(V,Tv3));
827 }
828 /********************************************************************************/
crystalloGetTv(GList * atoms,gdouble Tv[][3])829 gint crystalloGetTv(GList* atoms, gdouble Tv[][3])
830 {
831 GList *l = NULL;
832 gint nTv = 0;
833
834 if(!atoms) return 0;
835
836 for(l = g_list_first(atoms); l != NULL; l = l->next)
837 {
838 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
839 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"Tv"))
840 {
841 gint j;
842 for(j=0;j<3;j++) Tv[nTv][j]=crystalloAtom->C[j];
843 nTv++;
844 }
845 }
846 return nTv;
847 }
848 /********************************************************************************/
crystalloComputeLengthsAndAngles(Crystal * crystal)849 gboolean crystalloComputeLengthsAndAngles(Crystal* crystal)
850 {
851 gdouble calpha, cbeta, cgamma, sgamma;
852 gdouble cx,cy;
853 gdouble conv=M_PI/180.0;
854 gint nTv;
855 gdouble Tv[3][3];
856
857 if(!crystal) return FALSE;
858 if(!crystal->atoms) return FALSE;
859 nTv = crystalloGetTv(crystal->atoms, Tv);
860 if(nTv!=3) return FALSE;
861
862 /* Length of the basis vectors */
863 crystal->a = sqrt(dot(Tv[0],Tv[0]));
864 crystal->b = sqrt(dot(Tv[1],Tv[1]));
865 crystal->c = sqrt(dot(Tv[2],Tv[2]));
866 crystal->alpha = angle(Tv[1],Tv[2]);
867 crystal->beta = angle(Tv[0],Tv[2]);
868 crystal->gamma = angle(Tv[0],Tv[1]);
869 /* fprintf(stderr,"a=%f\n",crystal->a);*/
870 /* fprintf(stderr,"alpha=%f\n",crystal->alpha);*/
871 return TRUE;
872 }
873 /*****************************************************************************************************************/
getOptimalRadiusForCluster(GList * atoms,gint nSurfaces,gdouble ** surfaces,gdouble * layers)874 static double getOptimalRadiusForCluster(GList* atoms, gint nSurfaces, gdouble** surfaces, gdouble* layers)
875 {
876
877 gdouble Tv[3][3];
878 gint nTv = 0;
879 gint i,j=0;
880 gdouble Ts[3][3];
881
882 nTv = crystalloGetTv(atoms, Tv);
883 computeTs(Tv,Ts,nTv);
884 double radius = -1;
885
886 for(i=0;i<nSurfaces;i++)
887 {
888 gdouble V[3];
889 gint k;
890 gint c;
891 for(k=0;k<3;k++) V[k] = 0;
892 for(k=0;k<3;k++) for(c=0;c<3;c++) V[k] += surfaces[i][c]*Ts[c][k];
893 gdouble len=layers[i]/sqrt(dot(V,V));
894 if(len>radius) radius = len;
895 }
896 radius *=1.5;
897 fprintf(stderr,"Optimal radius for cluster = %f\n",radius);
898 return radius;
899
900
901 }
902 /********************************************************************************/
computeNSteps(GList * atoms,gint nSteps[],gdouble radius)903 static void computeNSteps(GList* atoms, gint nSteps[], gdouble radius)
904 {
905 gdouble Tv[3][3];
906 gint nTv = 0;
907 gint i,j=0;
908 gint k=0;
909
910 nTv = crystalloGetTv(atoms, Tv);
911
912 for(k=0;k<3;k++) nSteps[k]=1;
913 k=0;
914 //printf("Origin = %f %f %f\n", orig[0], orig[1], orig[2]);
915 for(i=0;i<nTv;i++)
916 {
917 gint c;
918 gdouble d= 0;
919 for(c=0;c<3;c++) d += Tv[i][c]*Tv[i][c];
920 d=sqrt(d);
921 nSteps[k]=(gint)(2*radius/d);
922 if(nSteps[k]*d<2*radius) nSteps[k]++;
923 /* R[k]=nSteps[k]*d;*/
924 k++;
925 }
926 fprintf(stderr,"nSteps= ");
927 for(k=0;k<3;k++) fprintf(stderr,"%d ",nSteps[k]);
928 fprintf(stderr,"\n");
929 /*
930 fprintf(stderr,"R=\n");
931 for(k=0;k<3;k++) fprintf(stderr,"%f ",R[k]);
932 fprintf(stderr,"\n");
933 */
934 }
935 /********************************************************************************/
crystalloCreateCellNano(GList ** patoms,gdouble radius)936 void crystalloCreateCellNano(GList** patoms, gdouble radius)
937 {
938 gint nSteps[3];
939 GList* atoms = *patoms;
940 gint k;
941 computeNSteps(atoms, nSteps, radius);
942 buildSuperCellSimple(patoms, nSteps[0], nSteps[1], nSteps[2]);
943 }
944 /*****************************************************************************************************/
crystalloCutPlane(GList ** patoms,gdouble direction[],gdouble layer)945 void crystalloCutPlane(GList** patoms, gdouble direction[], gdouble layer)
946 {
947 GList* atoms = *patoms;
948 GList* l = g_list_first(atoms);
949
950 while (l != NULL)
951 {
952 GList *next = l->next;
953 CrystalloAtom* a = (CrystalloAtom*)l->data;
954 gdouble pscal=dot(direction,a->C);
955 if(pscal>layer)
956 {
957 atoms = g_list_delete_link(atoms,l);
958 crystalloFreeAtom(a);
959 }
960 l = next;
961 }
962 *patoms = atoms;
963 }
964 /*****************************************************************************************************/
crystalloRemoveTv(GList ** patoms)965 void crystalloRemoveTv(GList** patoms)
966 {
967 GList* atoms = *patoms;
968 GList* l = g_list_first(atoms);
969
970 while (l != NULL)
971 {
972 GList *next = l->next;
973 CrystalloAtom* a = (CrystalloAtom*)l->data;
974 if(strstr(a->symbol,"Tv") || strstr(a->symbol,"Tv"))
975 {
976 atoms = g_list_delete_link(atoms,l);
977 crystalloFreeAtom(a);
978 }
979 l = next;
980 }
981 *patoms = atoms;
982 }
983 /********************************************************************************/
crystalloCenter(GList * atoms)984 gint crystalloCenter(GList* atoms)
985 {
986 GList *l = NULL;
987 gint nAtoms = 0;
988 gdouble C[] = {0,0,0};
989 gint j;
990
991 if(!atoms) return 0;
992
993 for(l = g_list_first(atoms); l != NULL; l = l->next)
994 {
995 CrystalloAtom* a = (CrystalloAtom*)l->data;
996 if(!(strstr(a->symbol,"Tv") || strstr(a->symbol,"Tv")))
997 {
998 for(j=0;j<3;j++) C[j] += a->C[j];
999 nAtoms++;
1000 }
1001 }
1002 if(nAtoms<1) return nAtoms;
1003 for(j=0;j<3;j++) C[j] /= nAtoms;
1004 for(l = g_list_first(atoms); l != NULL; l = l->next)
1005 {
1006 CrystalloAtom* a = (CrystalloAtom*)l->data;
1007 if(!(strstr(a->symbol,"Tv") || strstr(a->symbol,"Tv")))
1008 for(j=0;j<3;j++) a->C[j] -= C[j];;
1009 }
1010 return nAtoms;
1011 }
1012 /*********************************************************************************************************/
createWulffCluster(GList ** patoms,gint nSurfaces,gdouble ** surfaces,gdouble * layers)1013 void createWulffCluster(GList** patoms, gint nSurfaces, gdouble** surfaces, gdouble* layers)
1014 {
1015
1016 GList* atoms = *patoms;
1017 gdouble Ts[3][3];
1018 gdouble TvSmall[3][3];
1019 gdouble C[] = {0,0,0};
1020 gint i,j;
1021 gint nSteps[3];
1022 gint nAtoms = 0;
1023 GList* l = NULL;
1024
1025 gint nTvSmall = crystalloGetTv(atoms, TvSmall);
1026
1027
1028 gdouble radius = getOptimalRadiusForCluster(atoms, nSurfaces, surfaces, layers);
1029 computeNSteps(atoms, nSteps, radius);
1030 crystalloCreateCellNano(&atoms, radius);
1031
1032 computeTs(TvSmall,Ts,nTvSmall);
1033
1034 for(l = g_list_first(atoms); l != NULL; l = l->next)
1035 {
1036 CrystalloAtom* a = (CrystalloAtom*)l->data;
1037 if(strstr(a->symbol,"Tv") || strstr(a->symbol,"Tv")) continue;
1038 else {
1039 gint j;
1040 for(j=0;j<3;j++) C[j] += a->C[j];
1041 nAtoms++;
1042 }
1043 }
1044 if(nAtoms<1) return;
1045
1046 for(j=0;j<3;j++) C[j] /= nAtoms;
1047
1048 for(l = g_list_first(atoms); l != NULL; l = l->next)
1049 {
1050 CrystalloAtom* a = (CrystalloAtom*)l->data;
1051 if(strstr(a->symbol,"Tv") || strstr(a->symbol,"Tv")) continue;
1052 else {
1053 gint j;
1054 for(j=0;j<3;j++) a->C[j] -= C[j];
1055 }
1056 }
1057
1058 for(i=0;i<nSurfaces;i++)
1059 {
1060 gdouble V[3];
1061 gint k;
1062 gint c;
1063 for(k=0;k<3;k++) V[k] = 0;
1064 for(k=0;k<3;k++) for(c=0;c<3;c++) V[k] += surfaces[i][c]*Ts[c][k];
1065 crystalloCutPlane(&atoms, V, layers[i]);
1066 for(k=0;k<3;k++) V[k] = -V[k];
1067 crystalloCutPlane(&atoms, V, layers[i]);
1068 }
1069 for(l = g_list_first(atoms); l != NULL; l = l->next)
1070 {
1071 CrystalloAtom* a = (CrystalloAtom*)l->data;
1072 if(strstr(a->symbol,"Tv") || strstr(a->symbol,"Tv")) continue;
1073 else {
1074 gint j;
1075 for(j=0;j<3;j++) a->C[j] += C[j];
1076 }
1077 }
1078 crystalloRemoveTv(&atoms);
1079 crystalloCenter(atoms);
1080 fprintf(stderr,"Number of atoms in cluster.\n");
1081 crystalloPrintNumberOfAtoms(atoms);
1082 *patoms = atoms;
1083 }
1084 /********************************************************************************/
1085 /* Compute the greatest common divisor */
gcd(gint a,gint b)1086 static gint gcd(gint a, gint b)
1087 {
1088 a = abs(a);
1089 b = abs(b);
1090 if (a == 0 || b == 0) return 1;
1091
1092 while (a != b)
1093 {
1094 while (a < b) b -= a;
1095 while (b < a) a -= b;
1096 }
1097 return a;
1098 }
1099 /********************************************************************************/
1100 /* Extended Euclidean algorithm */
ext_gcd(gint * sp,gint * sq,gint a,gint b)1101 static gint ext_gcd(gint* sp, gint *sq, gint a, gint b)
1102 {
1103 a = abs(a);
1104 b = abs(b);
1105 gint aa[2]={1,0};
1106 gint bb[2]={0,1};
1107 gint q;
1108
1109 while( a!=0 && b!=0)
1110 {
1111 q = a / b;
1112 a = a % b;
1113 aa[0] = aa[0] - q*aa[1];
1114 bb[0] = bb[0] - q*bb[1];
1115 if (a == 0)
1116 {
1117 *sp = aa[1];
1118 *sq = bb[1];
1119 return b;
1120 };
1121 q = b / a;
1122 b = b % a;
1123 aa[1] = aa[1] - q*aa[0];
1124 bb[1] = bb[1] - q*bb[0];
1125 if (b == 0)
1126 {
1127 *sp = aa[0];
1128 *sq = bb[0];
1129 return a;
1130 }
1131 }
1132 if(a==0)
1133 {
1134 *sp = aa[1];
1135 *sq = bb[1];
1136 return b;
1137 }
1138 if (b == 0)
1139 {
1140 *sp = aa[0];
1141 *sq = bb[0];
1142 return a;
1143 }
1144 return a;
1145 }
1146 /********************************************************************************/
reduce3Int(gint tab[])1147 static gint reduce3Int(gint tab[])
1148 {
1149 gint g = gcd(tab[0], tab[1]);
1150 g = gcd(g, tab[2]);
1151 tab[0] /= g;
1152 tab[1] /= g;
1153 tab[2] /= g;
1154 return g;
1155 }
1156 /*********************************************************************************************************/
setTv(GList * atoms,gdouble Tv[][3])1157 static void setTv(GList* atoms, gdouble Tv[][3])
1158 {
1159
1160 GList *l = NULL;
1161 gint iTv;
1162 iTv = 0;
1163 l = g_list_first(atoms);
1164 while (l != NULL)
1165 {
1166 GList *next = l->next;
1167 CrystalloAtom* a = (CrystalloAtom*)l->data;
1168 if(strstr(a->symbol,"Tv") || strstr(a->symbol,"TV"))
1169 {
1170 if(iTv<=2)
1171 {
1172 gint j;
1173 for(j=0;j<3;j++) a->C[j] = Tv[iTv][j];
1174 iTv++;
1175 }
1176 if(iTv==3) break;
1177 }
1178 l=next;
1179 }
1180 }
1181 /****************************************************************************************/
getRotMatrixAtoB(gdouble rotMatrix[][3],gdouble A[],gdouble B[])1182 static void getRotMatrixAtoB(gdouble rotMatrix[][3], gdouble A[], gdouble B[])
1183 {
1184 /*
1185 v = cross(A,B);
1186 ssc = [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0];
1187 R = eye(3) + ssc + ssc^2*(1-dot(A,B))/(norm(v))^2;
1188
1189 */
1190 gdouble V1[3];
1191 gdouble V2[3];
1192 gdouble V[3];
1193 gdouble ssc[3][3];
1194 gint i,j;
1195 gdouble normV;
1196 gdouble norm;
1197
1198 for(i=0;i<3;i++) for(j=0;j<3;j++) rotMatrix[i][j] = 0;
1199 for(i=0;i<3;i++) rotMatrix[i][i] = 1;
1200
1201 for(i=0;i<3;i++) V1[i]=A[i];
1202 norm=0;
1203 for(i=0;i<3;i++) norm+=V1[i]*V1[i];
1204 norm=sqrt(norm);
1205 if(norm<=0) return;
1206 for(i=0;i<3;i++) V1[i] /= norm;
1207
1208 for(i=0;i<3;i++) V2[i]=B[i];
1209 norm=0;
1210 for(i=0;i<3;i++) norm+=V2[i]*V2[i];
1211 norm=sqrt(norm);
1212 if(norm<=0) return;
1213 for(i=0;i<3;i++) V2[i] /= norm;
1214
1215 cross(V1,V2,V);
1216 normV=0;
1217 for(i=0;i<3;i++) normV += V[i]*V[i];
1218 if(normV<=0) return;
1219
1220 ssc[0][0] = 0;
1221 ssc[0][1] = -V[2];
1222 ssc[0][2] = V[1];
1223
1224 ssc[1][0] = V[2];
1225 ssc[1][1] = 0;
1226 ssc[1][2] = -V[0];
1227
1228 ssc[2][0] = -V[1];
1229 ssc[2][1] = V[0];
1230 ssc[2][2] = 0;
1231
1232
1233 norm = (1-dot(V1,V2))/normV;
1234 for(i=0;i<3;i++) for(j=0;j<3;j++)
1235 {
1236 gint k;
1237 rotMatrix[i][j] = 0;
1238 for(k=0;k<3;k++) rotMatrix[i][j] += ssc[i][k]*ssc[k][j]*norm;
1239 }
1240 for(i=0;i<3;i++) for(j=0;j<3;j++) rotMatrix[i][j] += ssc[i][j];
1241 for(i=0;i<3;i++) rotMatrix[i][i] += 1;
1242 }
1243 /****************************************************************************************/
getRotMatrix2Toz(gdouble rotMatrix[][3],gdouble Tv[][3])1244 static void getRotMatrix2Toz(gdouble rotMatrix[][3], gdouble Tv[][3])
1245 {
1246 gdouble B[3] = {0,0,1};
1247 gint i = 2;
1248 gdouble A[3] = {Tv[i][0], Tv[i][1], Tv[i][2]};
1249 getRotMatrixAtoB(rotMatrix, A, B);
1250 }
1251 /****************************************************************************************/
getRotMatrix0Tox(gdouble rotMatrix[][3],gdouble Tv[][3])1252 static void getRotMatrix0Tox(gdouble rotMatrix[][3], gdouble Tv[][3])
1253 {
1254 gdouble B[3] = {1,0,0};
1255 gint i = 0;
1256 gdouble A[3] = {Tv[i][0], Tv[i][1], Tv[i][2]};
1257 getRotMatrixAtoB(rotMatrix, A, B);
1258 }
1259 /****************************************************************************************/
getZcutoff(gdouble * surface,gdouble zlayer,gint nTv,gdouble Ts[][3])1260 static gdouble getZcutoff(gdouble* surface, gdouble zlayer, gint nTv, gdouble Ts[][3])
1261 {
1262 gdouble V[3];
1263 gint k;
1264 gint c;
1265 gdouble zCutoff = -1;
1266
1267 for(c=0;c<3;c++) V[c] = 0;
1268 for(c=0;c<3;c++) for(k=0;k<nTv;k++) V[c] += surface[k]*Ts[k][c];
1269 zCutoff = zlayer/sqrt(dot(V,V));
1270 //fprintf(stderr,"zCutoff = %f\n",zCutoff);
1271 return zCutoff;
1272
1273 }
1274 /****************************************************************************************/
getNormalSurface(gdouble normal[],gdouble surface[],gint nTv,gdouble Ts[][3])1275 static void getNormalSurface(gdouble normal[], gdouble surface[], gint nTv, gdouble Ts[][3])
1276 {
1277 gint k;
1278 gint c;
1279
1280 for(c=0;c<3;c++) normal[c] = 0;
1281 for(c=0;c<3;c++) for(k=0;k<nTv;k++) normal[c] += surface[k]*Ts[k][c];
1282 normalize(normal);
1283 }
1284 /*********************************************************************************************************/
createSlab(GList ** patoms,gdouble surface[],gdouble layers[],gdouble emptySpaceSize,gboolean orientSurfaceXY)1285 void createSlab(GList** patoms, gdouble surface[], gdouble layers[], gdouble emptySpaceSize, gboolean orientSurfaceXY)
1286 {
1287
1288 gint hkl[] = {(gint)surface[0], (gint)surface[1], (gint)surface[2]};
1289 gint p,q;
1290 gint g;
1291 gint i;
1292 gint i0,i1,n0;
1293 gint nTv;
1294 gdouble Tv[3][3];
1295 gdouble surfaceVectors[3][3];
1296 gint j;
1297 gdouble V1[3];
1298 gdouble V2[3];
1299 GList* atoms = *patoms;
1300 GList *l = NULL;
1301 gint iTv;
1302 //gdouble Transf[3][3];
1303 gdouble normalSurface[3];
1304 gdouble Ts[3][3];
1305 gdouble zCutoff = 1.0;
1306
1307 gdouble zWithEmptySpace = 1+ fabs(emptySpaceSize);
1308
1309 nTv = crystalloGetTv(atoms, Tv);
1310 if(nTv!=3)
1311 {
1312 fprintf(stderr,"Error : I cannot build a slab with %d translation vectors\n",nTv);
1313 return;
1314 }
1315 computeTs(Tv, Ts, nTv);
1316 zCutoff = getZcutoff(surface, (gint)layers[2], nTv, Ts);
1317 getNormalSurface(normalSurface, surface, nTv, Ts);
1318 for(j=0;j<3;j++) normalSurface[j] *= zCutoff;
1319
1320 n0=0;
1321 i0=0;
1322 i1=1;
1323 for(i=0;i<3;i++)
1324 {
1325 if(hkl[i]==0) {n0++; i0=i;}
1326 else i1 = i;
1327 }
1328 if(n0==3)
1329 {
1330 fprintf(stderr,"Error : Miller indices cannot be all zero\n");
1331 return;
1332 }
1333 if(n0==2)
1334 {
1335 for(j=0;j<3;j++) surfaceVectors[2][j] = Tv[i1][j];
1336 for(j=0;j<3;j++) surfaceVectors[0][j] = Tv[(i1+1)%3][j];
1337 for(j=0;j<3;j++) surfaceVectors[1][j] = Tv[(i1+2)%3][j];
1338
1339 /*
1340 for(i=0;i<3;i++) for(j=0;j<3;j++) Transf[i][j] = 0;
1341 Transf[2][i1] = 1;
1342 Transf[0][(i1+1)%3] = 1;
1343 Transf[1][(i1+2)%3] = 1;
1344 */
1345 }
1346 else
1347 {
1348 /*
1349 p,q = ext_gcd(k,l)
1350 k1 = dot( p*(k*a1-h*a2)+q*(l*a1-h*a3) , l*a2-k*a3)
1351 k2 = dot( l*(k*a1-h*a2)-k*(l*a1-h*a3) , l*a2-k*a3)
1352 if abs(k2)>tol:
1353 c = -int(round(k1/k2))
1354 p,q = p+c*l, q-c*k
1355 v1 = p*array((k,-h,0))+q*array((l,0,-h))
1356 v2 = reduce(array((0,l,-k)))
1357 a,b = ext_gcd(p*k+q*l,h)
1358 v3 = array((b,a*p,a*q))
1359 */
1360 gint ih = i0;
1361 gint ik = (i0+1)%3;
1362 gint il = (i0+2)%3;
1363 gint d=reduce3Int(hkl);
1364 gint p,q;
1365 gint g = ext_gcd(&p, &q, hkl[ik], hkl[il]);
1366 gdouble k1,k2;
1367 gint a,b;
1368 gint c;
1369 gint h = hkl[ih];
1370 gint k = hkl[ik];
1371 gint l = hkl[il];
1372 //fprintf(stderr,"ihkl=%d %d %d\n",ih,ik,il);
1373 //fprintf(stderr,"hkl=%d %d %d\n",h,k,l);
1374 //fprintf(stderr,"pq=%d %d\n",p,q);
1375 for(j=0;j<3;j++) V1[j] = p*(k*Tv[ih][j]-h*Tv[ik][j])+q*(l*Tv[ih][j]-h*Tv[il][j]);
1376 for(j=0;j<3;j++) V2[j] = l*Tv[ik][j]-k*Tv[il][j];
1377 k1=dot(V1,V2);
1378 for(j=0;j<3;j++) V1[j] = l*(k*Tv[ih][j]-h*Tv[ik][j])-k*(l*Tv[ih][j]-h*Tv[il][j]);
1379 k2=dot(V1,V2);
1380 c = 0;
1381 //fprintf(stderr,"k1,k2=%f %f\n",k1,k2);
1382 if(fabs(k2)>1e-8) c = -(gint)(k1/k2);
1383 p += c*l;
1384 q -= c*k;
1385 g = gcd(l,k);
1386 //fprintf(stderr,"k,l,g=%d %d %d\n",k,l,g);
1387 /* v1*/
1388 /*
1389 Transf[0][0] = p*k+q*l;
1390 Transf[0][1] = -p*h;
1391 Transf[0][2] = -q*h;
1392 */
1393 for(j=0;j<3;j++) surfaceVectors[0][j] = p*(k*Tv[ih][j]-h*Tv[ik][j])+q*(l*Tv[ih][j]-h*Tv[il][j]);
1394 /* v2*/
1395 /*
1396 Transf[1][0] = 0;
1397 Transf[1][1] = l/g;
1398 Transf[1][2] = -k/g;
1399 */
1400 for(j=0;j<3;j++) surfaceVectors[1][j] = (l/g*Tv[ik][j]-k/g*Tv[il][j]);
1401 /* v3 */
1402 ext_gcd(&a, &b, p*k+q*l,h);
1403 //fprintf(stderr,"a,b=%d %d\n",a,b);
1404 //fprintf(stderr,"pq=%d %d\n",p,q);
1405 /*
1406 Transf[2][0] = b;
1407 Transf[2][1] = a*p;
1408 Transf[2][2] = a*q;
1409 */
1410 for(j=0;j<3;j++) surfaceVectors[2][j] = b*Tv[ih][j]+a*p*Tv[ik][j]+a*q*Tv[il][j];
1411
1412 }
1413 /*
1414 {
1415 gint i,j;
1416 fprintf(stderr,"%s :\n","Tv");
1417 for(i=0;i<3;i++)
1418 {
1419 for(j=0;j<3;j++) fprintf(stderr,"%f ",Tv[i][j]);
1420 fprintf(stderr,"\n");
1421 }
1422 }
1423
1424 {
1425 gint i,j;
1426 fprintf(stderr,"%s :\n","surfaceVectors");
1427 for(i=0;i<3;i++)
1428 {
1429 for(j=0;j<3;j++) fprintf(stderr,"%f ",surfaceVectors[i][j]);
1430 fprintf(stderr,"\n");
1431 }
1432 }
1433 {
1434 gint i,j;
1435 fprintf(stderr,"%s :\n","Transf matrix : Tv => surfaceVectors");
1436 for(i=0;i<3;i++)
1437 {
1438 for(j=0;j<3;j++) fprintf(stderr,"%f ",Transf[i][j]);
1439 fprintf(stderr,"\n");
1440 }
1441 }
1442 */
1443 {
1444 gint i,j;
1445 for(i=0;i<3;i++)
1446 for(j=0;j<3;j++) Tv[i][j] =surfaceVectors[i][j];
1447 }
1448 //fprintf(stderr,"VolumeTv=%f\n", crystalloGetVolume(atoms));
1449 setTv(atoms, surfaceVectors);
1450 crystalloSetCartnAtomsInBox(atoms);
1451 //fprintf(stderr,"VolumeSV=%f\n", crystalloGetVolume(atoms));
1452 buildSuperCellSimple(&atoms, (gint)layers[0], (gint)layers[1], (gint)layers[2]);
1453 fprintf(stderr,"Volume SuperCell befor third per to 2 others=%f\n", zWithEmptySpace*crystalloGetVolume(atoms));
1454 /* set third vector perpendocular to first and second ones*/
1455 if(orientSurfaceXY)
1456 {
1457 nTv = crystalloGetTv(atoms, Tv);
1458 //fprintf(stderr,"\n"); for(j=0;j<3;j++) fprintf(stderr,"%f ",Tv[2][j]);
1459 //fprintf(stderr,"\n"); for(j=0;j<3;j++) fprintf(stderr,"%f ",normalSurface[j]);
1460 for(j=0;j<3;j++) Tv[2][j] = normalSurface[j];
1461 setTv(atoms, Tv);
1462 nTv = crystalloGetTv(atoms, Tv);
1463 //for(j=0;j<3;j++) fprintf(stderr,"%f ",Tv[2][j]);
1464 crystalloSetCartnAtomsInBox(atoms);
1465 fprintf(stderr,"Volume SuperCell after third per to 2 others=%f\n", crystalloGetVolume(atoms));
1466 }
1467 /*
1468 if(orientSurfaceXY)
1469 {
1470 gdouble T[3];
1471 nTv = crystalloGetTv(atoms, Tv);
1472 //fprintf(stderr,"\n"); for(j=0;j<3;j++) fprintf(stderr,"%f ",Tv[2][j]);
1473 //fprintf(stderr,"\n"); for(j=0;j<3;j++) fprintf(stderr,"%f ",normalSurface[j]);
1474
1475 for(j=0;j<3;j++) T[j] = normalSurface[j]-Tv[2][j];
1476 for(j=0;j<3;j++) Tv[2][j] = normalSurface[j];
1477
1478 setTv(atoms, Tv);
1479 nTv = crystalloGetTv(atoms, Tv);
1480 //for(j=0;j<3;j++) fprintf(stderr,"%f ",Tv[2][j]);
1481 crystalloTranslate(atoms, T,TRUE);
1482 crystalloSetCartnAtomsInBox(atoms);
1483 fprintf(stderr,"Volume SuperCell after third per to 2 others=%f\n", crystalloGetVolume(atoms));
1484 }
1485 */
1486 /* set Tv0 on x and Tv2 on z */
1487 if(orientSurfaceXY)
1488 {
1489 gdouble rotMatrix[3][3];
1490 nTv = crystalloGetTv(atoms, Tv);
1491 getRotMatrix2Toz(rotMatrix, Tv);
1492 //{ gint i,j; fprintf(stderr,"%s :\n","Tv"); for(i=0;i<3;i++) { for(j=0;j<3;j++) fprintf(stderr,"%f ",Tv[i][j]); fprintf(stderr,"\n"); } }
1493 crystalloRotate(atoms, rotMatrix, FALSE);
1494 nTv = crystalloGetTv(atoms, Tv);
1495 //{ gint i,j; fprintf(stderr,"%s :\n","Tv"); for(i=0;i<3;i++) { for(j=0;j<3;j++) fprintf(stderr,"%f ",Tv[i][j]); fprintf(stderr,"\n"); } }
1496 getRotMatrix0Tox(rotMatrix, Tv);
1497 crystalloRotate(atoms, rotMatrix, FALSE);
1498 nTv = crystalloGetTv(atoms, Tv);
1499 //{ gint i,j; fprintf(stderr,"%s :\n","Tv"); for(i=0;i<3;i++) { for(j=0;j<3;j++) fprintf(stderr,"%f ",Tv[i][j]); fprintf(stderr,"\n"); } }
1500 }
1501 /* add empty space */
1502 {
1503 nTv = crystalloGetTv(atoms, Tv);
1504 for(j=0;j<3;j++) Tv[2][j] *= zWithEmptySpace;
1505 setTv(atoms, Tv);
1506 }
1507
1508 *patoms = atoms;
1509 }
1510 /********************************************************************************/
crystalloPrimitiveCell(GList ** patoms,gdouble symprec)1511 gboolean crystalloPrimitiveCell(GList** patoms, gdouble symprec)
1512 {
1513 gboolean ok = FALSE;
1514 GList* newAtoms = NULL;
1515 if(!patoms || !*patoms) return ok;
1516 ok = crystalloCartnToFract(*patoms);
1517 if(ok)
1518 {
1519 ok = FALSE;
1520 newAtoms = crystalloPrimitiveSPG(*patoms, symprec);
1521 if(newAtoms)
1522 {
1523 ok = TRUE;
1524 g_list_free_all(*patoms, crystalloFreeAtom);
1525 *patoms = newAtoms;
1526 }
1527 }
1528 crystalloFractToCartn(*patoms);
1529 return ok;
1530 }
1531 /********************************************************************************/
crystalloTranslate(GList * atoms,gdouble T[],gboolean cartn)1532 gboolean crystalloTranslate(GList* atoms, gdouble T[], gboolean cartn)
1533 {
1534 GList *l = NULL;
1535 gint j;
1536
1537 if(!atoms) return FALSE;
1538
1539 if(!cartn) crystalloFractToCartn(atoms);
1540
1541 for(l = g_list_first(atoms); l != NULL; l = l->next)
1542 {
1543 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
1544 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
1545 else
1546 {
1547 gint j;
1548 for(j=0;j<3;j++) crystalloAtom->C[j] += T[j];
1549 }
1550 }
1551 if(!cartn) crystalloCartnToFract(atoms);
1552 return TRUE;
1553 }
1554 /********************************************************************************/
1555 /* atoms coordinates must be in fract */
crystalloChangeCell(GList ** patoms,gdouble newTv[][3])1556 gboolean crystalloChangeCell(GList** patoms, gdouble newTv[][3])
1557 {
1558 GList *l = NULL;
1559 gdouble Tv[3][3];
1560 gdouble P[3][3];
1561 gdouble C[3][3];
1562 gint i,j;
1563 gint nTv;
1564
1565 if(!patoms || !*patoms) return FALSE;
1566
1567 nTv = crystalloGetTv(*patoms, Tv);
1568 if(nTv!=3) return FALSE;
1569
1570 nTv = 0;
1571 for(l = g_list_first(*patoms); l != NULL; l = l->next)
1572 {
1573 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
1574 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV"))
1575 {
1576 gint j;
1577 for(j=0;j<3;j++) crystalloAtom->C[j] = newTv[nTv][j];
1578 nTv++;
1579 }
1580 }
1581 if(nTv!=3) return FALSE;
1582
1583 //transposeMatrix3D(Tv);
1584 //transposeMatrix3D(newTv);
1585 /*
1586 Original basis vectors (abc)
1587 Final basis vectors (a′b′c′)
1588 The transformation matrix is obtained by P=(abc)(a′b′c′)^−1,
1589 */
1590 CInverse3(C, newTv);
1591 prodMatrix3D(P, Tv,C);
1592 for(l = g_list_first(*patoms); l != NULL; l = l->next)
1593 {
1594 CrystalloAtom* a = (CrystalloAtom*)l->data;
1595 if(strstr(a->symbol,"Tv") || strstr(a->symbol,"TV")) continue;
1596 for(i=0;i<3;i++) C[i][0] = a->C[i];
1597 for(i=0;i<3;i++) a->C[i] = 0;
1598 for(i=0;i<3;i++) for(j=0;j<3;j++) a->C[i] += P[j][i]*C[j][0];
1599 }
1600 crystalloSetAtomsInBox(*patoms);
1601 //transposeMatrix3D(newTv);
1602 return TRUE;
1603 }
1604 /********************************************************************************/
crystalloReduceCell(GList ** patoms,gdouble symprec,GabeditCrystalloReductionType type)1605 gboolean crystalloReduceCell(GList** patoms, gdouble symprec, GabeditCrystalloReductionType type)
1606 {
1607 gboolean ok = FALSE;
1608 gint nTv = 0;
1609 gdouble newTv[3][3];
1610 GList* atoms = *patoms;
1611
1612 if(!atoms) return ok;
1613
1614 if(type == GABEDIT_CRYSTALLO_REDUCTION_PRIMITIVE) return crystalloPrimitiveCell(patoms, symprec);
1615
1616 crystalloCartnToFract(atoms);
1617 if(type == GABEDIT_CRYSTALLO_REDUCTION_NIGGLI) ok = crystalloGetNiggli(atoms, newTv, symprec);
1618 else if(type == GABEDIT_CRYSTALLO_REDUCTION_DELAUNAY) ok = crystalloGetDelaunay(atoms, newTv, symprec);
1619 if(ok) ok = crystalloChangeCell(patoms, newTv);
1620 crystalloFractToCartn(atoms);
1621
1622 return ok;
1623 }
1624 /********************************************************************************/
generatePearsonSymbol(char pearsonSymbol[],char groupName[],gint numGroup)1625 static gboolean generatePearsonSymbol(char pearsonSymbol[], char groupName[], gint numGroup)
1626 {
1627 gchar crystalclass = ' ';
1628 gchar latticetype = ' ';
1629 if(numGroup <= 2) crystalclass = 'a';
1630 else if(numGroup <= 15) crystalclass = 'm';
1631 else if(numGroup <= 74) crystalclass = 'o';
1632 else if(numGroup <= 142) crystalclass = 't';
1633 else if(numGroup <= 194) crystalclass = 'h';
1634 else crystalclass = 'c';
1635 latticetype = groupName[0];
1636 if(latticetype=='A') latticetype = 'C';
1637 sprintf(pearsonSymbol,"%c%c", crystalclass, latticetype);
1638 return TRUE;
1639 }
1640 /********************************************************************************/
1641 /* Hinuma et al. Comput. Mat. Science 128 (2017) 140-184, tables 93 & 94 */
generateExtendedPearsonSymbol(char extendedPearsonSymbol[],Crystal * crystal,gdouble symprec)1642 gboolean generateExtendedPearsonSymbol(char extendedPearsonSymbol[], Crystal* crystal, gdouble symprec)
1643 {
1644 gchar crystalclass = ' ';
1645 gchar latticetype = ' ';
1646 gchar eType = ' ';
1647 gdouble a=crystal->a;
1648 gdouble b=crystal->b;
1649 gdouble c=crystal->c;
1650 gdouble a2=a*a;
1651 gdouble b2=b*b;
1652 gdouble c2=c*c;
1653 char groupName[100];
1654 gboolean ok = crystalloCartnToFract(crystal->atoms);
1655 gint numGroup = crystalloGetGroupName(crystal->atoms, groupName, symprec);
1656 if(ok) ok=crystalloFractToCartn(crystal->atoms);
1657 if(numGroup <= 2) crystalclass = 'a';
1658 else if(numGroup <= 15) crystalclass = 'm';
1659 else if(numGroup <= 74) crystalclass = 'o';
1660 else if(numGroup <= 142) crystalclass = 't';
1661 else if(numGroup <= 194) crystalclass = 'h';
1662 else crystalclass = 'c';
1663 latticetype = groupName[0];
1664 if(latticetype=='B') latticetype = 'C';
1665 eType='1';
1666 if(crystalclass=='t' && latticetype=='I' && c > a) eType='2';
1667 if(crystalclass=='c' && latticetype=='P' && numGroup>=207) eType='2';
1668 if(crystalclass=='c' && latticetype=='F' && numGroup>=207) eType='2';
1669 if(crystalclass=='o' && latticetype=='F' )
1670 {
1671 if (1.0/a2 > (1.0/b2+1.0/c2)) eType='1';
1672 else if(1.0/c2 > (1.0/a2+1.0/b2)) eType='2';
1673 else eType='3';
1674 }
1675 if(crystalclass=='o' && latticetype=='I' )
1676 {
1677 if (c>a && c>b) eType='1';
1678 else if(a>b && a>c) eType='2';
1679 else eType='3';
1680 }
1681 if(crystalclass=='o' && latticetype=='C' && a>b ) eType='2';
1682 if(crystalclass=='o' && latticetype=='A' && b>c ) eType='2';
1683 if(crystalclass=='h' && latticetype=='P' )
1684 {
1685 if (!( (numGroup>=143 && numGroup<=149) || numGroup==151 || numGroup==153 || numGroup==157 || (numGroup>=159 && numGroup<=163))) eType='2';
1686 }
1687 if(crystalclass=='h' && latticetype=='R' && a*sqrt(3.0)>c*sqrt(2.0) ) eType='2';
1688 if(crystalclass=='m' && latticetype=='C' )
1689 {
1690 gdouble sbeta=sin(crystal->beta/180.0*M_PI);
1691 gdouble cbeta=cos(crystal->beta/180.0*M_PI);
1692 /*
1693 gdouble r=a*a*sbeta*sbeta/b/b-a/c*cbeta;
1694 fprintf(stderr,"a=%f b=%f c=%f alpha=%f beta=%f gamma=%f r=%f\n",
1695 a,b,c,crystal->alpha,crystal->beta,crystal->gamma,r);
1696 */
1697 if(b<a*sbeta) eType='1';
1698 else if(a*a*sbeta*sbeta/b/b-a/c*cbeta<1) eType='2';
1699 else eType='3';
1700 }
1701 if(crystalclass=='a' && latticetype=='P' )
1702 {
1703 eType='2';
1704 if(!crystalloAllRecObtuse(crystal->atoms, symprec)) eType='3';
1705 }
1706 sprintf(extendedPearsonSymbol,"%c%c%c", crystalclass, latticetype,eType);
1707 return TRUE;
1708 }
1709 /********************************************************************************/
crystalloGetSpaceSymmetryGroup(GList * atoms,char groupName[],gdouble symprec)1710 gint crystalloGetSpaceSymmetryGroup(GList* atoms, char groupName[], gdouble symprec)
1711 {
1712 gboolean ok = crystalloCartnToFract(atoms);
1713 gint numGroup = crystalloGetGroupName(atoms, groupName, symprec);
1714 if(ok) ok=crystalloFractToCartn(atoms);
1715 return numGroup;
1716 }
1717 /********************************************************************************/
getWyckoffsLetter(gint i)1718 static gchar getWyckoffsLetter(gint i)
1719 {
1720 gchar w[] ={'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s','t','u','v','w','x','y','z'};
1721 gint size = sizeof(w)/sizeof(gchar);
1722 if(i>=0 && i<size) return w[i];
1723 return ' ';
1724 }
1725 /********************************************************************************/
1726 /*
1727 static gchar* getWyckoffsList(gint * wyckoffs, gint nAtoms)
1728 {
1729 gchar* tmp = NULL;
1730 gint i;
1731
1732 if(nAtoms<1) return NULL;
1733 tmp = g_strdup_printf("%c",getWyckoffsLetter(wyckoffs[0]));
1734 for(i=1;i<nAtoms;i++)
1735 {
1736 gchar* t = tmp;
1737 tmp = g_strdup_printf("%s-%c",t,getWyckoffsLetter(wyckoffs[i]));
1738 if(t) g_free(t);
1739 }
1740 return tmp;
1741 }
1742 */
1743 /********************************************************************************/
1744 /*
1745 static gchar* getStrListOfEquivAtoms(gint * equivalent_atoms, gint nAtoms)
1746 {
1747 gchar* tmp = NULL;
1748 gint i;
1749
1750 if(nAtoms<1) return NULL;
1751 tmp = g_strdup_printf("%d",equivalent_atoms[0]);
1752 for(i=1;i<nAtoms;i++)
1753 {
1754 gchar* t = tmp;
1755 tmp = g_strdup_printf("%s-%d",t,equivalent_atoms[i]);
1756 if(t) g_free(t);
1757 }
1758 return tmp;
1759 }
1760 */
1761 /********************************************************************************/
hasAnOldEquivalent(gint * equivalent_atoms,gint nAtoms,gint j)1762 static gboolean hasAnOldEquivalent(gint * equivalent_atoms, gint nAtoms, gint j)
1763 {
1764 gboolean ok=FALSE;
1765 gint i;
1766 if(nAtoms<1) return ok;
1767 for(i=0;i<j;i++)
1768 {
1769 if( equivalent_atoms[j] == equivalent_atoms[i]) { ok = TRUE; break; }
1770 }
1771 return ok;
1772 }
1773 /********************************************************************************/
getOneStrOp(gint v,gint i,gchar res[])1774 static void getOneStrOp(gint v, gint i, gchar res[])
1775 {
1776 gchar a[] ={'x','y','z'};
1777 gint size = sizeof(a)/sizeof(gchar);
1778
1779 if(i<0|| i>2) sprintf(res," ");
1780 else if(v==0) sprintf(res," ");
1781 else if(abs(v)>1000) sprintf(res," ");
1782 else if(v<-1) sprintf(res,"%d%c",v,a[i]);
1783 else if(v>1) sprintf(res,"+%d%c",v,a[i]);
1784 else if(v==1) sprintf(res,"+%c",a[i]);
1785 else if(v==-1) sprintf(res,"-%c",a[i]);
1786 else sprintf(res," ");
1787 }
1788 /********************************************************************************/
removeFirstSpaceAndPlus(gchar * str)1789 static void removeFirstSpaceAndPlus(gchar* str)
1790 {
1791 while(str[0]=='+' || str[0]==' ')
1792 {
1793 gint len=strlen(str);
1794 gint i;
1795 for(i=0;i<len;i++) str[i]=str[i+1];
1796 }
1797 }
1798 /********************************************************************************/
removeAllWhitespace(gchar * str)1799 static void removeAllWhitespace(gchar* str)
1800 {
1801 /* remove spaces */
1802 gint i,j;
1803 for(i=0;i<strlen(str);i++)
1804 {
1805 if(str[i]==' ')
1806 {
1807 gint len=strlen(str);
1808 for(j=i;j<len;j++) str[j]=str[j+1];
1809 i--;
1810 }
1811 }
1812 }
1813 /********************************************************************************/
1814 /*
1815 static gchar* getSymmetryPosxyz(int W[][3], double w[])
1816 {
1817 gchar rot[3][10];
1818 gchar t[10];
1819 gchar* tmp[3] = {NULL,NULL,NULL};
1820 gint i,j;
1821 gchar* res = NULL;
1822 gint fact = 1000;// 32
1823
1824 for(j=0;j<3;j++)
1825 {
1826 gint n,d,g;
1827 for(i=0;i<3;i++)
1828 {
1829 getOneStrOp(W[j][i], i, rot[i]);
1830 }
1831 n=(gint)(w[j]*fact+0.5);
1832 d = fact;
1833 g = gcd(n,d);
1834 n /= g;
1835 d /= g;
1836 if(abs(n)<1) sprintf(t," ");
1837 else if(d==1) sprintf(t," ");
1838 else sprintf(t,"%+d/%d",n,d);
1839 //sprintf(t,"%f",w[j]);
1840 tmp[j] = g_strdup_printf("%s%s%s%s",t,rot[0], rot[1],rot[2]);
1841 removeFirstSpaceAndPlus(tmp[j]);
1842 }
1843 res = g_strdup_printf("%s,%s,%s",tmp[0], tmp[1],tmp[2]);
1844 for(j=0;j<3;j++) if(tmp[j]) g_free(tmp[j]);
1845 removeAllWhitespace(res);
1846 return res;
1847 }
1848 */
1849 /********************************************************************************/
1850 /*
1851 static gchar* getAllSymmetryPosxyz(int W[][3][3], double w[][3], int nOp)
1852 {
1853 gint i;
1854 gchar* tmp = NULL;
1855 if(nOp<1) return NULL;
1856 tmp = g_strdup("loop_\n_symmetry_equiv_pos_as_xyz");
1857 for(i=0;i<nOp;i++)
1858 {
1859 gchar* t = tmp;
1860 gchar* t1 = getSymmetryPosxyz(W[i],w[i]);
1861 tmp = g_strdup_printf("%s\n%s",t,t1);
1862 if(t) g_free(t);
1863 if(t1) g_free(t1);
1864 }
1865 return tmp;
1866 }
1867 */
1868 /********************************************************************************/
crystalloGetVASPAtomsPositions(GList * atoms)1869 gchar* crystalloGetVASPAtomsPositions(GList* atoms)
1870 {
1871 gint i,j;
1872 gchar* tmp = NULL;
1873 GList* l;
1874 gint* types = NULL;
1875 gint nTypes = 0;
1876 gchar** listTypes = NULL;
1877 gint *nListTypes = NULL;
1878 gchar** symbols = NULL;
1879 gchar* t1 = NULL;
1880 gboolean ok = FALSE;
1881 gint nTv = 0;
1882 gint nAtoms = crystalloNumberOfAtoms(atoms)-3;
1883 if(nAtoms<1) fprintf(stderr,"Error nAtoms<1\n");
1884 if(nAtoms<1) return NULL;
1885 nTv = 0;
1886 for(l = g_list_first(atoms); l != NULL; l = l->next)
1887 {
1888 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
1889 if(strstr(crystalloAtom->symbol,"Tv") ||strstr(crystalloAtom->symbol,"TV")) nTv++;
1890 }
1891 if(nTv!=3) fprintf(stderr,"Error nTv!=3\n");
1892 if(nTv!=3) return NULL;
1893
1894 nAtoms += nTv;
1895
1896 types = g_malloc(nAtoms*sizeof(gint));
1897 for(i=0;i<nAtoms;i++) types[i] = -1;
1898
1899 symbols = g_malloc(nAtoms*sizeof(gchar*));
1900 for(i=0;i<nAtoms;i++) symbols[i] = NULL;;
1901
1902 nTypes = 0;
1903 i=-1;
1904 for(l = g_list_first(atoms); l != NULL; l = l->next)
1905 {
1906 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
1907 i++;
1908 symbols[i] = g_strdup(crystalloAtom->symbol);
1909 }
1910
1911
1912 for(i=0;i<nAtoms;i++)
1913 {
1914 gint j;
1915 ok = TRUE;
1916 if(strstr(symbols[i],"Tv")) continue;
1917 if(strstr(symbols[i],"TV")) continue;
1918
1919 for(j=0;j<i;j++)
1920 {
1921 if(!strcmp(symbols[i],symbols[j]))
1922 {
1923 types[i]= types[j];
1924 ok = FALSE;
1925 break;
1926 }
1927 }
1928 if(ok)
1929 {
1930 types[i]= nTypes;
1931 nTypes++;
1932 }
1933 }
1934 if(nTypes<1) fprintf(stderr,"nTypes = %d\n",nTypes);
1935 if(nTypes<1) return NULL;
1936 listTypes = g_malloc(nTypes*sizeof(gchar*));
1937 nListTypes = g_malloc(nTypes*sizeof(gint));
1938 for(i=0;i<nTypes;i++) listTypes[i] = NULL;
1939 for(i=0;i<nTypes;i++) nListTypes[i] = 0;
1940 for(i=0;i<nAtoms;i++)
1941 {
1942 if(types[i]==-1) continue;
1943 nListTypes[types[i]]++;
1944 if(!listTypes[types[i]]) listTypes[types[i]] = g_strdup(symbols[i]);
1945 }
1946 i=-1;
1947 tmp = g_strdup_printf("Coordinates with POSCAR format\n%s\n","1.0");
1948 for(l = g_list_first(atoms); l != NULL; l = l->next)
1949 {
1950 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
1951 if(strstr(crystalloAtom->symbol,"Tv") ||strstr(crystalloAtom->symbol,"TV"))
1952 {
1953 gchar* t1 = tmp;
1954 tmp = g_strdup_printf("%s%20.10lf %20.10lf %20.10lf\n",
1955 t1,
1956 crystalloAtom->C[0], crystalloAtom->C[1], crystalloAtom->C[2]
1957 );
1958 if(t1) g_free(t1);
1959 }
1960 }
1961 t1 = tmp;
1962 tmp = g_strdup_printf("%s%6s ",t1," ");
1963 if(t1) g_free(t1);
1964 for(i=0;i<nTypes;i++)
1965 {
1966 t1 = tmp;
1967 tmp = g_strdup_printf("%s%6s ",t1,listTypes[i]);
1968 if(t1) g_free(t1);
1969 }
1970 t1 = tmp; tmp = g_strdup_printf("%s\n",t1); if(t1) g_free(t1);
1971 t1 = tmp; tmp = g_strdup_printf("%s%6s ",t1," "); if(t1) g_free(t1);
1972 for(i=0;i<nTypes;i++)
1973 {
1974 t1 = tmp;
1975 tmp = g_strdup_printf("%s%6d ",t1,nListTypes[i]);
1976 if(t1) g_free(t1);
1977 }
1978 t1 = tmp; tmp = g_strdup_printf("%s\nDirect\n",t1); if(t1) g_free(t1);
1979
1980 for(j=0;j<nTypes;j++)
1981 {
1982 GList *l = NULL;
1983 i=-1;
1984 for(l = g_list_first(atoms); l != NULL; l = l->next)
1985 {
1986 gchar w=' ';
1987 i++;
1988 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
1989 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
1990 if(types[i]==-1) continue;
1991 if(strcmp(crystalloAtom->symbol,listTypes[j]))continue;
1992 t1 = tmp;
1993 tmp = g_strdup_printf("%s%6s %20.10f %20.10f %20.10f # %s %c\n",t1,
1994 " ",crystalloAtom->C[0], crystalloAtom->C[1],crystalloAtom->C[2],
1995 crystalloAtom->symbol,
1996 w);
1997 if(t1) g_free(t1);
1998 }
1999 }
2000 for(i=0;i<nTypes;i++) if(listTypes[i]) g_free(listTypes[i]);
2001 if(listTypes) g_free(listTypes);
2002 if(symbols) g_free(symbols);
2003 if(nListTypes) g_free(nListTypes);
2004 if(types) g_free(types);
2005 return tmp;
2006 }
2007 /********************************************************************************/
getCIFAtomsPositionsAndWyckoff(GList * atoms,gint * wyckoffs,gint * equivalent_atoms,gint nAtoms,gboolean allAtoms)2008 static gchar* getCIFAtomsPositionsAndWyckoff(GList* atoms, gint * wyckoffs, gint* equivalent_atoms, gint nAtoms, gboolean allAtoms)
2009 {
2010 gint i;
2011 gchar* tmp = NULL;
2012 GList* l;
2013 if(nAtoms<1) return NULL;
2014
2015 tmp = g_strdup_printf("%s",
2016 "loop_\n"
2017 "_atom_site_label\n"
2018 "_atom_site_type_symbol\n"
2019 "_atom_site_fract_x\n"
2020 "_atom_site_fract_y\n"
2021 "_atom_site_fract_z\n"
2022 "_atom_site_Wyckoff symbol\n");
2023
2024 i=-1;
2025 for(l = g_list_first(atoms); l != NULL; l = l->next)
2026 {
2027 gchar* t1;
2028 gchar* t2;
2029 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
2030 if(strstr(crystalloAtom->symbol,"Tv")) continue;
2031 if(strstr(crystalloAtom->symbol,"TV")) continue;
2032 i++;
2033 if(allAtoms || !hasAnOldEquivalent(equivalent_atoms, nAtoms, i))
2034 {
2035 t1 = tmp;
2036 tmp = g_strdup_printf("%s %s %s %lf %lf %lf %c\n",t1,
2037 crystalloAtom->mmType, crystalloAtom->symbol,
2038 crystalloAtom->C[0], crystalloAtom->C[1], crystalloAtom->C[2],
2039 getWyckoffsLetter(wyckoffs[i]));
2040 if(t1) g_free(t1);
2041 }
2042 }
2043 return tmp;
2044 }
2045 /********************************************************************************/
getCIFCell(Crystal * crystal)2046 static gchar* getCIFCell(Crystal* crystal)
2047 {
2048 gchar* tmp = g_strdup_printf(
2049 "_cell_length_a %0.8lf\n"
2050 "_cell_length_b %0.8lf\n"
2051 "_cell_length_c %0.8lf\n"
2052 "_cell_angle_alpha %0.8lf\n"
2053 "_cell_angle_beta %0.8lf\n"
2054 "_cell_angle_gamma %0.8lf\n",
2055 crystal->a, crystal->b, crystal->c,
2056 crystal->alpha, crystal->beta, crystal->gamma
2057 );
2058 return tmp;
2059 }
2060 /********************************************************************************/
getSymmetryStr(int W[][3],double w[],int j)2061 static gchar* getSymmetryStr(int W[][3], double w[], int j)
2062 {
2063 gchar rot[3][10];
2064 gchar t[10];
2065 gchar* tmp = NULL;
2066 gint i;
2067 gint fact = 1000;/* 32 ? */
2068
2069 {
2070 gint n,d,g;
2071 for(i=0;i<3;i++)
2072 {
2073 getOneStrOp(W[j][i], i, rot[i]);
2074 }
2075 n=(gint)(w[j]*fact+0.5);
2076 d = fact;
2077 g = gcd(n,d);
2078 n /= g;
2079 d /= g;
2080 if(abs(n)<1) sprintf(t," ");
2081 else if(d==1) sprintf(t," ");
2082 else sprintf(t,"%+d/%d",n,d);
2083 //sprintf(t,"%f",w[j]);
2084 tmp = g_strdup_printf("%s%s%s%s",t,rot[0], rot[1],rot[2]);
2085 removeFirstSpaceAndPlus(tmp);
2086 }
2087 removeAllWhitespace(tmp);
2088 return tmp;
2089 }
2090 /********************************************************************************/
crystalloBuildSymOperators(Crystal * crystal,int W[][3][3],double w[][3],int nOp)2091 gboolean crystalloBuildSymOperators(Crystal* crystal, int W[][3][3], double w[][3], int nOp)
2092 {
2093 gint numCol = -1;
2094 GList * l = NULL;
2095 gint i,j,k;
2096 if(nOp<1) return FALSE;
2097 if(crystal->operators) g_list_free_all(crystal->operators, crystalloFreeSymOp);
2098 crystal->operators = NULL;
2099 for(k=0;k<nOp;k++)
2100 {
2101 CrystalloSymOp* cifSymOp = g_malloc(sizeof(CrystalloSymOp));
2102 for(j=0;j<3;j++) cifSymOp->S[j] = getSymmetryStr(W[k], w[k], j);
2103 for(j=0;j<3;j++) cifSymOp->w[j] = w[k][j];
2104 for(i=0;i<3;i++) for(j=0;j<3;j++) cifSymOp->W[i][j] = cifSymOp->W[j][j] = W[k][i][j];
2105 crystal->operators=g_list_append(crystal->operators, (gpointer) cifSymOp);
2106 }
2107 return crystal->operators != NULL;
2108 }
2109 /********************************************************************************/
crystalloGetCIFOperators(GList * operators)2110 gchar* crystalloGetCIFOperators(GList* operators)
2111 {
2112 gchar* tmp = NULL;
2113 GList* l;
2114 if(!operators) return NULL;
2115 tmp = g_strdup("loop_\n_symmetry_equiv_pos_as_xyz\n");
2116
2117 for(l = g_list_first(operators); l != NULL; l = l->next)
2118 {
2119 gchar* t;
2120 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
2121 CrystalloSymOp* cifSymOp =(CrystalloSymOp*) l->data;
2122 t = tmp;
2123 tmp = g_strdup_printf("%s%s,%s,%s\n",t, cifSymOp->S[0], cifSymOp->S[1], cifSymOp->S[2]);
2124 if(t) g_free(t);
2125 }
2126 return tmp;
2127 }
2128 /********************************************************************************/
crystalloGetCIF(Crystal * crystal,gdouble symprec,gboolean withSymmetryOperators)2129 gchar* crystalloGetCIF(Crystal* crystal, gdouble symprec, gboolean withSymmetryOperators)
2130 {
2131 gchar* info = NULL;
2132 GList* atoms = crystal->atoms;
2133 gboolean ok = crystalloCartnToFract(atoms);
2134 SpglibDataset * spgDataSet = crystalloGetDataSet(atoms, symprec);
2135 if(!spgDataSet)
2136 info = g_strdup("Error : Sorry I cannnt find the Space group of this system\n");
2137 else{
2138 gchar* atomPositionsStr = NULL;
2139 gchar* cellInfoStr = NULL;
2140 gchar* cifOperators = NULL;
2141 gchar* nSymOps = NULL;
2142 char pearsonSymbol[5];
2143 crystalloComputeLengthsAndAngles(crystal);
2144 cellInfoStr = getCIFCell(crystal);
2145
2146 generatePearsonSymbol(pearsonSymbol, spgDataSet->international_symbol,spgDataSet->spacegroup_number);
2147
2148
2149 if(withSymmetryOperators)
2150 {
2151 atomPositionsStr = getCIFAtomsPositionsAndWyckoff(atoms, spgDataSet->wyckoffs, spgDataSet->equivalent_atoms, spgDataSet->n_atoms,FALSE);
2152 crystalloBuildSymOperators(crystal, spgDataSet->rotations, spgDataSet->translations, spgDataSet->n_operations);
2153 cifOperators = crystalloGetCIFOperators(crystal->operators);
2154 nSymOps = g_strdup_printf("# of symmetry operator = %d\n", spgDataSet->n_operations);
2155 }else
2156 {
2157 atomPositionsStr = getCIFAtomsPositionsAndWyckoff(atoms, spgDataSet->wyckoffs, spgDataSet->equivalent_atoms, spgDataSet->n_atoms,TRUE);
2158 cifOperators = g_strdup("\n");
2159 nSymOps = g_strdup("\n");
2160 }
2161
2162 info =g_strdup_printf(
2163 "# Pearson symbol %s \n"
2164 "# Space group name %s \n"
2165 "_space_group_IT_number\t\t %d\n"
2166 "_symmetry_space_group_name_H-M\t\t%s\n"
2167 "_symmetry_space_group_name_Hall\t\t%s\n"
2168 "%s\n"
2169 "%s\n"
2170 "%s\n"
2171 "%s\n"
2172 ,
2173 pearsonSymbol,
2174 spgDataSet->international_symbol,
2175 spgDataSet->spacegroup_number,
2176 spgDataSet->pointgroup_symbol,
2177 spgDataSet->hall_symbol,
2178 cellInfoStr,
2179 nSymOps,
2180 cifOperators,
2181 atomPositionsStr
2182 );
2183 if(cifOperators) g_free(cifOperators);
2184 if(cellInfoStr) g_free(cellInfoStr);
2185 if(atomPositionsStr) g_free(atomPositionsStr);
2186 if(nSymOps) g_free(nSymOps);
2187 }
2188 if(ok) ok=crystalloFractToCartn(atoms);
2189 /* fprintf(stderr,"%s\n", getStrListOfEquivAtoms(spgDataSet->equivalent_atoms, spgDataSet->n_atoms));*/
2190 spg_free_dataset(spgDataSet);
2191
2192 return info;
2193 }
2194 /********************************************************************************/
crystalloGetSymmetryInfo(Crystal * crystal,gdouble symprec)2195 gchar* crystalloGetSymmetryInfo(Crystal* crystal, gdouble symprec)
2196 {
2197 return crystalloGetCIF(crystal, symprec, TRUE);
2198 }
2199 /********************************************************************************/
crystalloStandardizeCell(GList ** patoms,gint to_primitive,gint no_idealize,gdouble symprec)2200 gboolean crystalloStandardizeCell(GList** patoms, gint to_primitive, gint no_idealize, gdouble symprec)
2201 {
2202 gboolean ok = FALSE;
2203 GList* newAtoms = NULL;
2204 if(!patoms || !*patoms) return ok;
2205 ok = crystalloCartnToFract(*patoms);
2206 if(ok)
2207 {
2208 ok = FALSE;
2209 newAtoms = crystalloStandardizeCellSPG(*patoms, to_primitive, no_idealize, symprec);
2210 if(newAtoms)
2211 {
2212 ok = TRUE;
2213 g_list_free_all(*patoms, crystalloFreeAtom);
2214 *patoms = newAtoms;
2215 }
2216 }
2217 crystalloFractToCartn(*patoms);
2218 return ok;
2219 }
2220 /*****************************************************************************************************************/
crystalloAllRecObtuse(GList * atoms,gdouble symprec)2221 gboolean crystalloAllRecObtuse(GList* atoms, gdouble symprec)
2222 {
2223
2224 gdouble Tv[3][3];
2225 gint nTv = 0;
2226 gdouble Ts[3][3];
2227 gdouble lattice[3][3];
2228 gdouble alpha, beta, gamma;
2229 gint i,j;
2230
2231 nTv = crystalloGetTv(atoms, Tv);
2232 computeTs(Tv,Ts,nTv);
2233 for(i=0;i<3;i++) for(j=0;j<3;j++) lattice[i][j] = Ts[j][i];
2234 if (0 != spg_niggli_reduce(lattice, (const double) symprec))
2235 {
2236 for(i=0;i<3;i++) for(j=0;j<3;j++) Ts[j][i] = lattice[i][j];
2237 }
2238
2239 alpha = angle(Ts[1],Ts[2]);
2240 beta = angle(Ts[0],Ts[2]);
2241 gamma = angle(Ts[0],Ts[1]);
2242 if(alpha>90 && beta >90 && gamma>90) return TRUE;
2243 return FALSE;
2244 }
2245 /*****************************************************************************************************************/
set3DMatrix(gdouble M[][3],gdouble f,gdouble x00,gdouble x01,gdouble x02,gdouble x10,gdouble x11,gdouble x12,gdouble x20,gdouble x21,gdouble x22)2246 static void set3DMatrix(gdouble M[][3], gdouble f,
2247 gdouble x00, gdouble x01, gdouble x02,
2248 gdouble x10, gdouble x11, gdouble x12,
2249 gdouble x20, gdouble x21, gdouble x22)
2250 {
2251 gint i,j;
2252 M[0][0] = x00; M[0][1] = x01; M[0][2] = x02;
2253 M[1][0] = x10; M[1][1] = x11; M[1][2] = x12;
2254 M[2][0] = x20; M[2][1] = x21; M[2][2] = x22;
2255 for(i=0;i<3;i++)
2256 for(j=0;j<3;j++) M[i][j] *= f;
2257 }
2258 /*****************************************************************************************************************/
getPAndPm1MAtrix(gdouble P[][3],gdouble Pm1[][3],gchar * pearsonSymbol)2259 static void getPAndPm1MAtrix(gdouble P[][3], gdouble Pm1[][3], gchar* pearsonSymbol)
2260 {
2261 gchar s[10];
2262 sprintf(s,"%c%c", pearsonSymbol[0], pearsonSymbol[1]);
2263
2264 if(!strcmp(s, "cP") || !strcmp(s, "tP") || !strcmp(s, "hP") || !strcmp(s, "oP") || !strcmp(s, "mP") || !strcmp(s, "aP") )
2265 {
2266 set3DMatrix(P,1.0,1, 0, 0, 0, 1, 0, 0, 0, 1);
2267 set3DMatrix(Pm1,1.0,1, 0, 0, 0, 1, 0, 0, 0, 1);
2268 }
2269 else if(!strcmp(s, "cF") || !strcmp(s, "oF"))
2270 {
2271 set3DMatrix(P, 1./2.,0, 1, 1, 1, 0, 1, 1, 1, 0);
2272 set3DMatrix(Pm1, 1.0, -1, 1, 1, 1, -1, 1, 1, 1, -1);
2273 }
2274 else if(!strcmp(s, "cI") || !strcmp(s, "tI") || !strcmp(s, "oI"))
2275 {
2276 set3DMatrix(P, 1./2. , -1, 1, 1, 1, -1, 1, 1, 1, -1);
2277 set3DMatrix(Pm1 , 1.0, 0, 1, 1, 1, 0, 1, 1, 1, 0);
2278 }
2279 else if(!strcmp(s, "hR"))
2280 {
2281 set3DMatrix(P , 1./3. ,2, -1, -1, 1, 1, -2, 1, 1, 1);
2282 set3DMatrix( Pm1, 1.0, 1, 0, 1, -1, 1, 1, 0, -1, 1);
2283 }
2284 else if(!strcmp(s, "oC"))
2285 {
2286 set3DMatrix(P, 1./2. , 1, 1, 0, -1, 1, 0, 0, 0, 2);
2287 set3DMatrix(Pm1 , 1.0, 1, -1, 0, 1, 1, 0, 0, 0, 1);
2288 }
2289 else if(!strcmp(s, "oA"))
2290 {
2291 set3DMatrix(P, 1./2. , 0, 0, 2, 1, 1, 0, -1, 1, 0);
2292 set3DMatrix( Pm1 , 1.0, 0, 1, -1, 0, 1, 1, 1, 0, 0);
2293 }
2294 else if(!strcmp(s, "mC"))
2295 {
2296 set3DMatrix(P, 1./2. , 1, -1, 0, 1, 1, 0, 0, 0, 2);
2297 set3DMatrix(Pm1 , 1.0,1, 1, 0, -1, 1, 0, 0, 0, 1);
2298 }
2299 else
2300 {
2301 set3DMatrix(P,1.0,1, 0, 0, 0, 1, 0, 0, 0, 1);
2302 set3DMatrix(Pm1,1.0,1, 0, 0, 0, 1, 0, 0, 0, 1);
2303 }
2304 }
2305 /********************************************************************************/
crystalloPrimitiveCellHinuma(GList ** patoms,gchar * pearsonSymbol)2306 gboolean crystalloPrimitiveCellHinuma(GList** patoms, gchar* pearsonSymbol)
2307 {
2308 gdouble P[3][3];
2309 gdouble Pm1[3][3];
2310 gdouble C[3];
2311 GList* l;
2312 gint i,j;
2313 if(!patoms || !*patoms) return FALSE;
2314 getPAndPm1MAtrix(P, Pm1, pearsonSymbol);
2315 crystalloCartnToFract(*patoms);
2316 for(l = g_list_first(*patoms); l != NULL; l = l->next)
2317 {
2318 CrystalloAtom* a = (CrystalloAtom*)l->data;
2319 for(j=0;j<3;j++) C[j] = a->C[j];
2320 for(j=0;j<3;j++) a->C[j] = 0;
2321
2322 if(strstr(a->symbol,"Tv") || strstr(a->symbol,"TV"))
2323 for(j=0;j<3;j++) for(i=0;i<3;i++) a->C[j] += P[i][j]*C[i];
2324 else
2325 for(j=0;j<3;j++) for(i=0;i<3;i++) a->C[j] += Pm1[j][i]*C[i];
2326 }
2327 crystalloSetAtomsInBox(*patoms);
2328 crystalloFractToCartn(*patoms);
2329 crystalloRemoveAtomsWithSmallDistance(patoms);
2330 return TRUE;
2331 }
2332