1 /* GabeditSPG.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 "../../spglib/spglib.h"
26 #include "../Common/Global.h"
27 #include "../Utils/Constants.h"
28 #include "../Utils/Utils.h"
29 #include "../Utils/AtomsProp.h"
30 #include "../Crystallography/Crystallo.h"
31
32
crystalloGetNiggli(GList * atoms,gdouble newTv[][3],gdouble symprec)33 gboolean crystalloGetNiggli(GList* atoms, gdouble newTv[][3], gdouble symprec)
34 {
35 gint i,j;
36 double lattice[3][3];
37 gint nTv = crystalloGetTv(atoms, newTv);
38 gboolean ok = FALSE;
39 if(nTv<3) return ok;
40 for(i=0;i<3;i++) for(j=0;j<3;j++) lattice[i][j] = newTv[j][i];
41 if (0 != spg_niggli_reduce(lattice, (const double) symprec))
42 {
43 fprintf(stderr," Ok new Tv vector\n");
44 fprintf(stderr," oldTV =\n"); for(i=0;i<3;i++) { for(j=0;j<3;j++) fprintf(stderr,"%f ",newTv[i][j]); fprintf(stderr,"\n"); }
45
46 for(i=0;i<3;i++) for(j=0;j<3;j++) newTv[j][i] = lattice[i][j];
47 fprintf(stderr," newTV =\n"); for(i=0;i<3;i++) { for(j=0;j<3;j++) fprintf(stderr,"%f ",newTv[i][j]); fprintf(stderr,"\n"); }
48 ok = TRUE;
49 }
50 return ok;
51 }
crystalloGetDelaunay(GList * atoms,gdouble newTv[][3],gdouble symprec)52 gboolean crystalloGetDelaunay(GList* atoms, gdouble newTv[][3], gdouble symprec)
53 {
54 gint i,j;
55 double lattice[3][3];
56 gint nTv = crystalloGetTv(atoms, newTv);
57 gboolean ok = FALSE;
58 if(nTv<3) return ok;
59 for(i=0;i<3;i++) for(j=0;j<3;j++) lattice[i][j] = newTv[j][i];
60 if (0 != spg_delaunay_reduce(lattice, (const double) symprec))
61 {
62 for(i=0;i<3;i++) for(j=0;j<3;j++) newTv[j][i] = lattice[i][j];
63 ok = TRUE;
64 }
65 return ok;
66 }
67 /*
68 gboolean crystalloGetPrimitive(GList* atoms, gdouble newTv[][3], gdouble symprec)
69 {
70 double (*position)[3];
71 int* types = NULL;
72 gint i,j;
73 double lattice[3][3];
74 gint nTv = crystalloGetTv(atoms, newTv);
75 GList *l = NULL;
76 int nAtoms = 0;
77 gboolean ok = FALSE;
78
79 if(nTv<3) return ok;
80 for(i=0;i<3;i++) for(j=0;j<3;j++) lattice[i][j] = newTv[j][i];
81
82 nAtoms = 0;
83 for(l = g_list_first(atoms); l != NULL; l = l->next)
84 {
85 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
86 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
87 nAtoms++;
88 }
89 if(nAtoms<1) return ok;
90 position = (double (*)[3]) malloc(sizeof(double[3]) * nAtoms);
91 types = (int*) malloc(sizeof(int) * nAtoms);
92 i=0;
93 for(l = g_list_first(atoms); l != NULL; l = l->next)
94 {
95 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
96 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
97 types[i] = (int)get_atomic_number_from_symbol(crystalloAtom->symbol);
98 for(j=0;j<3;j++) position[i][j] = crystalloAtom->C[j];
99 i++;
100 }
101
102 if( 0!= spg_find_primitive(lattice, position, types, nAtoms, (double)symprec))
103 {
104 for(i=0;i<3;i++) for(j=0;j<3;j++) newTv[j][i] = lattice[i][j];
105 ok=TRUE;
106 }
107
108 if(types) free(types);
109 if(position) free(position);
110 return ok;
111 }
112 */
crystalloGetDataSet(GList * atoms,gdouble symprec)113 SpglibDataset * crystalloGetDataSet(GList* atoms, gdouble symprec)
114 {
115 SpglibDataset * spgDataSet = NULL;
116 double (*position)[3];
117 int* types = NULL;
118 gint i,j;
119 double lattice[3][3];
120 gdouble newTv[3][3];
121 gint nTv = crystalloGetTv(atoms, newTv);
122 GList *l = NULL;
123 int nAtoms = 0;
124 gint numGroup = 0;
125
126 if(nTv<3) return spgDataSet;
127 for(i=0;i<3;i++) for(j=0;j<3;j++) lattice[i][j] = newTv[j][i];
128
129 nAtoms = 0;
130 for(l = g_list_first(atoms); l != NULL; l = l->next)
131 {
132 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
133 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
134 nAtoms++;
135 }
136 if(nAtoms<1) return spgDataSet;
137 position = (double (*)[3]) malloc(sizeof(double[3]) * nAtoms);
138 types = (int*) malloc(sizeof(int) * nAtoms);
139 i=0;
140 for(l = g_list_first(atoms); l != NULL; l = l->next)
141 {
142 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
143 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
144 types[i] = (int)get_atomic_number_from_symbol(crystalloAtom->symbol);
145 for(j=0;j<3;j++) position[i][j] = crystalloAtom->C[j];
146 i++;
147 }
148 spgDataSet = spg_get_dataset(lattice, position, types, nAtoms, (double)symprec);
149
150 if(types) free(types);
151 if(position) free(position);
152 return spgDataSet;
153 }
154
crystalloGetGroupName(GList * atoms,char groupName[],gdouble symprec)155 gint crystalloGetGroupName(GList* atoms, char groupName[], gdouble symprec)
156 {
157 double (*position)[3];
158 int* types = NULL;
159 gint i,j;
160 double lattice[3][3];
161 gdouble newTv[3][3];
162 gint nTv = crystalloGetTv(atoms, newTv);
163 GList *l = NULL;
164 int nAtoms = 0;
165 gint numGroup = 0;
166
167 if(nTv<3) return numGroup;
168 for(i=0;i<3;i++) for(j=0;j<3;j++) lattice[i][j] = newTv[j][i];
169
170 nAtoms = 0;
171 for(l = g_list_first(atoms); l != NULL; l = l->next)
172 {
173 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
174 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
175 nAtoms++;
176 }
177 if(nAtoms<1) return numGroup;
178 position = (double (*)[3]) malloc(sizeof(double[3]) * nAtoms);
179 types = (int*) malloc(sizeof(int) * nAtoms);
180 i=0;
181 for(l = g_list_first(atoms); l != NULL; l = l->next)
182 {
183 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
184 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
185 types[i] = (int)get_atomic_number_from_symbol(crystalloAtom->symbol);
186 for(j=0;j<3;j++) position[i][j] = crystalloAtom->C[j];
187 i++;
188 }
189
190 numGroup = (gint) spg_get_international(groupName, lattice, position, types, nAtoms, (double)symprec);
191 if(types) free(types);
192 if(position) free(position);
193 return numGroup;
194 }
setAtomTypes(CrystalloAtom * crystalAtom,gchar * type)195 static void setAtomTypes(CrystalloAtom* crystalAtom, gchar* type)
196 {
197 sprintf(crystalAtom->mmType,"%s",type);
198 sprintf(crystalAtom->pdbType,"%s",type);
199 sprintf(crystalAtom->residueName,"%s",type);
200 }
crystalloPrimitiveSPG(GList * atoms,gdouble symprec)201 GList* crystalloPrimitiveSPG(GList* atoms, gdouble symprec)
202 {
203 double (*position)[3];
204 int* types = NULL;
205 gint i,j;
206 double lattice[3][3];
207 gdouble newTv[3][3];
208 gint nTv = crystalloGetTv(atoms, newTv);
209 GList *l = NULL;
210 int nAtoms = 0;
211 gint numGroup = 0;
212 GList* newAtoms = NULL;
213 int nAtomsNew = 0;
214
215 if(nTv<3) return newAtoms;
216 for(i=0;i<3;i++) for(j=0;j<3;j++) lattice[i][j] = newTv[j][i];
217
218 nAtoms = 0;
219 for(l = g_list_first(atoms); l != NULL; l = l->next)
220 {
221 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
222 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
223 nAtoms++;
224 }
225 if(nAtoms<1) return newAtoms;
226 position = (double (*)[3]) malloc(sizeof(double[3]) * 4*nAtoms);
227 types = (int*) malloc(sizeof(int) * 4*nAtoms);
228 i=0;
229 for(l = g_list_first(atoms); l != NULL; l = l->next)
230 {
231 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
232 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
233 types[i] = (int)get_atomic_number_from_symbol(crystalloAtom->symbol);
234 for(j=0;j<3;j++) position[i][j] = crystalloAtom->C[j];
235 i++;
236 }
237
238 nAtomsNew = (gint) spg_find_primitive(lattice, position, types, nAtoms, (double)symprec);
239 newAtoms = NULL;
240 i=0;
241 for(i=0;i<nAtomsNew;i++)
242 {
243 CrystalloAtom* crystalAtom = g_malloc(sizeof(CrystalloAtom));
244 gchar* symbol = get_symbol_using_z(types[i]);
245 sprintf(crystalAtom->symbol,"%s", symbol);
246 g_free(symbol);
247 setAtomTypes(crystalAtom, crystalAtom->symbol);
248 crystalAtom->residueNumber=1;
249 crystalAtom->charge=0.0;
250
251 for(j=0;j<3;j++) crystalAtom->C[j]= position[i][j];
252 newAtoms = g_list_append(newAtoms, (gpointer) crystalAtom);
253 }
254 for(i=0;i<3;i++)
255 {
256 CrystalloAtom* crystalAtom = g_malloc(sizeof(CrystalloAtom));
257 sprintf(crystalAtom->symbol,"%s","Tv");
258 setAtomTypes(crystalAtom, crystalAtom->symbol);
259 crystalAtom->residueNumber=1;
260 crystalAtom->charge=0.0;
261
262 for(j=0;j<3;j++) crystalAtom->C[j]= lattice[j][i];
263 newAtoms = g_list_append(newAtoms, (gpointer) crystalAtom);
264 }
265
266 if(types) free(types);
267 if(position) free(position);
268 return newAtoms;
269 }
crystalloStandardizeCellSPG(GList * atoms,gint to_primitive,gint no_idealize,gdouble symprec)270 GList* crystalloStandardizeCellSPG(GList* atoms, gint to_primitive, gint no_idealize, gdouble symprec)
271 {
272 double (*position)[3];
273 int* types = NULL;
274 gint i,j;
275 double lattice[3][3];
276 gdouble newTv[3][3];
277 gint nTv = crystalloGetTv(atoms, newTv);
278 GList *l = NULL;
279 int nAtoms = 0;
280 gint numGroup = 0;
281 GList* newAtoms = NULL;
282 int nAtomsNew = 0;
283
284 if(nTv<3) return newAtoms;
285 for(i=0;i<3;i++) for(j=0;j<3;j++) lattice[i][j] = newTv[j][i];
286
287 nAtoms = 0;
288 for(l = g_list_first(atoms); l != NULL; l = l->next)
289 {
290 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
291 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
292 nAtoms++;
293 }
294 if(nAtoms<1) return newAtoms;
295 position = (double (*)[3]) malloc(sizeof(double[3]) * 4*nAtoms);
296 types = (int*) malloc(sizeof(int) * 4*nAtoms);
297 i=0;
298 for(l = g_list_first(atoms); l != NULL; l = l->next)
299 {
300 CrystalloAtom* crystalloAtom = (CrystalloAtom*)l->data;
301 if(strstr(crystalloAtom->symbol,"Tv") || strstr(crystalloAtom->symbol,"TV")) continue;
302 types[i] = (int)get_atomic_number_from_symbol(crystalloAtom->symbol);
303 for(j=0;j<3;j++) position[i][j] = crystalloAtom->C[j];
304 i++;
305 }
306
307 nAtomsNew = (gint) spg_standardize_cell(lattice, position, types, nAtoms, (const int) to_primitive, (const int)no_idealize, (const double) symprec);
308 newAtoms = NULL;
309 i=0;
310 for(i=0;i<nAtomsNew;i++)
311 {
312 CrystalloAtom* crystalAtom = g_malloc(sizeof(CrystalloAtom));
313 gchar* symbol = get_symbol_using_z(types[i]);
314 sprintf(crystalAtom->symbol,"%s", symbol);
315 g_free(symbol);
316 setAtomTypes(crystalAtom, crystalAtom->symbol);
317 crystalAtom->residueNumber=1;
318 crystalAtom->charge=0.0;
319
320 for(j=0;j<3;j++) crystalAtom->C[j]= position[i][j];
321 newAtoms = g_list_append(newAtoms, (gpointer) crystalAtom);
322 }
323 for(i=0;i<3;i++)
324 {
325 CrystalloAtom* crystalAtom = g_malloc(sizeof(CrystalloAtom));
326 sprintf(crystalAtom->symbol,"%s","Tv");
327 setAtomTypes(crystalAtom, crystalAtom->symbol);
328 crystalAtom->residueNumber=1;
329 crystalAtom->charge=0.0;
330
331 for(j=0;j<3;j++) crystalAtom->C[j]= lattice[j][i];
332 newAtoms = g_list_append(newAtoms, (gpointer) crystalAtom);
333 }
334
335 if(types) free(types);
336 if(position) free(position);
337 return newAtoms;
338 }
339 /********************************************************************************/
standardizeFromDataSetSPG(Crystal * crystal,gdouble symprec)340 SpglibDataset* standardizeFromDataSetSPG(Crystal* crystal, gdouble symprec)
341 {
342 gboolean ok = crystalloCartnToFract(crystal->atoms);
343 SpglibDataset * spgDataSet = crystalloGetDataSet(crystal->atoms, symprec);
344 if(spgDataSet)
345 {
346 GList* l;
347 gint i,j;
348 i=0;
349 GList* newAtoms = NULL;
350 i=0;
351 for(i=0;i<spgDataSet->n_std_atoms;i++)
352 {
353 CrystalloAtom* crystalAtom = g_malloc(sizeof(CrystalloAtom));
354 gchar* symbol = get_symbol_using_z(spgDataSet->std_types[i]);
355 sprintf(crystalAtom->symbol,"%s", symbol);
356 g_free(symbol);
357 setAtomTypes(crystalAtom, crystalAtom->symbol);
358 crystalAtom->residueNumber=1;
359 crystalAtom->charge=0.0;
360
361 for(j=0;j<3;j++) crystalAtom->C[j]= spgDataSet->std_positions[i][j];
362 newAtoms = g_list_append(newAtoms, (gpointer) crystalAtom);
363 }
364 for(i=0;i<3;i++)
365 {
366 CrystalloAtom* crystalAtom = g_malloc(sizeof(CrystalloAtom));
367 sprintf(crystalAtom->symbol,"%s","Tv");
368 setAtomTypes(crystalAtom, crystalAtom->symbol);
369 crystalAtom->residueNumber=1;
370 crystalAtom->charge=0.0;
371
372 for(j=0;j<3;j++) crystalAtom->C[j]= spgDataSet->std_lattice[j][i];
373 newAtoms = g_list_append(newAtoms, (gpointer) crystalAtom);
374 }
375 if(newAtoms)
376 {
377 g_list_foreach (crystal->atoms, (GFunc) crystalloFreeAtom, NULL);
378 g_list_free (crystal->atoms);
379 crystal->atoms = newAtoms;
380 }
381 }
382 if(ok) ok=crystalloFractToCartn(crystal->atoms);
383 return spgDataSet;
384 }
385