1 //
2 // point_group.c
3 // Symmetry
4 //
5 // Created by Marcus Johansson on 14/04/14.
6 // Copyright (c) 2014 Marcus Johansson.
7 //
8 // Distributed under the MIT License ( See LICENSE file or copy at http://opensource.org/licenses/MIT )
9 //
10
11
12 #include <math.h>
13 #include <float.h>
14 #include <stdio.h>
15 #include <string.h>
16 #include <stdlib.h>
17 #include "symmetry.h"
18 #include "linalg.h"
19 #include "point_group.h"
20 #include "permutation.h"
21
22 #include "debug.h"
23
24 #define PHI 1.618033988749894848204586834
25 #ifndef M_PI
26 #define M_PI 3.14159265358979323846264338327950288419716939937510582
27 #endif
28
29 #ifndef M_PI_2
30 #define M_PI_2 (M_PI/2)
31 #endif
32
33
34 #define CLASSIFICATION_THRESHOLD 0.01
35
36 msym_error_t determinePointGroup(int sopsl, msym_symmetry_operation_t *sops, msym_thresholds_t *thresholds, msym_point_group_t *pg);
37
38 msym_error_t generatePointGroup(msym_point_group_type_t type, int n, msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, msym_point_group_t **opg);
39 msym_error_t reorientAxes(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds);
40 msym_error_t transformAxes(msym_point_group_type_t type, int n, msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double transform[3][3]);
41 msym_error_t transformPrimary(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double transform[3][3]);
42 msym_error_t transformSecondary(msym_point_group_type_t type, msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double transform[3][3]);
43 msym_error_t getPointGroupOrder(msym_point_group_type_t type, int n, int *order);
44 msym_error_t getPointGroupName(msym_point_group_type_t type, int n, size_t max, char name[max]);
45
46 msym_error_t findSecondaryAxisSigma(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double r[3]);
47 msym_error_t findSecondaryAxisC2(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double r[3]);
48
49 msym_error_t findSecondaryAxisC4(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double r[3]);
50
51 msym_error_t findSecondaryAxisC2C5(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double r[3]);
52
53
54
55 msym_error_t generateSymmetryOperations(msym_point_group_type_t type, int n, int order, msym_symmetry_operation_t **osops);
56
57
58 msym_error_t generateSymmetryOperationsImpliedRot(int sopsl, msym_symmetry_operation_t sops[sopsl], int order, msym_thresholds_t *thresholds, int *osopsl);
59
60
61 msym_error_t generateSymmetryOperationsCs(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
62 msym_error_t generateSymmetryOperationsCi(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
63 msym_error_t generateSymmetryOperationsDn(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
64 msym_error_t generateSymmetryOperationsDnd(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
65 msym_error_t generateSymmetryOperationsDnh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
66 msym_error_t generateSymmetryOperationsSn(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
67 msym_error_t generateSymmetryOperationsCn(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
68 msym_error_t generateSymmetryOperationsCnv(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
69 msym_error_t generateSymmetryOperationsCnh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
70 msym_error_t generateSymmetryOperationsT(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
71 msym_error_t generateSymmetryOperationsTd(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
72 msym_error_t generateSymmetryOperationsTh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
73 msym_error_t generateSymmetryOperationsO(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
74 msym_error_t generateSymmetryOperationsOh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
75 msym_error_t generateSymmetryOperationsI(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
76 msym_error_t generateSymmetryOperationsIh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
77
78 msym_error_t generateSymmetryOperationsTetrahedral(int l, msym_symmetry_operation_t sops[l], int c2l, msym_symmetry_operation_t c2[c2l], int csl, msym_symmetry_operation_t cs[csl], int c3l, msym_symmetry_operation_t c3[c3l], int *pk);
79 msym_error_t generateSymmetryOperationsOctahedral(int l, msym_symmetry_operation_t sops[l], int c2l, msym_symmetry_operation_t c2[c2l], int c3l, msym_symmetry_operation_t c3[c3l], int c4l, msym_symmetry_operation_t c4[c4l], int *pk);
80 msym_error_t generateSymmetryOperationsIcosahedral(int l, msym_symmetry_operation_t sops[l], int c2l, msym_symmetry_operation_t c2[c2l], int c3l, msym_symmetry_operation_t c3[c3l], int c5l, msym_symmetry_operation_t c5[c5l], int *pk);
81 msym_error_t generateReflectionPlanes(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
82 msym_error_t generateC2Axes(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
83
84 msym_error_t generateSymmetryOperationsUnknown(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla);
85
86
87 msym_error_t pointGroupFromName(const char *name, msym_point_group_t *pg);
88 msym_error_t pointGroupFromType(msym_point_group_type_t type, int n, msym_point_group_t *pg);
89 msym_error_t generatePointGroupFromStruct(msym_point_group_t *pg, double transform[3][3], msym_thresholds_t *thresholds);
90
91 int classifySymmetryOperations(msym_point_group_t *pg);
92 void sortSymmetryOperations(msym_point_group_t *pg, int classes);
93
generatePointGroupFromType(msym_point_group_type_t type,int n,double transform[3][3],msym_thresholds_t * thresholds,msym_point_group_t ** opg)94 msym_error_t generatePointGroupFromType(msym_point_group_type_t type, int n, double transform[3][3], msym_thresholds_t *thresholds, msym_point_group_t **opg){
95 msym_error_t ret = MSYM_SUCCESS;
96 msym_point_group_t *pg = calloc(1,sizeof(msym_point_group_t));
97 if(MSYM_SUCCESS != (ret = pointGroupFromType(type,n,pg))) goto err;
98 if(MSYM_SUCCESS != (ret = generatePointGroupFromStruct(pg, transform, thresholds))) goto err;
99 *opg = pg;
100 return ret;
101 err:
102 free(pg);
103 return ret;
104 }
105
generatePointGroupFromName(const char * name,double transform[3][3],msym_thresholds_t * thresholds,msym_point_group_t ** opg)106 msym_error_t generatePointGroupFromName(const char *name, double transform[3][3], msym_thresholds_t *thresholds, msym_point_group_t **opg){
107 msym_error_t ret = MSYM_SUCCESS;
108 msym_point_group_t *pg = calloc(1,sizeof(msym_point_group_t));
109 if(MSYM_SUCCESS != (ret = pointGroupFromName(name,pg))) goto err;
110 if(MSYM_SUCCESS != (ret = generatePointGroupFromStruct(pg, transform, thresholds))) goto err;
111 *opg = pg;
112 return ret;
113 err:
114 free(pg);
115 return ret;
116 }
117
generatePointGroupFromStruct(msym_point_group_t * pg,double transform[3][3],msym_thresholds_t * thresholds)118 msym_error_t generatePointGroupFromStruct(msym_point_group_t *pg, double transform[3][3], msym_thresholds_t *thresholds){
119 msym_error_t ret = MSYM_SUCCESS;
120
121 if(MSYM_SUCCESS != (ret = generateSymmetryOperations(pg->type, pg->n, pg->order, &pg->sops))) goto err;
122
123 if(isLinearPointGroup(pg)){
124 pg->perm = NULL;
125 } else {
126 if(MSYM_SUCCESS != (ret = findSymmetryOperationPermutations(pg->order,pg->sops, thresholds, &pg->perm))) goto err;
127 }
128
129 memcpy(pg->transform, transform, sizeof(double[3][3]));
130
131 double T[3][3];
132 minv(pg->transform, T);
133
134 for(msym_symmetry_operation_t *s = pg->sops;s < (pg->sops + pg->order);s++){
135 if(pg->primary == NULL || (s->type == PROPER_ROTATION && s->order > pg->primary->order)) pg->primary = s;
136 mvmul(s->v,T,s->v);
137 }
138
139 return ret;
140
141 err:
142 free(pg->sops);
143 pg->sops = NULL;
144 // Need to free findSymmetryOperationPermutations if there is ever a possibility of an error after that call
145 return ret;
146 }
147
pointGroupFromType(msym_point_group_type_t type,int n,msym_point_group_t * pg)148 msym_error_t pointGroupFromType(msym_point_group_type_t type, int n, msym_point_group_t *pg){
149 msym_error_t ret = MSYM_SUCCESS;
150 pg->type = type;
151 switch (pg->type) {
152 case MSYM_POINT_GROUP_TYPE_Cs:
153 case MSYM_POINT_GROUP_TYPE_Ci:
154 n = 1;
155 break;
156 case MSYM_POINT_GROUP_TYPE_T:
157 case MSYM_POINT_GROUP_TYPE_Td:
158 case MSYM_POINT_GROUP_TYPE_Th:
159 n = 3;
160 break;
161 case MSYM_POINT_GROUP_TYPE_O:
162 case MSYM_POINT_GROUP_TYPE_Oh:
163 n = 4;
164 break;
165 case MSYM_POINT_GROUP_TYPE_I:
166 case MSYM_POINT_GROUP_TYPE_Ih:
167 n = 5;
168 break;
169 default:
170 break;
171 }
172 pg->n = n;
173 if(MSYM_SUCCESS != (ret = getPointGroupOrder(pg->type, pg->n, &pg->order))) goto err;
174 if(MSYM_SUCCESS != (ret = getPointGroupName(pg->type, pg->n, sizeof(pg->name)/sizeof(char), pg->name))) goto err;
175
176 err:
177 return ret;
178 }
179
pointGroupFromName(const char * name,msym_point_group_t * pg)180 msym_error_t pointGroupFromName(const char *name, msym_point_group_t *pg){
181 msym_error_t ret = MSYM_SUCCESS;
182 int n = 0, gi = 0, ri = 0;;
183 char g = 0, r = 0;
184
185 int map[7][6];
186
187
188 struct _pg_map {
189 int i;
190 msym_point_group_type_t type;
191 } pg_map[] = {
192 {1, MSYM_POINT_GROUP_TYPE_Cn},
193 {2, MSYM_POINT_GROUP_TYPE_Cnv},
194 {3, MSYM_POINT_GROUP_TYPE_Cnh},
195 {4, MSYM_POINT_GROUP_TYPE_Ci},
196 {5, MSYM_POINT_GROUP_TYPE_Cs},
197 {6, MSYM_POINT_GROUP_TYPE_Dn},
198 {7, MSYM_POINT_GROUP_TYPE_Dnh},
199 {8, MSYM_POINT_GROUP_TYPE_Dnd},
200 {9, MSYM_POINT_GROUP_TYPE_Sn},
201 {10, MSYM_POINT_GROUP_TYPE_T},
202 {11, MSYM_POINT_GROUP_TYPE_Th},
203 {12, MSYM_POINT_GROUP_TYPE_Td},
204 {13, MSYM_POINT_GROUP_TYPE_O},
205 {14, MSYM_POINT_GROUP_TYPE_Oh},
206 {15, MSYM_POINT_GROUP_TYPE_I},
207 {16, MSYM_POINT_GROUP_TYPE_Ih},
208 {17, MSYM_POINT_GROUP_TYPE_K},
209 {18, MSYM_POINT_GROUP_TYPE_Kh}
210
211 };
212
213 memset(map,0,sizeof(int[7][6]));
214
215 map[0][0] = 1;
216 map[0][1] = 2;
217 map[0][2] = 3;
218 map[0][4] = 4;
219 map[0][5] = 5;
220 map[1][0] = 6;
221 map[1][2] = 7;
222 map[1][3] = 8;
223 map[2][0] = 9;
224 map[3][0] = 10;
225 map[3][2] = 11;
226 map[3][3] = 12;
227 map[4][0] = 13;
228 map[4][2] = 14;
229 map[5][0] = 15;
230 map[5][2] = 16;
231 map[6][0] = 17;
232 map[6][2] = 18;
233
234 if(sscanf(name,"%c%d%c",&g,&n,&r) < 2 && sscanf(name,"%c%c",&g,&r) < 1){
235 ret = MSYM_INVALID_POINT_GROUP;
236 msymSetErrorDetails("Invalid point group name %s",name);
237 goto err;
238 }
239
240 if(n < 0) {
241 ret = MSYM_INVALID_POINT_GROUP;
242 msymSetErrorDetails("Invalid point group order %d",n);
243 goto err;
244 }
245
246 switch(g){
247 case 'C' : gi = 0; break;
248 case 'D' : gi = 1; break;
249 case 'S' : {
250 if(n < 4 || n % 2 != 0){
251 ret = MSYM_INVALID_POINT_GROUP;
252 msymSetErrorDetails("Improper rotation order (%d) must be even and >= 4",n);
253 goto err;
254 }
255 gi = 2;
256 break;
257 }
258 case 'T' : gi = 3; break;
259 case 'O' : gi = 4; break;
260 case 'I' : gi = 5; break;
261 case 'K' : gi = 6; break;
262 default :
263 ret = MSYM_INVALID_POINT_GROUP;
264 msymSetErrorDetails("Invalid point group type %c",g);
265 goto err;
266 }
267
268 switch(r){
269 case 0 : ri = 0; break;
270 case 'v' : ri = 1; break;
271 case 'h' : ri = 2; break;
272 case 'd' : ri = 3; break;
273 case 'i' : ri = 4; break;
274 case 's' : ri = 5; break;
275 default :
276 ret = MSYM_INVALID_POINT_GROUP;
277 msymSetErrorDetails("Invalid point group subtype %c",r);
278 goto err;
279 }
280
281 int fi, fil = sizeof(pg_map)/sizeof(pg_map[0]);
282 for(fi = 0;fi < fil;fi++){
283 if(pg_map[fi].i == map[gi][ri]) break;
284 }
285
286 if(fi == fil){
287 ret = MSYM_INVALID_POINT_GROUP;
288 msymSetErrorDetails("Cannot find point group %s",name);
289 goto err;
290 }
291
292 return pointGroupFromType(pg_map[fi].type, n, pg);
293 err:
294 return ret;
295
296 }
297
298
getPointGroupName(msym_point_group_type_t type,int n,size_t max,char name[max])299 msym_error_t getPointGroupName(msym_point_group_type_t type, int n, size_t max, char name[max]){
300 msym_error_t ret = MSYM_SUCCESS;
301 switch(type) {
302 case MSYM_POINT_GROUP_TYPE_Ci : snprintf(name,max,"Ci"); break;
303 case MSYM_POINT_GROUP_TYPE_Cs : snprintf(name,max,"Cs"); break;
304 case MSYM_POINT_GROUP_TYPE_Cn : snprintf(name,max,"C%d",n); break;
305 case MSYM_POINT_GROUP_TYPE_Cnh : snprintf(name,max,"C%dh",n); break;
306 case MSYM_POINT_GROUP_TYPE_Cnv : snprintf(name,max,"C%dv",n); break;
307 case MSYM_POINT_GROUP_TYPE_Dn : snprintf(name,max,"D%d",n); break;
308 case MSYM_POINT_GROUP_TYPE_Dnh : snprintf(name,max,"D%dh",n); break;
309 case MSYM_POINT_GROUP_TYPE_Dnd : snprintf(name,max,"D%dd",n); break;
310 case MSYM_POINT_GROUP_TYPE_Sn : snprintf(name,max,"S%d",n); break;
311 case MSYM_POINT_GROUP_TYPE_T : snprintf(name,max,"T"); break;
312 case MSYM_POINT_GROUP_TYPE_Td : snprintf(name,max,"Td"); break;
313 case MSYM_POINT_GROUP_TYPE_Th : snprintf(name,max,"Th"); break;
314 case MSYM_POINT_GROUP_TYPE_O : snprintf(name,max,"O"); break;
315 case MSYM_POINT_GROUP_TYPE_Oh : snprintf(name,max,"Oh"); break;
316 case MSYM_POINT_GROUP_TYPE_I : snprintf(name,max,"I"); break;
317 case MSYM_POINT_GROUP_TYPE_Ih : snprintf(name,max,"Ih"); break;
318 case MSYM_POINT_GROUP_TYPE_K : snprintf(name,max,"K"); break;
319 case MSYM_POINT_GROUP_TYPE_Kh : snprintf(name,max,"Kh"); break;
320 default :
321 msymSetErrorDetails("Unknown point group when determining name");
322 ret = MSYM_POINT_GROUP_ERROR;
323 goto err;
324 }
325 err:
326 return ret;
327 }
328
329
generatePointGroup(msym_point_group_type_t type,int n,msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,msym_point_group_t ** opg)330 msym_error_t generatePointGroup(msym_point_group_type_t type, int n, msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, msym_point_group_t **opg){
331 msym_error_t ret = MSYM_SUCCESS;
332 msym_point_group_t *pg = calloc(1,sizeof(msym_point_group_t));
333 pg->type = type;
334 pg->n = n;
335 if(MSYM_SUCCESS != (ret = getPointGroupName(type,n,sizeof(pg->name)/sizeof(char),pg->name))) goto err;
336 if(MSYM_SUCCESS != (ret = getPointGroupOrder(type, n, &pg->order))) goto err;
337 if(pg->order < sopsl){
338 ret = MSYM_POINT_GROUP_ERROR;
339 msymSetErrorDetails("More symmetry operations than order of point group (%s). Order: %d Number of operations: %d",pg->name, pg->order,sopsl);
340 goto err;
341 }
342
343
344 //TODO: this may still be needed, in broken symmetry
345 //if(MSYM_SUCCESS != (ret = generateSymmetryOperationsImpliedRot(sopsl, sops, pg->order, thresholds, &sopsl))) goto err;
346
347 if(MSYM_SUCCESS != (ret = transformAxes(type, n, primary, sopsl, sops, thresholds, pg->transform))) goto err;
348
349 if(MSYM_SUCCESS != (ret = generateSymmetryOperations(type, n, pg->order, &pg->sops))) goto err;
350
351 if(isLinearPointGroup(pg)){
352 pg->perm = NULL;
353 } else {
354 if(MSYM_SUCCESS != (ret = findSymmetryOperationPermutations(pg->order,pg->sops, thresholds, &pg->perm))) goto err;
355 }
356
357 double T[3][3];
358 minv(pg->transform, T);
359
360 for(int i = 0; i < pg->order;i++){
361 mvmul(pg->sops[i].v,T,pg->sops[i].v);
362 if(pg->sops[i].type == PROPER_ROTATION){
363 if(pg->primary == NULL || pg->sops[i].order > pg->primary->order) pg->primary = &(pg->sops[i]);
364 }
365 }
366
367 for(int i = 0;i < pg->order;i++){
368 printSymmetryOperation(&pg->sops[i]);
369 }
370
371 *opg = pg;
372 return ret;
373
374 err:
375 if(pg) free(pg->sops);
376 free(pg);
377 return ret;
378 }
379
getPointGroupOrder(msym_point_group_type_t type,int n,int * order)380 msym_error_t getPointGroupOrder(msym_point_group_type_t type, int n, int *order){
381
382 msym_error_t ret = MSYM_SUCCESS;
383 switch (type) {
384 case (MSYM_POINT_GROUP_TYPE_Cs) :
385 case (MSYM_POINT_GROUP_TYPE_Ci) : *order = 2; break;
386 case (MSYM_POINT_GROUP_TYPE_Cn) :
387 case (MSYM_POINT_GROUP_TYPE_Sn) : *order = n; break;
388 case (MSYM_POINT_GROUP_TYPE_Cnh) :
389 case (MSYM_POINT_GROUP_TYPE_Dn) : *order = 2*n; break;
390 case (MSYM_POINT_GROUP_TYPE_Cnv) : *order = (n == 0 ? 2 : 2*n); break;
391 case (MSYM_POINT_GROUP_TYPE_Dnh) : *order = (n == 0 ? 4 : 4*n); break;
392 case (MSYM_POINT_GROUP_TYPE_Dnd) : *order = 4*n; break;
393 case (MSYM_POINT_GROUP_TYPE_T) : *order = 12; break;
394 case (MSYM_POINT_GROUP_TYPE_Td) :
395 case (MSYM_POINT_GROUP_TYPE_Th) :
396 case (MSYM_POINT_GROUP_TYPE_O) : *order = 24; break;
397 case (MSYM_POINT_GROUP_TYPE_Oh) : *order = 48; break;
398 case (MSYM_POINT_GROUP_TYPE_I) : *order = 60; break;
399 case (MSYM_POINT_GROUP_TYPE_Ih) : *order = 120; break;
400 case (MSYM_POINT_GROUP_TYPE_K) :
401 case (MSYM_POINT_GROUP_TYPE_Kh) : *order = 0;
402 default :
403 msymSetErrorDetails("Point group has no primary axis for reorientation");
404 goto err;
405 }
406 err:
407 return ret;
408 }
409
findPointGroup(int sopsl,msym_symmetry_operation_t * sops,msym_thresholds_t * thresholds,msym_point_group_t ** pg)410 msym_error_t findPointGroup(int sopsl, msym_symmetry_operation_t *sops, msym_thresholds_t *thresholds, msym_point_group_t **pg){
411
412 msym_error_t ret = MSYM_SUCCESS;
413
414 msym_point_group_t tpg = {.type = MSYM_POINT_GROUP_TYPE_Kh, .n = 0, .primary = NULL};
415
416 msym_symmetry_operation_t *tsops = NULL;
417
418 int tsopsl = sopsl;
419
420 if(MSYM_SUCCESS != (ret = determinePointGroup(sopsl, sops, thresholds, &tpg))) goto err;
421 if(MSYM_SUCCESS != (ret = getPointGroupOrder(tpg.type, tpg.n, &tpg.order))) goto err;
422
423 if(tpg.order < sopsl){
424 int length = 2*sopsl > 121 ? 2*sopsl : 121;
425 tsops = calloc(length, sizeof(msym_symmetry_operation_t));
426 memcpy(tsops, sops, sizeof(msym_symmetry_operation_t[sopsl]));
427 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsImpliedRot(tsopsl, tsops, length, thresholds, &tsopsl))) goto err;
428 tpg.primary = NULL;
429 tpg.n = 0;
430 tpg.type = MSYM_POINT_GROUP_TYPE_Kh;
431 if(MSYM_SUCCESS != (ret = determinePointGroup(tsopsl, tsops, thresholds, &tpg))) goto err;
432 if(MSYM_SUCCESS != (ret = getPointGroupOrder(tpg.type, tpg.n, &tpg.order))) goto err;
433
434 if(tpg.order < tsopsl){
435 char buf[4] = {0,0,0,0};
436 getPointGroupName(tpg.type, tpg.n, sizeof(buf),buf);
437 ret = MSYM_POINT_GROUP_ERROR;
438 msymSetErrorDetails("More symmetry operations than order of point group (%s). Order: %d Number of operations: %d", buf, tpg.order,tsopsl);
439 goto err;
440 }
441
442 if(MSYM_SUCCESS != (ret = generatePointGroup(tpg.type, tpg.n, tpg.primary, tsopsl, tsops, thresholds, pg))) goto err;
443
444 } else {
445 if(MSYM_SUCCESS != (ret = generatePointGroup(tpg.type, tpg.n, tpg.primary, sopsl, sops, thresholds, pg))) goto err;
446 }
447
448 err:
449 free(tsops);
450 return ret;
451
452 }
453
determinePointGroup(int sopsl,msym_symmetry_operation_t * sops,msym_thresholds_t * thresholds,msym_point_group_t * pg)454 msym_error_t determinePointGroup(int sopsl, msym_symmetry_operation_t *sops, msym_thresholds_t *thresholds, msym_point_group_t *pg){
455 msym_error_t ret = MSYM_SUCCESS;
456 int inversion = 0, sigma = 0, nC[6] = {0,0,0,0,0,0}, linear = 0;
457 msym_symmetry_operation_t *primary = NULL;
458 msym_symmetry_operation_t *s = NULL;
459
460 for(int i = 0;i < sopsl;i++){
461 if(sops[i].type == PROPER_ROTATION && sops[i].order == 0){
462 linear = 1;
463 break;
464 } else if (sops[i].type == PROPER_ROTATION && sops[i].order > 2){
465 break;
466 }
467 }
468 if(!linear){
469 for(int i = 0;i < sopsl;i++){
470 switch(sops[i].type){
471 case PROPER_ROTATION :
472 if(primary == NULL || sops[i].order > primary->order) primary = &(sops[i]);
473 if(sops[i].order <= 5) nC[sops[i].order]++;
474 break;
475 case INVERSION :
476 inversion = 1;
477 break;
478 case REFLECTION :
479 sigma = 1;
480 break;
481 case IMPROPER_ROTATION :
482 if(s == NULL || sops[i].order > s->order) s = &(sops[i]);
483 break;
484 default :
485 break;
486
487 }
488 }
489
490 pg->n = NULL == primary ? 1 : primary->order;
491 pg->primary = primary;
492
493 if(nC[3] >= 2) {
494 if(nC[5] >= 2){
495 if(inversion){
496 pg->type = MSYM_POINT_GROUP_TYPE_Ih;
497 } else {
498 pg->type = MSYM_POINT_GROUP_TYPE_I;
499 }
500 } else if (nC[4] >= 2) {
501 if(inversion){
502 pg->type = MSYM_POINT_GROUP_TYPE_Oh;
503 } else {
504 pg->type = MSYM_POINT_GROUP_TYPE_O;
505 }
506
507 } else if (sigma){
508 if(inversion){
509 pg->type = MSYM_POINT_GROUP_TYPE_Th;
510 } else {
511 pg->type = MSYM_POINT_GROUP_TYPE_Td;
512 }
513 } else {
514 pg->type = MSYM_POINT_GROUP_TYPE_T;
515 }
516
517 } else if (primary == NULL){
518 if(sigma){
519 pg->type = MSYM_POINT_GROUP_TYPE_Cs;
520 pg->n = 1;
521 } else if(inversion){
522 pg->type = MSYM_POINT_GROUP_TYPE_Ci;
523 pg->n = 1;
524 } else {
525 pg->type = MSYM_POINT_GROUP_TYPE_Cn;
526 pg->n = 1;
527 pg->sops = NULL;
528 pg->order = 0;
529 }
530 } else {
531 int nC2 = 0;
532 int sigma_h = 0;
533 int nsigma_v = 0;
534
535 //Special case for D2d
536 if(primary->order == 2 && s != NULL && !vparallel(primary->v, s->v, thresholds->angle)){
537 for(int i = 0; i < sopsl;i++){
538 if(sops[i].type == PROPER_ROTATION && sops[i].order == 2 && vparallel(sops[i].v, s->v, thresholds->angle)){
539 primary = &sops[i];
540 break;
541 }
542 }
543 }
544
545 for(int i = 0; i < sopsl;i++){
546 nC2 += sops[i].type == PROPER_ROTATION && sops[i].order == 2 && vperpendicular(primary->v,sops[i].v, thresholds->angle);
547 sigma_h = sigma_h || (sops[i].type == REFLECTION && vparallel(primary->v,sops[i].v,thresholds->angle));
548 nsigma_v += (sops[i].type == REFLECTION && vperpendicular(primary->v,sops[i].v,thresholds->angle));
549 }
550
551 pg->n = primary->order;
552 pg->primary = primary;
553
554 if(nC2){ //actually nC2 == primary->order but less is acceptable here since we can generate the rest
555 if(sigma_h){
556 pg->type = MSYM_POINT_GROUP_TYPE_Dnh;
557 } else if (nsigma_v){ //actually nsigma_v == primary->order but less is acceptable here since we can generate the rest
558 pg->type = MSYM_POINT_GROUP_TYPE_Dnd;
559 } else {
560 pg->type = MSYM_POINT_GROUP_TYPE_Dn;
561 }
562
563 } else if (sigma_h) {
564 pg->type = MSYM_POINT_GROUP_TYPE_Cnh;
565 } else if(nsigma_v) { //actually nsigma_v == primary->order but less is acceptable here since we can generate the rest
566 pg->type = MSYM_POINT_GROUP_TYPE_Cnv;
567 } else if(s != NULL){
568 pg->n = s->order;
569 pg->type = MSYM_POINT_GROUP_TYPE_Sn;
570 } else {
571 pg->type = MSYM_POINT_GROUP_TYPE_Cn;
572 }
573 }
574 } else {
575 for(int i = 0; i < sopsl;i++){
576 inversion = inversion || sops[i].type == INVERSION;
577
578 if(sops[i].type == PROPER_ROTATION && (sops[i].order == 0 || primary == NULL)){ primary = &(sops[i]);
579 }
580 }
581
582 pg->n = primary->order;
583 pg->primary = primary;
584
585 if(inversion){
586 pg->type = MSYM_POINT_GROUP_TYPE_Dnh;
587
588 } else {
589 pg->type = MSYM_POINT_GROUP_TYPE_Cnv;
590 }
591 }
592
593 return ret;
594
595 //err:
596 // return ret;
597
598 }
599
600
findSubgroup(msym_subgroup_t * subgroup,msym_thresholds_t * thresholds)601 msym_error_t findSubgroup(msym_subgroup_t *subgroup, msym_thresholds_t *thresholds){
602 msym_error_t ret = MSYM_SUCCESS;
603 int inversion = 0, sigma = 0, nC[6] = {0,0,0,0,0,0}, linear = 0;
604 msym_symmetry_operation_t *primary = NULL;
605 msym_symmetry_operation_t *s = NULL, *sop = NULL;;
606
607
608 for(int i = 0;i < subgroup->order;i++){
609 if(subgroup->sops[i]->type == PROPER_ROTATION && subgroup->sops[i]->order == 0){
610 linear = 1;
611 break;
612 } else if (subgroup->sops[i]->type == PROPER_ROTATION && subgroup->sops[i]->order > 2){
613 break;
614 }
615 }
616 if(!linear){
617 for(int i = 0;i < subgroup->order;i++){
618 sop = subgroup->sops[i];
619 if(sop->power > 1) continue;
620 switch(subgroup->sops[i]->type){
621 case PROPER_ROTATION :
622 if(primary == NULL || sop->order > primary->order) primary = sop;
623 if(sop->order <= 5) nC[sop->order]++;
624 break;
625 case INVERSION :
626 inversion = 1;
627 break;
628 case REFLECTION :
629 sigma = 1;
630 break;
631 case IMPROPER_ROTATION :
632 if(s == NULL || sop->order > s->order) s = sop;
633 break;
634 default :
635 break;
636
637 }
638 }
639 if(nC[3] >= 2) {
640 if(nC[5] >= 2){
641 if(inversion){
642 subgroup->type = MSYM_POINT_GROUP_TYPE_Ih;
643 subgroup->n = primary->order;
644 } else {
645 subgroup->type = MSYM_POINT_GROUP_TYPE_I;
646 subgroup->n = primary->order;
647 }
648 } else if (nC[4] >= 2) {
649 if(inversion){
650 subgroup->type = MSYM_POINT_GROUP_TYPE_Oh;
651 subgroup->n = primary->order;
652 } else {
653 subgroup->type = MSYM_POINT_GROUP_TYPE_O;
654 subgroup->n = primary->order;
655 }
656
657 } else if (sigma){
658 if(inversion){
659 subgroup->type = MSYM_POINT_GROUP_TYPE_Th;
660 subgroup->n = primary->order;
661 } else {
662 subgroup->type = MSYM_POINT_GROUP_TYPE_Td;
663 subgroup->n = primary->order;
664 }
665 } else {
666 subgroup->type = MSYM_POINT_GROUP_TYPE_T;
667 subgroup->n = primary->order;
668 }
669
670 } else if (primary == NULL){
671 if(sigma){
672 subgroup->type = MSYM_POINT_GROUP_TYPE_Cs;
673 subgroup->n = 1;
674 } else if(inversion){
675 subgroup->type = MSYM_POINT_GROUP_TYPE_Ci;
676 subgroup->n = 1;
677 } else {
678 subgroup->type = MSYM_POINT_GROUP_TYPE_Cn;
679 subgroup->n = 1;
680 }
681 } else {
682 int nC2 = 0;
683 int sigma_h = 0;
684 int nsigma_v = 0;
685
686 if(primary->order == 2 && s != NULL && !vparallel(primary->v, s->v, thresholds->angle)){
687 for(int i = 0; i < subgroup->order;i++){
688 sop = subgroup->sops[i];
689 if(sop->power > 1) continue;
690 if(sop->type == PROPER_ROTATION && sop->order == 2 && vparallel(sop->v, s->v, thresholds->angle)){
691 primary = sop;
692 break;
693 }
694 }
695 }
696
697 for(int i = 0; i < subgroup->order;i++){
698 sop = subgroup->sops[i];
699 if(sop->power > 1) continue;
700 nC2 += sop->type == PROPER_ROTATION && sop->order == 2 && vperpendicular(primary->v,sop->v, thresholds->angle);
701 sigma_h = sigma_h || (sop->type == REFLECTION && vparallel(primary->v,sop->v,thresholds->angle));
702 nsigma_v += (sop->type == REFLECTION && vperpendicular(primary->v,sop->v,thresholds->angle));
703 }
704 if(nC2 == primary->order){
705 if(sigma_h){
706 subgroup->type = MSYM_POINT_GROUP_TYPE_Dnh;
707 subgroup->n = primary->order;
708 } else if (nsigma_v == primary->order){
709 subgroup->type = MSYM_POINT_GROUP_TYPE_Dnd;
710 subgroup->n = primary->order;
711 } else {
712 subgroup->type = MSYM_POINT_GROUP_TYPE_Dn;
713 subgroup->n = primary->order;
714 }
715
716 } else if (sigma_h) {
717 subgroup->type = MSYM_POINT_GROUP_TYPE_Cnh;
718 subgroup->n = primary->order;
719 } else if(nsigma_v == primary->order) {
720 subgroup->type = MSYM_POINT_GROUP_TYPE_Cnv;
721 subgroup->n = primary->order;
722 } else if(s != NULL){
723 subgroup->type = MSYM_POINT_GROUP_TYPE_Sn;
724 subgroup->n = s->order;
725 } else {
726 subgroup->type = MSYM_POINT_GROUP_TYPE_Cn;
727 subgroup->n = primary->order;
728 }
729 }
730 } else {
731 for(int i = 0; i < subgroup->order;i++){
732 inversion = inversion || subgroup->sops[i]->type == INVERSION;
733
734 if(subgroup->sops[i]->type == PROPER_ROTATION && (subgroup->sops[i]->order == 0 || primary == NULL)){
735 primary = subgroup->sops[i];
736 }
737 }
738
739 if(inversion){
740 subgroup->type = MSYM_POINT_GROUP_TYPE_Dnh;
741 subgroup->n = primary->order;
742 } else {
743 subgroup->type = MSYM_POINT_GROUP_TYPE_Cnv;
744 subgroup->n = primary->order;
745 }
746 }
747
748 subgroup->primary = primary;
749 if(MSYM_SUCCESS != (ret = getPointGroupName(subgroup->type, subgroup->n, sizeof(subgroup->name)/sizeof(char), subgroup->name))) goto err;
750 return ret;
751
752
753 err:
754 return ret;
755
756 }
757
transformAxes(msym_point_group_type_t type,int n,msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,double transform[3][3])758 msym_error_t transformAxes(msym_point_group_type_t type, int n, msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double transform[3][3]){
759 msym_error_t ret = MSYM_SUCCESS;
760 switch (type){
761 case (MSYM_POINT_GROUP_TYPE_Ci) :
762 case (MSYM_POINT_GROUP_TYPE_K) :
763 case (MSYM_POINT_GROUP_TYPE_Kh) :
764 break;
765 case (MSYM_POINT_GROUP_TYPE_Cs) :
766 for(primary = sops; primary < (primary + sopsl) && primary->type != REFLECTION; primary++){};
767 case (MSYM_POINT_GROUP_TYPE_Cn) :
768 case (MSYM_POINT_GROUP_TYPE_Cnh) :
769 case (MSYM_POINT_GROUP_TYPE_Sn) :
770 if(MSYM_SUCCESS != (ret = reorientAxes(primary, sopsl, sops, thresholds))) goto err;
771 if(MSYM_SUCCESS != (ret = transformPrimary(primary, sopsl, sops, thresholds, transform))) goto err;
772 break;
773 case (MSYM_POINT_GROUP_TYPE_Cnv) :
774 case (MSYM_POINT_GROUP_TYPE_Dnh) :
775 if(MSYM_SUCCESS != (ret = reorientAxes(primary, sopsl, sops, thresholds))) goto err;
776 if(MSYM_SUCCESS != (ret = transformPrimary(primary, sopsl, sops, thresholds, transform))) goto err;
777 if(n > 0){
778 if(MSYM_SUCCESS != (ret = transformSecondary(type, primary, sopsl, sops, thresholds, transform))) goto err;
779 }
780 break;
781 case (MSYM_POINT_GROUP_TYPE_Dn) :
782 case (MSYM_POINT_GROUP_TYPE_Dnd) :
783 case (MSYM_POINT_GROUP_TYPE_O) :
784 case (MSYM_POINT_GROUP_TYPE_Oh) :
785 if(MSYM_SUCCESS != (ret = reorientAxes(primary, sopsl, sops, thresholds))) goto err;
786 if(MSYM_SUCCESS != (ret = transformPrimary(primary, sopsl, sops, thresholds, transform))) goto err;
787 if(MSYM_SUCCESS != (ret = transformSecondary(type, primary, sopsl, sops, thresholds, transform))) goto err;
788 break;
789 case (MSYM_POINT_GROUP_TYPE_T) :
790 case (MSYM_POINT_GROUP_TYPE_Td) :
791 case (MSYM_POINT_GROUP_TYPE_Th) :
792 case (MSYM_POINT_GROUP_TYPE_I) :
793 case (MSYM_POINT_GROUP_TYPE_Ih) : {
794 msym_symmetry_operation_t *dprimary = NULL;
795 msym_symmetry_operation_t *sop;
796 double z = -2.0;
797 for(sop = sops; sop < (sops + sopsl); sop++){
798 if(sop->type == PROPER_ROTATION && sop->order == 2){
799 double v[3];
800 vnorm2(sop->v,v);
801 if(v[2] > z || dprimary == NULL) {
802 dprimary = sop;
803 z = sop->v[2];
804 }
805 }
806 }
807 if(dprimary == NULL) {
808 msymSetErrorDetails("Cannot find C2 axis for point group symmetrization of polyhedral point group");
809 ret = MSYM_POINT_GROUP_ERROR;
810 goto err;
811 }
812 if(MSYM_SUCCESS != (ret = reorientAxes(dprimary, sopsl, sops, thresholds))) goto err;
813 if(MSYM_SUCCESS != (ret = transformPrimary(dprimary, sopsl, sops, thresholds, transform))) goto err;
814 if(MSYM_SUCCESS != (ret = transformSecondary(type, dprimary, sopsl, sops, thresholds, transform))) goto err;
815 break;
816 }
817
818 }
819
820 err:
821 return ret;
822 }
823
isLinearPointGroup(msym_point_group_t * pg)824 int isLinearPointGroup(msym_point_group_t *pg){
825 return 0 == pg->n && (MSYM_POINT_GROUP_TYPE_Dnh == pg->type || MSYM_POINT_GROUP_TYPE_Cnv == pg->type);
826 }
827
isLinearSubgroup(msym_point_group_t * pg)828 int isLinearSubgroup(msym_point_group_t *pg){
829 int sub = 0;
830 switch(pg->type){
831 case MSYM_POINT_GROUP_TYPE_Cnv: sub = pg->n == 0 && pg->order > 2; break;
832 case MSYM_POINT_GROUP_TYPE_Dnh: sub = pg->n == 0 && pg->order > 4; break;
833 default: break;
834 }
835 return sub;
836 }
837
reduceLinearPointGroup(msym_point_group_t * pg,int n,msym_thresholds_t * thresholds)838 msym_error_t reduceLinearPointGroup(msym_point_group_t *pg, int n, msym_thresholds_t *thresholds){
839 msym_error_t ret = MSYM_SUCCESS;
840 int order = 0;
841 msym_permutation_t *perm = NULL;
842 msym_symmetry_operation_t *sops = NULL;
843 msym_symmetry_operation_t *primary = NULL;
844 if(!isLinearPointGroup(pg)){
845 msymSetErrorDetails("Trying to reduce non linear point group");
846 ret = MSYM_POINT_GROUP_ERROR;
847 goto err;
848 }
849
850 if(n == 0) n = 2;
851
852 if(MSYM_SUCCESS != (ret = getPointGroupOrder(pg->type, n, &order))) goto err;
853
854 if(MSYM_SUCCESS != (ret = generateSymmetryOperations(pg->type, 0, order, &sops))) goto err;
855
856 for(int i = 0;i < order;i++){
857 if(PROPER_ROTATION == sops[i].type && n == sops[i].order && HORIZONTAL == sops[i].orientation && 1 == sops[i].power){
858 primary = &sops[i];
859 break;
860 }
861 }
862
863 if(NULL == primary){
864 msymSetErrorDetails("Cannot find primary axis when reducing linear group");
865 ret = MSYM_POINT_GROUP_ERROR;
866 goto err;
867 }
868
869
870 double T[3][3];
871 minv(pg->transform, T);
872
873 for(int i = 0; i < order;i++){
874 mvmul(sops[i].v,T,sops[i].v);
875 }
876
877
878
879 if(MSYM_SUCCESS != (ret = findSymmetryOperationPermutations(order,sops, thresholds, &perm))) goto err;
880
881 for(int i = 0;i < pg->order && pg->perm != NULL;i++){
882 freePermutationData(&pg->perm[i]);
883 }
884
885 free(pg->sops);
886
887
888 pg->sops = sops;
889 pg->order = order;
890 pg->perm = perm;
891 pg->primary = primary;
892
893 return ret;
894 err:
895 free(sops);
896 return ret;
897
898 }
899
900 /* Really should split the orientation and class from the symmetry operation
901 * and move it to the group so we don't have to regenerate */
pointGroupFromSubgroup(const msym_subgroup_t * sg,msym_thresholds_t * thresholds,msym_point_group_t ** opg)902 msym_error_t pointGroupFromSubgroup(const msym_subgroup_t *sg, msym_thresholds_t *thresholds, msym_point_group_t **opg){
903 msym_error_t ret = MSYM_SUCCESS;
904 *opg = calloc(1,sizeof(msym_point_group_t));
905 msym_point_group_t *pg = *opg;
906 pg->type = sg->type;
907 pg->n = sg->n;
908 pg->sops = malloc(sizeof(msym_symmetry_operation_t[sg->order]));
909 memcpy(pg->name,sg->name,sizeof(pg->name));
910
911 if(MSYM_SUCCESS != (ret = getPointGroupOrder(pg->type, pg->n, &pg->order))) goto err;
912 if(pg->order != sg->order){
913 msymSetErrorDetails("Point group order %d does not equal nuber of operations in subgroup %d for point gropu %s",pg->order,sg->order,pg->name);
914 ret = MSYM_POINT_GROUP_ERROR;
915 goto err;
916 }
917
918 for(int i = 0;i < sg->order;i++){
919 if(sg->primary == sg->sops[i]) pg->primary = &pg->sops[i];
920 memcpy(&pg->sops[i], sg->sops[i], sizeof(msym_symmetry_operation_t));
921 }
922
923 mleye(3, pg->transform);
924
925 if(MSYM_SUCCESS != (ret = transformAxes(pg->type, pg->n, pg->primary, pg->order, pg->sops, thresholds, pg->transform))) goto err;
926
927 /* Unfortunately we need to regenerate these as they need a specific
928 * class ordering for orbital symmetrization */
929 free(pg->sops);
930 pg->sops = NULL;
931 pg->primary = NULL;
932
933 if(MSYM_SUCCESS != (ret = generateSymmetryOperations(pg->type, pg->n, pg->order, &pg->sops))) goto err;
934
935 if(isLinearPointGroup(pg)){
936 pg->perm = NULL;
937 } else {
938 if(MSYM_SUCCESS != (ret = findSymmetryOperationPermutations(pg->order,pg->sops, thresholds, &pg->perm))) goto err;
939 }
940
941 double T[3][3];
942 minv(pg->transform, T);
943
944 for(int i = 0; i < pg->order;i++){
945 // vector and orientation may not be equal
946 if(NULL != sg->primary &&
947 NULL == pg->primary &&
948 pg->sops[i].type == sg->primary->type &&
949 pg->sops[i].order == sg->primary->order &&
950 pg->sops[i].power == sg->primary->power){
951
952 pg->primary = &pg->sops[i];
953 }
954 mvmul(pg->sops[i].v,T,pg->sops[i].v);
955 }
956
957 return ret;
958 err:
959 *opg = NULL;
960 free(pg->sops);
961 free(pg);
962 return ret;
963 }
964
965
966 /* Point primary axes above xy plane and all symops above the primary plane.
967 * Not really needed but we may not have to deal with transforms that flip things over
968 * This could be improved with some thresholds, since some things lie in the plane
969 * and get flipped. Should try to align with x as well
970 */
reorientAxes(msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds)971 msym_error_t reorientAxes(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds
972 ){
973 double x[3] = {1.0,0.0,0.0}, y[3] = {0.0,1.0,0.0}, z[3] = {0.0, 0.0, 1.0};
974
975 if(primary == NULL) goto err;
976
977 if(vdot(primary->v,z) < 0.0) vinv(primary->v);
978
979 for(msym_symmetry_operation_t *sop = sops; sop < (sops + sopsl); sop++){
980 if(vperpendicular(sop->v, z, thresholds->angle)){
981 double proj = vdot(sop->v, y)/vabs(sop->v);
982 if(fabs(fabs(proj)-1.0) < thresholds->angle && proj < 0.0) { //along y axis
983 vinv(sop->v);
984 } else if (vdot(sop->v,x) < 0.0){ //in xy-plane not in y axis, reorient to positive x
985 vinv(sop->v);
986 }
987
988 } else if(vdot(primary->v,sop->v) < 0.0) vinv(sop->v); //not in xy-plane reorient to positive along primary axiz
989 }
990
991 return MSYM_SUCCESS;
992
993 err:
994 msymSetErrorDetails("Point group has no primary axis for reorientation");
995 return MSYM_POINT_GROUP_ERROR;
996 }
997
998
999
1000
transformPrimary(msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,double transform[3][3])1001 msym_error_t transformPrimary(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double transform[3][3]){
1002 msym_error_t ret = MSYM_SUCCESS;
1003 if(primary != NULL){
1004 double z[3] = {0.0, 0.0, 1.0};
1005 malign(primary->v,z,transform);
1006 for(msym_symmetry_operation_t *sop = sops; sop < (sops + sopsl); sop++){
1007 mvmul(sop->v,transform,sop->v);
1008 }
1009 vcopy(z,primary->v); // get rid of small errors
1010 } else {
1011 msymSetErrorDetails("Point group has no primary axis for transformation");
1012 ret = MSYM_POINT_GROUP_ERROR;
1013 }
1014 return ret;
1015 }
1016
transformSecondary(msym_point_group_type_t type,msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,double transform[3][3])1017 msym_error_t transformSecondary(msym_point_group_type_t type, msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double transform[3][3]){
1018 msym_error_t ret = MSYM_SUCCESS;
1019 double axis[3], x[3] = {1.0,0.0,0.0};
1020
1021 switch(type){
1022 case (MSYM_POINT_GROUP_TYPE_Cnv) :
1023 if(MSYM_SUCCESS != (ret = findSecondaryAxisSigma(primary, sopsl, sops, thresholds, axis))) goto err;
1024 break;
1025 case (MSYM_POINT_GROUP_TYPE_O) :
1026 case (MSYM_POINT_GROUP_TYPE_Oh) :
1027 if(MSYM_SUCCESS != (ret = findSecondaryAxisC4(primary, sopsl, sops, thresholds, axis))) goto err;
1028 break;
1029 case (MSYM_POINT_GROUP_TYPE_Dn) :
1030 case (MSYM_POINT_GROUP_TYPE_Dnh) :
1031 case (MSYM_POINT_GROUP_TYPE_Dnd) :
1032 case (MSYM_POINT_GROUP_TYPE_T) :
1033 case (MSYM_POINT_GROUP_TYPE_Td) :
1034 case (MSYM_POINT_GROUP_TYPE_Th) :
1035 if(MSYM_SUCCESS != (ret = findSecondaryAxisC2(primary, sopsl, sops, thresholds, axis))) goto err;
1036 break;
1037 case (MSYM_POINT_GROUP_TYPE_I) :
1038 case (MSYM_POINT_GROUP_TYPE_Ih) :
1039 if(MSYM_SUCCESS != (ret = findSecondaryAxisC2C5(primary, sopsl, sops, thresholds, axis))) goto err;
1040 break;
1041 default :
1042 msymSetErrorDetails("Unknown point group when determining secondary axis");
1043 ret = MSYM_POINT_GROUP_ERROR;
1044 goto err;
1045 }
1046
1047 double m[3][3];
1048
1049 malign(axis, x, m);
1050
1051 for(msym_symmetry_operation_t *sop = sops; sop < (sops + sopsl); sop++){
1052 mvmul(sop->v,m,sop->v);
1053 }
1054 mmmul(m,transform,transform);
1055
1056 err:
1057 return ret;
1058 }
1059
1060
1061
1062 /* For point groups where we use a perpendicular reflection plane to indicate direction.
1063 We use a vector where the the xy-plane and reflection plane cross
1064 */
findSecondaryAxisSigma(msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,double r[3])1065 msym_error_t findSecondaryAxisSigma(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double r[3]){
1066 msym_error_t ret = MSYM_SUCCESS;
1067 msym_symmetry_operation_t *sop = NULL;
1068 for(sop = sops; sop < (sops + sopsl); sop++){
1069 if(sop->type == REFLECTION){
1070 vcross(sop->v, primary->v, r);
1071 vnorm(r);
1072 break;
1073 }
1074 }
1075 if(sop == (sops + sopsl)){
1076 msymSetErrorDetails("Can't find secondary reflection plane when symmetrizing point group");
1077 ret = MSYM_POINT_GROUP_ERROR;
1078 goto err;
1079 }
1080 err:
1081 return ret;
1082 }
1083
1084 /* For point groups where we use a perpendicular C2 axis to indicate direction.
1085 Adjusted to make sure it's perfectly in the xy-plane.
1086 */
findSecondaryAxisC2(msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,double r[3])1087 msym_error_t findSecondaryAxisC2(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double r[3]){
1088 msym_error_t ret = MSYM_SUCCESS;
1089 msym_symmetry_operation_t *sop = NULL;
1090 for(msym_symmetry_operation_t *sop = sops; sop < (sops + sopsl); sop++){
1091 //Choose a C2 perpendicular to the primary axis, it'll make things relatively easy
1092 if(sop != primary && sop->type == PROPER_ROTATION && sop->order == 2 && vperpendicular(sop->v, primary->v,thresholds->angle)){
1093 vproj_plane(sop->v, primary->v, r);
1094 vnorm(r);
1095 break;
1096 }
1097 }
1098 if(sop == (sops + sopsl)){
1099 msymSetErrorDetails("Can't find secondary C2 axis when symmetrizing point group");
1100 ret = MSYM_POINT_GROUP_ERROR;
1101 goto err;
1102 }
1103 err:
1104 return ret;
1105 }
1106
1107
findSecondaryAxisC2C5(msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,double r[3])1108 msym_error_t findSecondaryAxisC2C5(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double r[3]){
1109 msym_error_t ret = MSYM_SUCCESS;
1110 msym_symmetry_operation_t *c2[2], *c5 = NULL;
1111 int c2i = 0;
1112
1113 for(msym_symmetry_operation_t *sop = sops; sop < (sops + sopsl) && (c5 == NULL || c2i < 2); sop++){
1114 if(vperpendicular(sop->v, primary->v,thresholds->angle)){
1115 if(sop->type == PROPER_ROTATION && sop->order == 2){
1116 //printf("Found perpendicular C2\n");
1117 c2[c2i++] = sop;
1118 } else if (sop->type == PROPER_ROTATION && sop->order == 5){
1119 //printf("Found perpendicular C5\n");
1120 c5 = sop;
1121 }
1122 }
1123 }
1124
1125 if(c5 == NULL || c2i < 2) {
1126 msymSetErrorDetails("Can't find secondary C2 axis when symmetrizing point group: (%s %s)",c5 == NULL ? "C5 axis missing" : "", c2i < 2 ? "C2 axis missing" : "");
1127 ret = MSYM_POINT_GROUP_ERROR;
1128 goto err;
1129 }
1130
1131 if(fabs(vdot(c5->v, c2[0]->v)) > fabs(vdot(c5->v, c2[1]->v))){
1132 vproj_plane(c2[0]->v, primary->v, r);
1133 } else {
1134 vproj_plane(c2[1]->v, primary->v, r);
1135 }
1136
1137 err:
1138 return ret;
1139 }
1140
1141
1142
1143 // Same as C2
findSecondaryAxisC4(msym_symmetry_operation_t * primary,int sopsl,msym_symmetry_operation_t sops[sopsl],msym_thresholds_t * thresholds,double r[3])1144 msym_error_t findSecondaryAxisC4(msym_symmetry_operation_t *primary, int sopsl, msym_symmetry_operation_t sops[sopsl], msym_thresholds_t *thresholds, double r[3]){
1145 msym_error_t ret = MSYM_SUCCESS;
1146 msym_symmetry_operation_t *sop = NULL;
1147 for(sop = sops; sop < (sops + sopsl); sop++){
1148 //Choose a C2 perpendicular to the primary axis, it'll make things relatively easy
1149 if(sop != primary && sop->type == PROPER_ROTATION && sop->order == 4 && vperpendicular(sop->v, primary->v,thresholds->angle)){
1150 vproj_plane(sop->v, primary->v, r);
1151 vnorm(r);
1152 break;
1153 }
1154 }
1155 if(sop == (sops + sopsl)){
1156 msymSetErrorDetails("Can't find secondary C4 axis when symmetrizing point group");
1157 ret = MSYM_POINT_GROUP_ERROR;
1158 goto err;
1159 }
1160 err:
1161 return ret;
1162 }
1163
generateSymmetryOperationsImpliedRot(int sopsl,msym_symmetry_operation_t sops[sopsl],int order,msym_thresholds_t * thresholds,int * osopsl)1164 msym_error_t generateSymmetryOperationsImpliedRot(int sopsl, msym_symmetry_operation_t sops[sopsl], int order, msym_thresholds_t *thresholds, int *osopsl){
1165 int isopsl = sopsl;
1166 for(msym_symmetry_operation_t *sopi = sops; sopi < (sops + sopsl) && isopsl < order; sopi++){
1167 if(sopi->type == PROPER_ROTATION){
1168 for(msym_symmetry_operation_t *sopj = sops; sopj < (sops + sopsl); sopj++){
1169 int istype = (sopj->type == REFLECTION || sopj->type == IMPROPER_ROTATION || (sopj->type == PROPER_ROTATION));
1170 if(sopi != sopj && istype && !vparallel(sopi->v, sopj->v,thresholds->angle)){
1171 copySymmetryOperation(&(sops[isopsl]), sopj);
1172 applySymmetryOperation(sopi,sops[isopsl].v,sops[isopsl].v);
1173 isopsl += !findSymmetryOperation(&(sops[isopsl]), sops, isopsl,thresholds);
1174 if(isopsl > order) goto err;
1175 }
1176 }
1177 }
1178 }
1179 *osopsl = isopsl;
1180 return MSYM_SUCCESS;
1181 err:
1182 msymSetErrorDetails("Generation of implied symmetry operations by rotation resulted in more operations than point group order");
1183 return MSYM_POINT_GROUP_ERROR;
1184 }
1185
generateSymmetryOperations(msym_point_group_type_t type,int n,int order,msym_symmetry_operation_t ** osops)1186 msym_error_t generateSymmetryOperations(msym_point_group_type_t type, int n, int order, msym_symmetry_operation_t **osops){
1187 msym_error_t ret = MSYM_SUCCESS;
1188 msym_symmetry_operation_t *sops = calloc(order, sizeof(msym_symmetry_operation_t));
1189 sops[0].cla = 0;
1190 sops[0].type = IDENTITY;
1191 sops[0].power = 1;
1192 sops[0].order = 1;
1193 sops[0].orientation = NONE;
1194 int cla = 1, gsopsl = 1;
1195
1196 const struct _fmap {
1197 msym_point_group_type_t type;
1198 msym_error_t (*f)(int, int, msym_symmetry_operation_t *, int *, int *);
1199 } fmap[18] = {
1200
1201 [ 0] = {MSYM_POINT_GROUP_TYPE_Ci, generateSymmetryOperationsCi},
1202 [ 1] = {MSYM_POINT_GROUP_TYPE_Cs, generateSymmetryOperationsCs},
1203 [ 2] = {MSYM_POINT_GROUP_TYPE_Cn, generateSymmetryOperationsCn},
1204 [ 3] = {MSYM_POINT_GROUP_TYPE_Cnh, generateSymmetryOperationsCnh},
1205 [ 4] = {MSYM_POINT_GROUP_TYPE_Cnv, generateSymmetryOperationsCnv},
1206 [ 5] = {MSYM_POINT_GROUP_TYPE_Dn, generateSymmetryOperationsDn},
1207 [ 6] = {MSYM_POINT_GROUP_TYPE_Dnh, generateSymmetryOperationsDnh},
1208 [ 7] = {MSYM_POINT_GROUP_TYPE_Dnd, generateSymmetryOperationsDnd},
1209 [ 8] = {MSYM_POINT_GROUP_TYPE_Sn, generateSymmetryOperationsSn},
1210 [ 9] = {MSYM_POINT_GROUP_TYPE_T, generateSymmetryOperationsT},
1211 [10] = {MSYM_POINT_GROUP_TYPE_Td, generateSymmetryOperationsTd},
1212 [11] = {MSYM_POINT_GROUP_TYPE_Th, generateSymmetryOperationsTh},
1213 [12] = {MSYM_POINT_GROUP_TYPE_O, generateSymmetryOperationsO},
1214 [13] = {MSYM_POINT_GROUP_TYPE_Oh, generateSymmetryOperationsOh},
1215 [14] = {MSYM_POINT_GROUP_TYPE_I, generateSymmetryOperationsI},
1216 [15] = {MSYM_POINT_GROUP_TYPE_Ih, generateSymmetryOperationsIh},
1217 [16] = {MSYM_POINT_GROUP_TYPE_K, generateSymmetryOperationsUnknown},
1218 [17] = {MSYM_POINT_GROUP_TYPE_Kh, generateSymmetryOperationsUnknown}
1219 };
1220
1221 int fi, fil = sizeof(fmap)/sizeof(fmap[0]);
1222 for(fi = 0; fi < fil;fi++){
1223 if(fmap[fi].type == type) {
1224 if(MSYM_SUCCESS != (ret = fmap[fi].f(n,order,sops,&gsopsl,&cla))) goto err;
1225 break;
1226 }
1227 }
1228
1229 if(fi == fil){
1230 msymSetErrorDetails("Unknown point group when generating symmetry operations");
1231 ret = MSYM_POINT_GROUP_ERROR;
1232 goto err;
1233 }
1234
1235 if(gsopsl != order){
1236 msymSetErrorDetails("Generated incorrect number of symmetry operations %d != %d",gsopsl,order);
1237 ret = MSYM_INVALID_POINT_GROUP;
1238 goto err;
1239 }
1240
1241 for(int i = 0;i < order;i++) printSymmetryOperation(&sops[i]);
1242
1243 *osops = sops;
1244 return ret;
1245 err:
1246 free(sops);
1247 return ret;
1248
1249 }
1250
generateSymmetryOperationsUnknown(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1251 msym_error_t generateSymmetryOperationsUnknown(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1252 msymSetErrorDetails("Generating symmetry operations for unknown point group");
1253 return MSYM_POINT_GROUP_ERROR;
1254 }
generateSymmetryOperationsSn(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1255 msym_error_t generateSymmetryOperationsSn(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1256 msym_error_t ret = MSYM_SUCCESS;
1257 int k = *pk, cla = *pcla, m = (n << (n & 1));
1258 double z[3] = {0.0,0.0,1.0};
1259 msym_symmetry_operation_t sn = {.type = IMPROPER_ROTATION, .order = n, .power = 1, .orientation = HORIZONTAL};
1260 vcopy(z,sn.v);
1261 if(k + m - 1 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating S%d symmetry operations",n); goto err;}
1262 for(int i = 1;i <= m >> 1;i++){
1263 int index = k + ((i-1) << 1);
1264 symopPow(&sn, i, &sops[index]);
1265 sops[index].cla = cla + i - 1;
1266 clean_debug_printf("i = %d m = %d index = %d ",i,m,index);
1267 printSymmetryOperation(&sops[index]);
1268 }
1269 /*
1270 int ri = k + (((m >> 1)-1) << 1);
1271 sops[ri].cla = cla + (m >> 1) - 1;
1272 sops[ri].power = 1;
1273 sops[ri].p.orientation = HORIZONTAL;
1274 sops[ri].type = REFLECTION;
1275 vcopy(z,sops[ri].v);
1276 clean_debug_printf("replacing symmetry operation %d\n",ri);
1277 printSymmetryOperation(&sops[ri]);*/
1278
1279 for(int i = 1;i < m >> 1;i++){
1280 int index = k + 1 + ((i-1) << 1);
1281 symopPow(&sn, m-i, &sops[index]);
1282 sops[index].cla = cla + i - 1;
1283 clean_debug_printf("i = %d m = %d index = %d ",i,m,index);
1284 printSymmetryOperation(&sops[index]);
1285
1286 }
1287 k += m - 1;
1288 cla += m >> 1;
1289
1290 clean_debug_printf("------ Sn %d operations %d classes------\n",k-*pk, cla-*pcla);
1291 *pk = k; *pcla = cla;
1292
1293 return ret;
1294 err:
1295 return ret;
1296 }
1297
generateSymmetryOperationsCs(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1298 msym_error_t generateSymmetryOperationsCs(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1299 msym_error_t ret = MSYM_SUCCESS;
1300 int k = *pk, cla = *pcla;
1301
1302 if(k > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating Cs symmetry operations"); goto err;}
1303
1304
1305 msym_symmetry_operation_t sigma = {.type = REFLECTION, .orientation = HORIZONTAL, .order = 1, .power = 1, .cla = cla, .v = {0,0,1}};
1306
1307 memcpy(&(sops[k]), &sigma, sizeof(msym_symmetry_operation_t));
1308
1309 k++;
1310 cla++;
1311
1312 clean_debug_printf("------ Cs %d operations %d classes------\n",k-*pk, cla-*pcla);
1313 *pk = k; *pcla = cla;
1314
1315 return ret;
1316 err:
1317 return ret;
1318
1319 }
1320
generateSymmetryOperationsCi(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1321 msym_error_t generateSymmetryOperationsCi(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1322 msym_error_t ret = MSYM_SUCCESS;
1323 int k = *pk, cla = *pcla;
1324
1325 if(k > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating Ci symmetry operations"); goto err;}
1326
1327
1328 msym_symmetry_operation_t inv = {.type = INVERSION, .orientation = NONE, .order = 1, .power = 1, .cla = cla, .v = {0,0,0}};
1329
1330 memcpy(&(sops[k]), &inv, sizeof(msym_symmetry_operation_t));
1331
1332 k++;
1333 cla++;
1334
1335 clean_debug_printf("------ Ci %d operations %d classes------\n",k-*pk, cla-*pcla);
1336 *pk = k; *pcla = cla;
1337
1338 return ret;
1339 err:
1340 return ret;
1341
1342 }
1343
generateSymmetryOperationsCn(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1344 msym_error_t generateSymmetryOperationsCn(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1345 msym_error_t ret = MSYM_SUCCESS;
1346 int k = *pk, cla = *pcla;
1347 double z[3] = {0.0,0.0,1.0};
1348 msym_symmetry_operation_t cn = {.type = PROPER_ROTATION, .order = n, .power = 1, .orientation = HORIZONTAL};
1349 if(k + n - 1 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating C%d symmetry operations",n); goto err;}
1350 vcopy(z,cn.v);
1351
1352 for(int i = 1;i <= (n >> 1);i++){
1353 int index = k + (i << 1) - 2;
1354 symopPow(&cn, i, &sops[index]);
1355 sops[index].cla = cla + (index >> 1);
1356 clean_debug_printf("i = %d index = %d ",i,index);
1357 printSymmetryOperation(&sops[index]);
1358 }
1359
1360 for(int i = 1;i < (n >> 1) + (n & 1);i++){
1361 int index = k + (i << 1) - 1;
1362 symopPow(&cn, n-i, &sops[index]);
1363 sops[index].cla = cla + ((index - 1) >> 1);
1364 clean_debug_printf("i = %d index = %d ",i,index);
1365 printSymmetryOperation(&sops[index]);
1366 }
1367
1368 k += n - 1;
1369 cla += n >> 1;
1370
1371 clean_debug_printf("------ Cn %d operations %d classes------\n",k-*pk, cla-*pcla);
1372 *pk = k; *pcla = cla;
1373
1374 return ret;
1375 err:
1376 return ret;
1377 }
1378
generateSymmetryOperationsCnh(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1379 msym_error_t generateSymmetryOperationsCnh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1380 msym_error_t ret = MSYM_SUCCESS;
1381 int k = *pk, cla = *pcla, s = 0;
1382 double z[3] = {0.0,0.0,1.0};
1383 msym_symmetry_operation_t cn = {.type = PROPER_ROTATION, .order = n, .power = 1, .orientation = HORIZONTAL};
1384 msym_symmetry_operation_t sn = {.type = IMPROPER_ROTATION, .order = n, .power = 1, .orientation = HORIZONTAL};
1385 if(k + (n << 1) - 1 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating C%dh symmetry operations",n); goto err;}
1386 vcopy(z,cn.v); vcopy(z,sn.v);
1387 for(s = n;s % 2 == 0;s = s >> 1){
1388 cn.order = s;
1389 for(int i = 1;i <= s >> 1;i += 2){
1390 int index = k + ((i >> 1) << 1);
1391 symopPow(&cn, i, &sops[index]);
1392 sops[index].cla = cla + (i >> 1);
1393 printSymmetryOperation(&sops[index]);
1394 }
1395 for(int i = 1;i < s >> 1;i += 2){
1396 int index = k + 1 + ((i >> 1) << 1);
1397 symopPow(&cn, s-i, &sops[index]);
1398 sops[index].cla = cla + (i >> 1);
1399 printSymmetryOperation(&sops[index]);
1400 }
1401
1402 k += (s >> 1);
1403 cla += (s >> 2) + ((s >> 1) & 1);
1404
1405 sn.order = s;
1406 for(int i = 1;i <= s >> 1;i += 2){
1407 int index = k + ((i >> 1) << 1);
1408 symopPow(&sn, i, &sops[index]);
1409 sops[index].cla = cla + (i >> 1);
1410 printSymmetryOperation(&sops[index]);
1411 }
1412
1413 for(int i = 1;i < s >> 1;i += 2){
1414 int index = k + 1 + ((i >> 1) << 1);
1415 symopPow(&sn, s-i, &sops[index]);
1416 sops[index].cla = cla + (i >> 1);
1417 printSymmetryOperation(&sops[index]);
1418 }
1419
1420 k += (s >> 1);
1421 cla += (s >> 2) + ((s >> 1) & 1);
1422
1423 }
1424
1425 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsSn(s,l,sops,&k,&cla))) goto err;
1426
1427 clean_debug_printf("------ Cnh %d operations %d classes------\n",k-*pk, cla-*pcla);
1428 *pk = k; *pcla = cla;
1429
1430 return ret;
1431 err:
1432 return ret;
1433
1434 }
1435
generateSymmetryOperationsCnv(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1436 msym_error_t generateSymmetryOperationsCnv(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1437 msym_error_t ret = MSYM_SUCCESS;
1438 int k = *pk, cla = *pcla;
1439
1440 if(n == 0 && l == 2){ // normal c0v
1441 if(k + 1 > l){
1442 ret = MSYM_POINT_GROUP_ERROR;
1443 msymSetErrorDetails("Too many operations when generating C%dv symmetry operations",n);
1444 goto err;
1445 }
1446 msym_symmetry_operation_t c0 = {.type = PROPER_ROTATION, .order = n, .power = 1, .orientation = HORIZONTAL, .cla = cla, .v = {0,0,1}};
1447 memcpy(&sops[k], &c0, sizeof(msym_symmetry_operation_t));
1448 k += 1;
1449 cla += 1;
1450
1451 } else if(n == 0){ // c0v represented by a subclass cnv where n is even
1452 int n0 = l/2;
1453 double z[3] = {0,0,1};
1454 if(n0 & 1){
1455 ret = MSYM_POINT_GROUP_ERROR;
1456 msymSetErrorDetails("Cannot generate an odd representative (C%dv) of C0v",n0);
1457 goto err;
1458 }
1459
1460 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsCn(n0,l,sops,&k,&cla))) goto err;
1461
1462 if(k + n0 > l){
1463 ret = MSYM_POINT_GROUP_ERROR;
1464 msymSetErrorDetails("Too many operations when generating D0h (D%dh) symmetry operations",n0);
1465 goto err;
1466 }
1467
1468 // Can't use generateReflectionPlanes they'll generate vertical and dihedral
1469 msym_symmetry_operation_t sigma = {.type = REFLECTION, .power = 1, .order = 1, .orientation = VERTICAL, .v = {0,1,0}};
1470 for(int i = 0;i < n0;i++){
1471 int index = k+i;
1472 memcpy(&(sops[index]), &sigma, sizeof(msym_symmetry_operation_t));
1473 vrotate(i*M_PI/n0, sigma.v, z, sops[index].v);
1474 sops[index].cla = cla;
1475 }
1476
1477 k += n0;
1478 cla += 1;
1479
1480 clean_debug_printf("\n------ C0v %d operations %d classes------\n",k-*pk, cla-*pcla);
1481 } else {
1482 if(k + (n << 1) - 1 > l){
1483 ret = MSYM_POINT_GROUP_ERROR;
1484 msymSetErrorDetails("Too many operations when generating C%dv symmetry operations",n);
1485 goto err;
1486 }
1487 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsCn(n,l,sops,&k,&cla))) goto err;
1488 if(MSYM_SUCCESS != (ret = generateReflectionPlanes(n,l,sops,&k,&cla))) goto err;
1489 clean_debug_printf("------ Cnv %d operations %d classes------\n",k-*pk, cla-*pcla);
1490
1491 }
1492
1493 *pk = k; *pcla = cla;
1494
1495 return ret;
1496 err:
1497 return ret;
1498 }
1499
generateSymmetryOperationsDn(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1500 msym_error_t generateSymmetryOperationsDn(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1501 msym_error_t ret = MSYM_SUCCESS;
1502 int k = *pk, cla = *pcla;
1503 if(k + (n << 1) - 1 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating D%d symmetry operations",n); goto err;}
1504 //k += (n << 1) - 1;
1505 //cla = sops[k-1].cla + 1;
1506 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsCn(n,l,sops,&k,&cla))) goto err;
1507 //k += n;
1508 //cla = sops[k-1].cla + 1;
1509 if(MSYM_SUCCESS != (ret = generateC2Axes(n,l,sops,&k,&cla))) goto err;
1510 //k += n;
1511 //cla = sops[k-1].cla + 1;
1512
1513 clean_debug_printf("\n------ Dn %d operations %d classes------\n",k-*pk, cla-*pcla);
1514 *pk = k; *pcla = cla;
1515
1516 return ret;
1517 err:
1518 return ret;
1519 }
1520
generateSymmetryOperationsDnh(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1521 msym_error_t generateSymmetryOperationsDnh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1522 msym_error_t ret = MSYM_SUCCESS;
1523 int k = *pk, cla = *pcla;
1524
1525 if(n == 0 && l == 4){ // standard d0h
1526 if(k + 3 > l){
1527 ret = MSYM_POINT_GROUP_ERROR;
1528 msymSetErrorDetails("Too many operations when generating D%dh symmetry operations",n);
1529 goto err;
1530 }
1531 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsCnv(n,2,sops,&k,&cla))) goto err;
1532 msym_symmetry_operation_t sigma = {.type = REFLECTION, .order = 1, .power = 1, .orientation = HORIZONTAL, .cla = cla, .v = {0,0,1}};
1533 msym_symmetry_operation_t inversion = {.type = INVERSION, .order = 1, .power = 1, .orientation = NONE, .cla = cla+1, .v = {0,0,1}};
1534 memcpy(&sops[k], &sigma, sizeof(msym_symmetry_operation_t));
1535 memcpy(&sops[k+1], &inversion, sizeof(msym_symmetry_operation_t));
1536 k += 2;
1537 cla += 2;
1538 }
1539 else if(n == 0){ // d0h represented by a subclass cnv where n is even
1540 int n0 = l/4;
1541 double z[3] = {0,0,1};
1542 if(n0 & 1){
1543 ret = MSYM_POINT_GROUP_ERROR;
1544 msymSetErrorDetails("Cannot generate an odd representative (D%dh) of D0h",n0);
1545 goto err;
1546 }
1547
1548 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsCnh(n0,l,sops,&k,&cla))) goto err;
1549
1550 if(k + (n0 << 1) > l){
1551 ret = MSYM_POINT_GROUP_ERROR;
1552 msymSetErrorDetails("Too many operations when generating D0h (D%dh) symmetry operations",n0);
1553 goto err;
1554 }
1555
1556 // Can't use generateReflectionPlanes/generateC2Axes they'll generate vertical and dihedral
1557 msym_symmetry_operation_t c2 = {.type = PROPER_ROTATION, .power = 1, .order = 2, .orientation = VERTICAL, .v = {1,0,0}};
1558 msym_symmetry_operation_t sigma = {.type = REFLECTION, .power = 1, .order = 1, .orientation = VERTICAL, .v = {0,1,0}};
1559 for(int i = 0;i < n0;i++){
1560 int index = k+i;
1561 memcpy(&(sops[index]), &sigma, sizeof(msym_symmetry_operation_t));
1562 vrotate(i*M_PI/n0, sigma.v, z, sops[index].v);
1563 sops[index].cla = cla;
1564 clean_debug_printf("generated sigma at %d\n",index);
1565 index += n0;
1566 memcpy(&(sops[index]), &c2, sizeof(msym_symmetry_operation_t));
1567 vrotate(i*M_PI/n0, c2.v, z, sops[index].v);
1568 sops[index].cla = cla+1;
1569 }
1570
1571 k += n0 << 1;
1572 cla += 2;
1573
1574 clean_debug_printf("\n------ D0h %d operations %d classes------\n",k-*pk, cla-*pcla);
1575 } else {
1576
1577 if(k + (n << 2) - 1 > l){
1578 ret = MSYM_POINT_GROUP_ERROR;
1579 msymSetErrorDetails("Too many operations when generating D%dh symmetry operations",n);
1580 goto err;
1581 }
1582 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsCnh(n,l,sops,&k,&cla))) goto err;
1583 if(MSYM_SUCCESS != (ret = generateReflectionPlanes(n,l,sops,&k,&cla))) goto err;
1584 if(MSYM_SUCCESS != (ret = generateC2Axes(n,l,sops,&k,&cla))) goto err;
1585 clean_debug_printf("\n------ Dnh %d operations %d classes------\n",k-*pk, cla-*pcla);
1586
1587 }
1588
1589 *pk = k; *pcla = cla;
1590
1591 return ret;
1592 err:
1593 return ret;
1594 }
1595
generateSymmetryOperationsDnd(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1596 msym_error_t generateSymmetryOperationsDnd(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1597 msym_error_t ret = MSYM_SUCCESS;
1598 int k = *pk, cla = *pcla;
1599 double x[3] = {1.0,0.0,0.0}, y[3] = {0.0,1.0,0.0}, z[3] = {0.0,0.0,1.0};
1600 msym_symmetry_operation_t sigma = {.type = REFLECTION, .orientation = DIHEDRAL, .power = 1};
1601 msym_symmetry_operation_t c2 = {.type = PROPER_ROTATION, .order = 2, .power = 1, .orientation = VERTICAL};
1602
1603 if(k + (n << 2) - 1 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating D%dd symmetry operations",n); goto err;}
1604
1605 vcopy(x,c2.v); vrotate(M_PI_2/n, y, z, sigma.v);
1606
1607 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsSn((n << 1),l,sops, &k, &cla))) goto err;
1608
1609 for(int i = 0;i < n;i++){
1610 memcpy(&(sops[k+i]), &sigma, sizeof(msym_symmetry_operation_t));
1611 vrotate(i*M_PI/n, sigma.v, z, sops[k+i].v);
1612 sops[k+i].cla = cla;
1613 }
1614 k += n;
1615 cla += 1;
1616
1617 for(int i = 0;i < n;i++){
1618 memcpy(&(sops[k+i]), &c2, sizeof(msym_symmetry_operation_t));
1619 vrotate(i*M_PI/n, c2.v, z, sops[k+i].v);
1620 sops[k+i].cla = cla;
1621 }
1622 k += n;
1623 cla += 1;
1624
1625
1626 clean_debug_printf("\n------ Dnd %d operations %d classes------\n",k-*pk, cla-*pcla);
1627 *pk = k; *pcla = cla;
1628
1629 return ret;
1630 err:
1631 return ret;
1632 }
1633
generateSymmetryOperationsT(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1634 msym_error_t generateSymmetryOperationsT(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1635 msym_error_t ret = MSYM_SUCCESS;
1636 int k = *pk, cla = *pcla;
1637
1638 msym_symmetry_operation_t c2[1] = {
1639 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla, .orientation = HORIZONTAL}
1640 };
1641
1642 msym_symmetry_operation_t c3[2] = {
1643 [0] = {.type = PROPER_ROTATION, .order = 3, .power = 1, .cla = cla+1, .orientation = NONE},
1644 [1] = {.type = PROPER_ROTATION, .order = 3, .power = 2, .cla = cla+1, .orientation = NONE}
1645 };
1646
1647 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsTetrahedral(l,sops, 1, c2, 0, NULL, 2, c3, &k))) goto err;
1648
1649 cla += 2;
1650
1651 clean_debug_printf("------ T %d operations %d classes------\n",k-*pk, cla-*pcla);
1652 *pk = k; *pcla = cla;
1653
1654 return ret;
1655 err:
1656 return ret;
1657 }
1658
generateSymmetryOperationsTd(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1659 msym_error_t generateSymmetryOperationsTd(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1660 msym_error_t ret = MSYM_SUCCESS;
1661 int k = *pk, cla = *pcla;
1662
1663 msym_symmetry_operation_t c2[3] = {
1664 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla, .orientation = HORIZONTAL},
1665 [1] = {.type = IMPROPER_ROTATION, .order = 4, .power = 1, .cla = cla+1, .orientation = HORIZONTAL},
1666 [2] = {.type = IMPROPER_ROTATION, .order = 4, .power = 3, .cla = cla+1, .orientation = HORIZONTAL}
1667 };
1668
1669 msym_symmetry_operation_t cs[4] = {
1670 [0] = {.type = REFLECTION, .order = 1, .power = 1, .cla = cla+2, .orientation = DIHEDRAL},
1671 };
1672
1673 msym_symmetry_operation_t c3[2] = {
1674 [0] = {.type = PROPER_ROTATION, .order = 3, .power = 1, .cla = cla+3, .orientation = NONE},
1675 [1] = {.type = PROPER_ROTATION, .order = 3, .power = 2, .cla = cla+3, .orientation = NONE}
1676 };
1677
1678 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsTetrahedral(l,sops, 3, c2, 1, cs, 2, c3, &k))) goto err;
1679
1680 cla += 4;
1681
1682 clean_debug_printf("------ Td %d operations %d classes------\n",k-*pk, cla-*pcla);
1683 *pk = k; *pcla = cla;
1684
1685 return ret;
1686 err:
1687 return ret;
1688 }
1689
generateSymmetryOperationsTh(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1690 msym_error_t generateSymmetryOperationsTh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1691 msym_error_t ret = MSYM_SUCCESS;
1692 int k = *pk, cla = *pcla;
1693
1694 msym_symmetry_operation_t c2[2] = {
1695 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla, .orientation = HORIZONTAL},
1696 [1] = {.type = REFLECTION, .order = 1, .power = 1, .cla = cla+1, .orientation = HORIZONTAL}
1697 };
1698
1699 msym_symmetry_operation_t c3[4] = {
1700 [0] = {.type = PROPER_ROTATION, .order = 3, .power = 1, .cla = cla+2, .orientation = NONE},
1701 [1] = {.type = PROPER_ROTATION, .order = 3, .power = 2, .cla = cla+2, .orientation = NONE},
1702 [2] = {.type = IMPROPER_ROTATION, .order = 6, .power = 1, .cla = cla+3, .orientation = NONE},
1703 [3] = {.type = IMPROPER_ROTATION, .order = 6, .power = 5, .cla = cla+3, .orientation = NONE}
1704 };
1705
1706 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsTetrahedral(l,sops, 2, c2, 0, NULL, 4, c3, &k))) goto err;
1707
1708 if(k - 1 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating Th operations %d >= %d",k, l); goto err;}
1709
1710 sops[k].type = INVERSION;
1711 sops[k].power = 1;
1712 sops[k].order = 1;
1713 sops[k].orientation = NONE;
1714 sops[k].cla = cla+4;
1715 k++;
1716
1717 cla += 5;
1718
1719 clean_debug_printf("------ Th %d operations %d classes------\n",k-*pk, cla-*pcla);
1720 *pk = k; *pcla = cla;
1721
1722 return ret;
1723 err:
1724 return ret;
1725 }
1726
generateSymmetryOperationsO(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1727 msym_error_t generateSymmetryOperationsO(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1728 msym_error_t ret = MSYM_SUCCESS;
1729 int k = *pk, cla = *pcla;
1730
1731 msym_symmetry_operation_t c2[1] = {
1732 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla, .orientation = VERTICAL}
1733 };
1734
1735 msym_symmetry_operation_t c3[2] = {
1736 [0] = {.type = PROPER_ROTATION, .order = 3, .power = 1, .cla = cla+1, .orientation = NONE},
1737 [1] = {.type = PROPER_ROTATION, .order = 3, .power = 2, .cla = cla+1, .orientation = NONE}
1738 };
1739
1740 msym_symmetry_operation_t c4[4] = {
1741 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla+2, .orientation = HORIZONTAL},
1742 [1] = {.type = PROPER_ROTATION, .order = 4, .power = 1, .cla = cla+3, .orientation = HORIZONTAL},
1743 [2] = {.type = PROPER_ROTATION, .order = 4, .power = 3, .cla = cla+3, .orientation = HORIZONTAL}
1744 };
1745
1746 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsOctahedral(l,sops, 1, c2, 2, c3, 3, c4, &k))) goto err;
1747
1748 cla += 4;
1749
1750 clean_debug_printf("------ O %d operations %d classes------\n",k-*pk, cla-*pcla);
1751 *pk = k; *pcla = cla;
1752
1753 return ret;
1754 err:
1755 return ret;
1756 }
1757
generateSymmetryOperationsOh(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1758 msym_error_t generateSymmetryOperationsOh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1759 msym_error_t ret = MSYM_SUCCESS;
1760 int k = *pk, cla = *pcla;
1761
1762 msym_symmetry_operation_t c2[2] = {
1763 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla, .orientation = VERTICAL},
1764 [1] = {.type = REFLECTION, .order = 1, .power = 1, .cla = cla+1, .orientation = DIHEDRAL}
1765 };
1766
1767 msym_symmetry_operation_t c3[4] = {
1768 [0] = {.type = PROPER_ROTATION, .order = 3, .power = 1, .cla = cla+2, .orientation = NONE},
1769 [1] = {.type = PROPER_ROTATION, .order = 3, .power = 2, .cla = cla+2, .orientation = NONE},
1770 [2] = {.type = IMPROPER_ROTATION, .order = 6, .power = 1, .cla = cla+3, .orientation = NONE},
1771 [3] = {.type = IMPROPER_ROTATION, .order = 6, .power = 5, .cla = cla+3, .orientation = NONE}
1772 };
1773
1774 msym_symmetry_operation_t c4[6] = {
1775 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla+4, .orientation = HORIZONTAL},
1776 [1] = {.type = PROPER_ROTATION, .order = 4, .power = 1, .cla = cla+5, .orientation = HORIZONTAL},
1777 [2] = {.type = PROPER_ROTATION, .order = 4, .power = 3, .cla = cla+5, .orientation = HORIZONTAL},
1778 [3] = {.type = IMPROPER_ROTATION, .order = 4, .power = 1, .cla = cla+6, .orientation = HORIZONTAL},
1779 [4] = {.type = IMPROPER_ROTATION, .order = 4, .power = 3, .cla = cla+6, .orientation = HORIZONTAL},
1780 [5] = {.type = REFLECTION, .order = 1, .power = 1, .cla = cla+7, .orientation = HORIZONTAL}
1781 };
1782
1783 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsOctahedral(l,sops, 2, c2, 4, c3, 6, c4, &k))) goto err;
1784
1785 if(k - 1 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating Oh operations %d >= %d",k, l); goto err;}
1786
1787 sops[k].type = INVERSION;
1788 sops[k].power = 1;
1789 sops[k].order = 1;
1790 sops[k].orientation = NONE;
1791 sops[k].cla = cla+8;
1792 k++;
1793
1794 cla += 9;
1795
1796 clean_debug_printf("------ O %d operations %d classes------\n",k-*pk, cla-*pcla);
1797 *pk = k; *pcla = cla;
1798
1799 return ret;
1800 err:
1801 return ret;
1802 }
1803
generateSymmetryOperationsI(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1804 msym_error_t generateSymmetryOperationsI(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1805 msym_error_t ret = MSYM_SUCCESS;
1806 int k = *pk, cla = *pcla;
1807
1808 msym_symmetry_operation_t c2[1] = {
1809 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla, .orientation = NONE}
1810 };
1811
1812 msym_symmetry_operation_t c3[2] = {
1813 [0] = {.type = PROPER_ROTATION, .order = 3, .power = 1, .cla = cla+1, .orientation = NONE},
1814 [1] = {.type = PROPER_ROTATION, .order = 3, .power = 2, .cla = cla+1, .orientation = NONE}
1815 };
1816
1817 msym_symmetry_operation_t c4[4] = {
1818 [0] = {.type = PROPER_ROTATION, .order = 5, .power = 1, .cla = cla+2, .orientation = NONE},
1819 [1] = {.type = PROPER_ROTATION, .order = 5, .power = 4, .cla = cla+2, .orientation = NONE},
1820 [2] = {.type = PROPER_ROTATION, .order = 5, .power = 2, .cla = cla+3, .orientation = NONE},
1821 [3] = {.type = PROPER_ROTATION, .order = 5, .power = 3, .cla = cla+3, .orientation = NONE}
1822 };
1823
1824 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsIcosahedral(l,sops, 1, c2, 2, c3, 4, c4, &k))) goto err;
1825
1826 cla += 4;
1827
1828 clean_debug_printf("------ I %d operations %d classes------\n",k-*pk, cla-*pcla);
1829 *pk = k; *pcla = cla;
1830
1831 return ret;
1832 err:
1833 return ret;
1834 }
1835
1836
generateSymmetryOperationsIh(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)1837 msym_error_t generateSymmetryOperationsIh(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
1838 msym_error_t ret = MSYM_SUCCESS;
1839 int k = *pk, cla = *pcla;
1840
1841 msym_symmetry_operation_t c2[2] = {
1842 [0] = {.type = PROPER_ROTATION, .order = 2, .power = 1, .cla = cla, .orientation = NONE},
1843 [1] = {.type = REFLECTION, .order = 1, .power = 1, .cla = cla+1, .orientation = NONE}
1844 };
1845
1846 msym_symmetry_operation_t c3[4] = {
1847 [0] = {.type = PROPER_ROTATION, .order = 3, .power = 1, .cla = cla+2, .orientation = NONE},
1848 [1] = {.type = PROPER_ROTATION, .order = 3, .power = 2, .cla = cla+2, .orientation = NONE},
1849 [2] = {.type = IMPROPER_ROTATION, .order = 6, .power = 1, .cla = cla+3, .orientation = NONE},
1850 [3] = {.type = IMPROPER_ROTATION, .order = 6, .power = 5, .cla = cla+3, .orientation = NONE}
1851 };
1852
1853 msym_symmetry_operation_t c5[8] = {
1854 [0] = {.type = PROPER_ROTATION, .order = 5, .power = 1, .cla = cla+4, .orientation = NONE},
1855 [1] = {.type = PROPER_ROTATION, .order = 5, .power = 4, .cla = cla+4, .orientation = NONE},
1856 [2] = {.type = PROPER_ROTATION, .order = 5, .power = 2, .cla = cla+5, .orientation = NONE},
1857 [3] = {.type = PROPER_ROTATION, .order = 5, .power = 3, .cla = cla+5, .orientation = NONE},
1858 [4] = {.type = IMPROPER_ROTATION, .order = 10, .power = 1, .cla = cla+6, .orientation = NONE},
1859 [5] = {.type = IMPROPER_ROTATION, .order = 10, .power = 9, .cla = cla+6, .orientation = NONE},
1860 [6] = {.type = IMPROPER_ROTATION, .order = 10, .power = 3, .cla = cla+7, .orientation = NONE},
1861 [7] = {.type = IMPROPER_ROTATION, .order = 10, .power = 7, .cla = cla+7, .orientation = NONE},
1862 };
1863
1864 if(MSYM_SUCCESS != (ret = generateSymmetryOperationsIcosahedral(l,sops, 2, c2, 4, c3, 8, c5, &k))) goto err;
1865
1866 if(k - 1 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating Ih operations %d >= %d",k + 12, l); goto err;}
1867
1868 sops[k].type = INVERSION;
1869 sops[k].power = 1;
1870 sops[k].order = 1;
1871 sops[k].orientation = NONE;
1872 sops[k].cla = cla+8;
1873 k++;
1874
1875 cla += 9;
1876
1877 clean_debug_printf("------ Ih %d operations %d classes------\n",k-*pk, cla-*pcla);
1878 *pk = k; *pcla = cla;
1879
1880 return ret;
1881 err:
1882 return ret;
1883 }
1884
generateSymmetryOperationsTetrahedral(int l,msym_symmetry_operation_t sops[l],int c2l,msym_symmetry_operation_t c2[c2l],int csl,msym_symmetry_operation_t cs[csl],int c3l,msym_symmetry_operation_t c3[c3l],int * pk)1885 msym_error_t generateSymmetryOperationsTetrahedral(int l, msym_symmetry_operation_t sops[l], int c2l, msym_symmetry_operation_t c2[c2l], int csl, msym_symmetry_operation_t cs[csl], int c3l, msym_symmetry_operation_t c3[c3l], int *pk){
1886 msym_error_t ret = MSYM_SUCCESS;
1887 int k = *pk;
1888
1889 if(k + c2l*3 + csl*6 + c3l*4 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating tetrahedral operations %d >= %d",k + c2l*3 + csl*6 + c3l*4, l); goto err;}
1890
1891 const double v2[3][3] = {
1892 { 1, 0, 0},
1893 { 0, 1, 0},
1894 { 0, 0, 1},
1895 };
1896
1897 const double vs[6][3] = {
1898 { 1, 0, 1},
1899 { 0, 1, 1},
1900 {-1, 0, 1},
1901 { 0,-1, 1},
1902 { 1, 1, 0},
1903 { 1,-1, 0}
1904 };
1905
1906 const double v3[4][3] = {
1907 {1,1,1},
1908 {-1,1,1},
1909 {1,-1,1},
1910 {-1,-1,1}
1911 };
1912
1913
1914
1915 for(int i = 0;i < c2l;i++){
1916 for(int j = 0; j < 3;j++){
1917 memcpy(&sops[k], &c2[i], sizeof(msym_symmetry_operation_t));
1918 vnorm2(v2[j],sops[k].v);
1919 k++;
1920 }
1921 }
1922
1923 for(int i = 0;i < csl;i++){
1924 for(int j = 0; j < 6;j++){
1925 memcpy(&sops[k], &cs[i], sizeof(msym_symmetry_operation_t));
1926 vnorm2(vs[j],sops[k].v);
1927 k++;
1928 }
1929 }
1930
1931 for(int i = 0;i < c3l;i++){
1932 for(int j = 0; j < 4;j++){
1933 memcpy(&sops[k], &c3[i], sizeof(msym_symmetry_operation_t));
1934 vnorm2(v3[j],sops[k].v);
1935 k++;
1936 }
1937 }
1938
1939 clean_debug_printf("------ Tetra %d operations 0 classes------\n",k-*pk);
1940 *pk = k;
1941
1942 return ret;
1943 err:
1944 return ret;
1945
1946 }
1947
generateSymmetryOperationsOctahedral(int l,msym_symmetry_operation_t sops[l],int c2l,msym_symmetry_operation_t c2[c2l],int c3l,msym_symmetry_operation_t c3[c3l],int c4l,msym_symmetry_operation_t c4[c4l],int * pk)1948 msym_error_t generateSymmetryOperationsOctahedral(int l, msym_symmetry_operation_t sops[l], int c2l, msym_symmetry_operation_t c2[c2l], int c3l, msym_symmetry_operation_t c3[c3l], int c4l, msym_symmetry_operation_t c4[c4l], int *pk){
1949 msym_error_t ret = MSYM_SUCCESS;
1950 int k = *pk;
1951
1952 if(k + c2l*6 + c3l*4 + c4l*3 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating octahedral operations %d >= %d",k + c2l*6 + c3l*4 + c4l*3, l); goto err;}
1953
1954 const double v2[6][3] = {
1955 { 1, 0, 1},
1956 { 0, 1, 1},
1957 {-1, 0, 1},
1958 { 0,-1, 1},
1959 { 1, 1, 0},
1960 { 1,-1, 0}
1961 };
1962
1963 const double v3[4][3] = {
1964 {1,1,1},
1965 {-1,1,1},
1966 {1,-1,1},
1967 {-1,-1,1}
1968 };
1969
1970 const double v4[3][3] = {
1971 { 1, 0, 0},
1972 { 0, 1, 0},
1973 { 0, 0, 1},
1974 };
1975
1976 for(int i = 0;i < c2l;i++){
1977 for(int j = 0; j < 6;j++){
1978 memcpy(&sops[k], &c2[i], sizeof(msym_symmetry_operation_t));
1979 vnorm2(v2[j],sops[k].v);
1980 k++;
1981 }
1982 }
1983
1984 for(int i = 0;i < c3l;i++){
1985 for(int j = 0; j < 4;j++){
1986 memcpy(&sops[k], &c3[i], sizeof(msym_symmetry_operation_t));
1987 vnorm2(v3[j],sops[k].v);
1988 k++;
1989 }
1990 }
1991
1992 for(int i = 0;i < c4l;i++){
1993 for(int j = 0; j < 3;j++){
1994 memcpy(&sops[k], &c4[i], sizeof(msym_symmetry_operation_t));
1995 vnorm2(v4[j],sops[k].v);
1996 k++;
1997 }
1998 }
1999
2000 clean_debug_printf("------ Octa %d operations 0 classes------\n",k-*pk);
2001 *pk = k;
2002
2003 return ret;
2004 err:
2005 return ret;
2006
2007 }
2008
2009
generateSymmetryOperationsIcosahedral(int l,msym_symmetry_operation_t sops[l],int c2l,msym_symmetry_operation_t c2[c2l],int c3l,msym_symmetry_operation_t c3[c3l],int c5l,msym_symmetry_operation_t c5[c5l],int * pk)2010 msym_error_t generateSymmetryOperationsIcosahedral(int l, msym_symmetry_operation_t sops[l], int c2l, msym_symmetry_operation_t c2[c2l], int c3l, msym_symmetry_operation_t c3[c3l], int c5l, msym_symmetry_operation_t c5[c5l], int *pk){
2011 msym_error_t ret = MSYM_SUCCESS;
2012 int k = *pk;
2013
2014 if(k + c2l*15 + c3l*10 + c5l*6 > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating icosahedral operations %d >= %d",k + c2l*15 + c3l*10 + c5l*6, l); goto err;}
2015
2016 const double v2[15][3] = {
2017 {0,0,1},
2018 {1,0,0},
2019 {0,1,0},
2020 {PHI,-(PHI+1),1},
2021 {(PHI+1),1,-PHI},
2022 {1,PHI,(PHI+1)},
2023 {PHI,(PHI+1),1},
2024 {(PHI+1),-1,-PHI},
2025 {-1,PHI,-(PHI+1)},
2026 {(PHI+1),1,PHI},
2027 {1,PHI,-(PHI+1)},
2028 {-PHI,(PHI+1),1},
2029 {1,-PHI,-(PHI+1)},
2030 {PHI,(PHI+1),-1},
2031 {(PHI+1),-1,PHI}
2032 };
2033
2034 const double v3[10][3] = {
2035 {1,1,1},
2036 {-1,1,1},
2037 {1,-1,1},
2038 {-1,-1,1},
2039 {(PHI+1),0,1},
2040 {0,-1,(PHI+1)},
2041 {-1,-(PHI+1),0},
2042 {-1,(PHI+1),0},
2043 {0,1,(PHI+1)},
2044 {(PHI+1),0,-1}
2045 };
2046
2047 const double v5[6][3] = {
2048 {PHI,1,0},
2049 {-PHI,1,0},
2050 {0,PHI,1},
2051 {0,-PHI,1},
2052 {1,0,PHI},
2053 {1,0,-PHI}
2054 };
2055
2056 for(int i = 0;i < c2l;i++){
2057 for(int j = 0; j < 15;j++){
2058 memcpy(&sops[k], &c2[i], sizeof(msym_symmetry_operation_t));
2059 vnorm2(v2[j],sops[k].v);
2060 k++;
2061 }
2062 }
2063
2064 for(int i = 0;i < c3l;i++){
2065 for(int j = 0; j < 10;j++){
2066 memcpy(&sops[k], &c3[i], sizeof(msym_symmetry_operation_t));
2067 vnorm2(v3[j],sops[k].v);
2068 k++;
2069 }
2070 }
2071
2072 for(int i = 0;i < c5l;i++){
2073 for(int j = 0; j < 6;j++){
2074 memcpy(&sops[k], &c5[i], sizeof(msym_symmetry_operation_t));
2075 vnorm2(v5[j],sops[k].v);
2076 k++;
2077 }
2078 }
2079
2080 clean_debug_printf("------ Icosa %d operations 0 classes------\n",k-*pk);
2081 *pk = k;
2082
2083 return ret;
2084 err:
2085 return ret;
2086
2087 }
2088
generateReflectionPlanes(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)2089 msym_error_t generateReflectionPlanes(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
2090 msym_error_t ret = MSYM_SUCCESS;
2091 int k = *pk, cla = *pcla;
2092 double z[3] = {0.0,0.0,1.0}, y[3] = {0.0,1.0,0.0};
2093 msym_symmetry_operation_t sigma = {.type = REFLECTION, .power = 1, .order = 1};
2094 if(k + n > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating reflection planes"); goto err;}
2095 vcopy(y,sigma.v);
2096 enum _msym_symmetry_operation_orientation orientation[2] = {VERTICAL, DIHEDRAL};
2097 for(int i = 0;i < n;i++){
2098 int e = 1 & ~n, ie = ((i & e)), index = k + (i >> e) + (ie ? (n >> 1) : 0); //power = 1 - (ie << 1),
2099 memcpy(&(sops[index]), &sigma, sizeof(msym_symmetry_operation_t));
2100 vrotate(i*M_PI/n, sigma.v, z, sops[index].v);
2101 //sops[index].power = power; // Used the finding of eigenvalues for character tables
2102 sops[index].orientation = orientation[ie];
2103 sops[index].cla = cla + ie;
2104 }
2105
2106 k += n;
2107 cla += 1 << (~n & 1); //1 or 2 added classes
2108 clean_debug_printf("------ R %d operations %d classes------\n",k-*pk, cla-*pcla);
2109 *pk = k; *pcla = cla;
2110
2111 return ret;
2112 err:
2113 return ret;
2114 }
2115
generateC2Axes(int n,int l,msym_symmetry_operation_t sops[l],int * pk,int * pcla)2116 msym_error_t generateC2Axes(int n, int l, msym_symmetry_operation_t sops[l], int *pk, int *pcla){
2117 msym_error_t ret = MSYM_SUCCESS;
2118 int k = *pk, cla = *pcla;
2119 double z[3] = {0.0,0.0,1.0}, x[3] = {1.0,0.0,0.0};
2120 msym_symmetry_operation_t c2 = {.type = PROPER_ROTATION, .order = 2, .power = 1};
2121 if(k + n > l){ret = MSYM_POINT_GROUP_ERROR; msymSetErrorDetails("Too many operations when generating C2 axes"); goto err;}
2122 vcopy(x,c2.v);
2123 enum _msym_symmetry_operation_orientation orientation[2] = {VERTICAL, DIHEDRAL};
2124 for(int i = 0;i < n;i++){
2125 int e = 1 & ~n, ie = ((i & e)), index = k + (i >> e) + (ie ? (n >> 1) : 0);
2126 memcpy(&(sops[index]), &c2, sizeof(msym_symmetry_operation_t));
2127 vrotate(i*M_PI/n, c2.v, z, sops[index].v);
2128 //sops[index].power = power; // Used the finding of eigenvalues for character tables
2129 sops[index].orientation = orientation[ie];
2130 sops[index].cla = cla + ie;
2131 }
2132
2133 k += n;
2134 cla += 1 << (~n & 1); //1 or 2 added classes
2135 clean_debug_printf("------ C2 %d operations %d classes------\n",k-*pk, cla-*pcla);
2136 *pk = k; *pcla = cla;
2137
2138 return ret;
2139 err:
2140 return ret;
2141 }
2142
2143
classifySymmetryOperations(msym_point_group_t * pg)2144 int classifySymmetryOperations(msym_point_group_t *pg){
2145 int c = 1;
2146 double (*mop)[3][3] = malloc(sizeof(double[pg->order][3][3]));
2147 double (*imop)[3][3] = malloc(sizeof(double[pg->order][3][3]));
2148
2149 //There may be a better way to do this
2150 for(int i = 0; i < pg->order;i++){
2151 if(pg->sops[i].type == IDENTITY){
2152 pg->sops[i].cla = 0;
2153 } else {
2154 pg->sops[i].cla = -1;
2155 }
2156 msym_symmetry_operation_t isop;
2157 invertSymmetryOperation(&(pg->sops[i]), &isop);
2158 symmetryOperationMatrix(&(pg->sops[i]), mop[i]);
2159 symmetryOperationMatrix(&isop, imop[i]);
2160 }
2161
2162 for(int i = 0; i < pg->order;i++){
2163 if(pg->sops[i].cla < 0){
2164 pg->sops[i].cla = c;
2165 for(int j = 0; j < pg->order;j++){
2166 double m[3][3];
2167 mmmul(mop[i], imop[j], m);
2168 mmmul(mop[j],m,m);
2169 for(int k = 0; k < pg->order;k++){
2170 if(mequal(mop[k],m,CLASSIFICATION_THRESHOLD)){ //Don't need to be dynamic, this is done on generated point groups (always the same)
2171 pg->sops[k].cla = c;
2172 }
2173 }
2174 }
2175 c++;
2176 }
2177 }
2178
2179 free(mop);
2180 free(imop);
2181
2182 return c;
2183
2184
2185 }
2186
2187 //cant be bothereed writing an efficient sorting alg for this
sortSymmetryOperations(msym_point_group_t * pg,int classes)2188 void sortSymmetryOperations(msym_point_group_t *pg, int classes){
2189 msym_symmetry_operation_t *tmp = malloc(pg->order*sizeof(msym_symmetry_operation_t));
2190 int n = 0;
2191
2192 for(int c = 0; c < classes;c++){
2193 for(int i = 0; i < pg->order;i++){
2194 if(pg->sops[i].cla == c){
2195 copySymmetryOperation(&tmp[n], &pg->sops[i]);
2196 n++;
2197 }
2198 }
2199 }
2200 memcpy(pg->sops, tmp,pg->order*sizeof(msym_symmetry_operation_t));
2201
2202 free(tmp);
2203 }
2204
numberOfSubgroups(msym_point_group_t * pg)2205 int numberOfSubgroups(msym_point_group_t *pg){
2206
2207 int n = pg->n;
2208 int size = 0, ndiv = n >= 2, sdiv = 0, nodd = 0, sodd = 0, neven = 0, seven = 0;
2209
2210 if(isLinearSubgroup(pg)){
2211 switch (pg->type) {
2212 case MSYM_POINT_GROUP_TYPE_Cnv : n = pg->order / 4; break;
2213 case MSYM_POINT_GROUP_TYPE_Dnh : n = pg->order / 2; break;
2214 default: break;
2215 }
2216 }
2217
2218 switch (pg->type) {
2219 case MSYM_POINT_GROUP_TYPE_Kh : size = -1; break;
2220 case MSYM_POINT_GROUP_TYPE_K : size = -1; break;
2221 case MSYM_POINT_GROUP_TYPE_Ih : size = 162; break;
2222 case MSYM_POINT_GROUP_TYPE_I : size = 57; break;
2223 case MSYM_POINT_GROUP_TYPE_Oh : size = 96; break;
2224 case MSYM_POINT_GROUP_TYPE_O : size = 28; break;
2225 case MSYM_POINT_GROUP_TYPE_Th : size = 24; break;
2226 case MSYM_POINT_GROUP_TYPE_Td : size = 28; break;
2227 case MSYM_POINT_GROUP_TYPE_T : size = 9; break;
2228 case MSYM_POINT_GROUP_TYPE_Ci : size = 0; break;
2229 case MSYM_POINT_GROUP_TYPE_Cs : size = 0; break;
2230 default: {
2231 for(int i = 2; i < n;i++){
2232 if(n % i == 0){ ndiv++; sdiv += i; }
2233 }
2234 for(int i = 3; i < n;i += 2){
2235 if(n % i == 0){ nodd++; sodd += i;}
2236 }
2237 for(int i = 4; i <= n;i += 2){
2238 if(n % i == 0){ neven++; seven += (n << 1)/i; }
2239 }
2240 switch (pg->type) {
2241 case MSYM_POINT_GROUP_TYPE_Cnh : {
2242 size = 2*ndiv;
2243 if(n % 2 == 0){
2244 int n2 = n >> 1;
2245 for(int i = 2; i < n2;i++){
2246 if(n2 % i == 0){ size++;}
2247 }
2248 size += 1 + (n2 >= 2);
2249 }
2250 break;
2251 }
2252 case MSYM_POINT_GROUP_TYPE_Dn :
2253 case MSYM_POINT_GROUP_TYPE_Cnv : size = n + ndiv + sdiv; break;
2254 case MSYM_POINT_GROUP_TYPE_Cn : size = ndiv - 1; break;
2255 case MSYM_POINT_GROUP_TYPE_Dnh : {
2256 if(n % 2 == 0) size = 4*n + 2*ndiv + 3*sdiv + 4 + neven + seven;
2257 else size = 3*(n+sdiv+1) + 2*ndiv;
2258 break;
2259 }
2260 case MSYM_POINT_GROUP_TYPE_Dnd : {
2261 if(n % 2 == 0) size = 2*n + 3 + ndiv + 2*sdiv + nodd + sodd;
2262 else size = 3*(n+sdiv+1) + 2*ndiv;
2263 break;
2264 }
2265 case MSYM_POINT_GROUP_TYPE_Sn : size = ndiv - 1; break;
2266 default : break;
2267 }
2268 }
2269 }
2270
2271 return size;
2272 }
2273
2274