1 //
2 // permutation.c
3 // Symmetry
4 //
5 // Created by Marcus Johansson on 01/02/15.
6 // Copyright (c) 2015 Marcus Johansson.
7 //
8 // Distributed under the MIT License ( See LICENSE file or copy at http://opensource.org/licenses/MIT )
9 //
10
11 #include <stdlib.h>
12 #include <string.h>
13
14 #include "msym.h"
15 #include "permutation.h"
16 #include "linalg.h"
17
18 #include "debug.h"
19
20 msym_error_t setPermutationCycles(msym_permutation_t *perm);
21
freePermutationData(msym_permutation_t * perm)22 void freePermutationData(msym_permutation_t *perm){
23 if(perm != NULL){
24 free(perm->c);
25 free(perm->p);
26 }
27 }
28
findPermutation(msym_symmetry_operation_t * sop,int l,double (* v[l])[3],msym_thresholds_t * t,msym_permutation_t * perm)29 msym_error_t findPermutation(msym_symmetry_operation_t *sop, int l, double (*v[l])[3], msym_thresholds_t *t, msym_permutation_t *perm){
30 msym_error_t ret = MSYM_SUCCESS;
31 double m[3][3];
32 symmetryOperationMatrix(sop, m);
33
34 perm->p = malloc(sizeof(int[l]));
35 memset(perm->p, -1, sizeof(int[l])); //TODO: 2s complement
36 perm->p_length = l;
37
38 for(int i = 0; i < l;i++){
39 int j;
40 double r[3];
41 mvmul(*v[i], m, r);
42 for(j = 0;j < l;j++){
43 if(vequal(r, *v[j],t->permutation)){
44 perm->p[i] = j;
45 break;
46 }
47 }
48 if(j == l) {
49 char buf[16];
50 symmetryOperationName(sop, sizeof(buf), buf);
51 msymSetErrorDetails("Unable to determine permutation for symmetry operation %s",buf);
52 ret = MSYM_PERMUTATION_ERROR;
53 goto err;
54 }
55 }
56 if(MSYM_SUCCESS != (ret = setPermutationCycles(perm))) goto err;
57
58 return ret;
59
60 err:
61 free(perm->p);
62 return ret;
63 }
64
65
66 typedef struct _perm_subgroup {
67 int sopsl;
68 int *sops;
69 int subgroup[2];
70 } perm_subgroup_t;
71
72
73
findPermutationSubgroups(int l,msym_permutation_t perm[l],int sgmax,msym_symmetry_operation_t * sops,int * subgroupl,msym_subgroup_t ** subgroup)74 msym_error_t findPermutationSubgroups(int l, msym_permutation_t perm[l], int sgmax, msym_symmetry_operation_t *sops, int *subgroupl, msym_subgroup_t **subgroup){
75 msym_error_t ret = MSYM_SUCCESS;
76 perm_subgroup_t *group = calloc(l, sizeof(perm_subgroup_t));
77
78 int *isops = malloc(sizeof(int[l]));
79 int *msops = malloc(sizeof(int[l]));
80 int gl = 0, glm = 0;
81
82 for(int i = 0;i < l;i++){
83 if((sops[i].power == 1 && (sops[i].type == PROPER_ROTATION || sops[i].type == IMPROPER_ROTATION)) || sops[i].type == INVERSION || sops[i].type == REFLECTION){
84 msym_permutation_cycle_t* c = perm[i].c;
85 glm = gl;
86 memset(msops, 0, sizeof(int[l]));
87 group[gl].sopsl = c->l;
88 group[gl].sops = realloc(group[gl].sops, sizeof(int[c->l]));
89 memset(group[gl].sops,0,sizeof(int[c->l]));
90 //group[gl].sops = calloc(c->l, sizeof(int));
91 group[gl].subgroup[0] = group[gl].subgroup[1] = -1;
92 for(int next = c->s, j = 0;j < c->l;j++){
93 msops[next] = 1;
94 group[gl].sops[j] = next;
95 next = perm[i].p[next];
96 }
97
98 int n = 0;
99 for(int k = 0;k < l && n < l;k++){
100 if(msops[k]){
101 group[gl].sops[n] = k;
102 n++;
103 }
104 }
105 gl += n < l;
106 }
107 }
108
109 if(glm == gl){
110 free(group[gl].sops);
111 }
112
113 for(int i = 0;i < gl && gl < sgmax;i++){
114 for(int j = i+1;j < gl && gl < sgmax;j++){
115 int minl = group[i].sopsl < group[j].sopsl ? group[i].sopsl : group[j].sopsl;
116 if(0 == memcmp(group[i].sops,group[j].sops,sizeof(int[minl]))) continue;
117
118 int n = 0;
119 memset(isops, 0, sizeof(int[l]));
120 memset(msops, 0, sizeof(int[l]));
121
122 for(int k = 0;k < group[i].sopsl;k++){
123 int s = group[i].sops[k];
124 msops[s] = 1;
125 isops[k] = s;
126 }
127 n = group[i].sopsl;
128 for(int k = 0;k < group[j].sopsl;k++){
129 int s = group[j].sops[k];
130 if(msops[s] == 0){
131 msops[s] = 1;
132 isops[n] = s;
133 n++;
134 }
135 }
136 for(int p = 0;p < n && n < l;p++){
137 for(int q = 0;q < n && n < l;q++){
138 int next = perm[isops[p]].p[isops[q]];
139 if(msops[next] == 0){
140 msops[next] = 1;
141 isops[n] = next;
142 n++;
143 }
144 }
145 }
146
147 if(n < l && n > 1) {
148 n = 0;
149 memset(isops, 0, sizeof(int[l]));
150
151 for(int k = 0;k < l;k++){
152 if(msops[k]){
153 isops[n] = k;
154 n++;
155 }
156 }
157 int f;
158 for(f = 0;f < gl;f++){
159 if(group[f].sopsl == n && 0 == memcmp(group[f].sops, isops, sizeof(int[n]))){
160 break;
161 }
162 }
163 if(f == gl){
164 group = realloc(group, sizeof(perm_subgroup_t[gl+1]));
165 group[gl].sopsl = n;
166 group[gl].sops = malloc(sizeof(int[n]));
167 memcpy(group[gl].sops, isops, sizeof(int[n]));
168 group[gl].subgroup[0] = i;
169 group[gl].subgroup[1] = j;
170 gl++;
171 }
172 }
173 }
174 }
175
176 msym_subgroup_t *mgroup = calloc(gl, sizeof(msym_subgroup_t));
177 for(int i = 0;i < gl;i++){
178 mgroup[i].sops = calloc(group[i].sopsl, sizeof(msym_symmetry_operation_t *));
179 mgroup[i].order = group[i].sopsl;
180 mgroup[i].generators[0] = group[i].subgroup[0] < 0 ? NULL : &mgroup[group[i].subgroup[0]];
181 mgroup[i].generators[1] = group[i].subgroup[1] < 0 ? NULL : &mgroup[group[i].subgroup[1]];
182
183 for(int j = 0;j < group[i].sopsl;j++){
184 mgroup[i].sops[j] = &sops[group[i].sops[j]];
185 }
186 }
187
188 *subgroup = mgroup;
189 *subgroupl = gl;
190
191 //err:
192 for(int i = 0;i < gl;i++){
193 free(group[i].sops);
194 }
195 free(group);
196 free(isops);
197 free(msops);
198 return ret;
199 }
200
findSymmetryOperationPermutations(int l,msym_symmetry_operation_t sops[l],msym_thresholds_t * t,msym_permutation_t ** rperm)201 msym_error_t findSymmetryOperationPermutations(int l, msym_symmetry_operation_t sops[l], msym_thresholds_t *t, msym_permutation_t **rperm){
202
203 msym_error_t ret = MSYM_SUCCESS;
204 //Don't block allocate this, it's a pain to keep track of the pointers
205 msym_permutation_t *permutations = malloc(sizeof(msym_permutation_t[l]));
206
207 for(int i = 0; i < l;i++){
208 permutations[i].p = malloc(sizeof(int[l]));
209 memset(permutations[i].p, -1, sizeof(int[l]));
210 permutations[i].p_length = l;
211 }
212
213 double (*msops)[3][3] = malloc(sizeof(double[l][3][3]));
214
215 for(int i = 0; i < l;i++){
216 symmetryOperationMatrix(&sops[i], msops[i]);
217 }
218
219 for(int i = 0; i < l;i++){
220 if((sops[i].type == PROPER_ROTATION && sops[i].order == 0) || sops[i].type == IDENTITY){
221 for(int j = 0;j < l;j++) permutations[i].p[j] = j;
222 } else {
223 for(int j = 0; j < l;j++){
224 int k;
225 double rsop[3][3];
226 mmmul(msops[i], msops[j], rsop);
227 for(k = 0;k < l;k++){
228 if(mequal(rsop,msops[k],t->permutation)){
229 permutations[i].p[j] = k;
230 break;
231 }
232 }
233 if(k == l){
234 char buf1[16];
235 char buf2[16];
236 symmetryOperationName(&sops[i], sizeof(buf1), buf1);
237 symmetryOperationName(&sops[j], sizeof(buf2), buf2);
238 msymSetErrorDetails("Unable to determine permutation for symmetry operation %s and %s",buf1, buf2);
239 ret = MSYM_PERMUTATION_ERROR;
240 goto err;
241 }
242 }
243 }
244 }
245
246 for(int i = 0; i < l;i++){
247 if(MSYM_SUCCESS != (ret = setPermutationCycles(&permutations[i]))) goto err;
248 }
249
250 free(msops);
251 *rperm = permutations;
252
253 return ret;
254
255 err:
256 free(msops);
257 for(int i = 0; i < l;i++){
258 free(permutations[i].p);
259 }
260 free(permutations);
261 *rperm = NULL;
262 return ret;
263
264 }
265
setPermutationCycles(msym_permutation_t * perm)266 msym_error_t setPermutationCycles(msym_permutation_t *perm){
267 msym_error_t ret = MSYM_SUCCESS;
268 int l = perm->p_length;
269 int *icycle = malloc(sizeof(int[l]));
270 int *pcycle = malloc(sizeof(int[l]));
271 int *lcycle = malloc(sizeof(int[l]));
272
273 int cl = 0;
274 memset(icycle, -1,sizeof(int[l])); //TODO: 2s complement
275 memset(lcycle, 0,sizeof(int[l]));
276
277 perm->c = NULL;
278 perm->c_length = 0;
279
280 for(int i = 0; i < l;i++){
281 if(icycle[i] >= 0) continue;
282 lcycle[cl] = 1;
283 pcycle[cl] = i;
284 icycle[i] = cl;
285 for(int next = perm->p[i], loop = 0; next != i;next = perm->p[next]){
286 if(loop++ > l) {
287 msymSetErrorDetails("Encountered loop when determining permutation cycle");
288 ret = MSYM_PERMUTATION_ERROR;
289 goto err;
290 }
291 icycle[next] = cl;
292 lcycle[cl]++;
293 };
294 cl++;
295 }
296 perm->c_length = cl;
297 perm->c = malloc(sizeof(msym_permutation_cycle_t[cl]));
298 for(int c = 0; c < cl;c++){
299 perm->c[c].l = lcycle[c];
300 perm->c[c].s = pcycle[c];
301 }
302
303
304 err:
305 free(icycle);
306 free(pcycle);
307 free(lcycle);
308 return ret;
309 }
310
311
312 //We need these as doubles later, so might as well, even though we could represent these with ALOT less memory
permutationMatrix(msym_permutation_t * perm,double m[perm->p_length][perm->p_length])313 void permutationMatrix(msym_permutation_t *perm, double m[perm->p_length][perm->p_length]){
314 memset(m, 0, sizeof(double[perm->p_length][perm->p_length]));
315 for(int i = 0;i < perm->p_length;i++){
316 //m[i][perm->p[i]] = 1.0;
317 m[perm->p[i]][i] = 1.0;
318 }
319 }
320
321