1 //
2 //  msym.c
3 //  libmsym
4 //
5 //  Created by Marcus Johansson on 30/01/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 <stdio.h>
12 #include <stdlib.h>
13 #include <string.h>
14 #include "msym.h"
15 #include "context.h"
16 #include "symmetry.h"
17 #include "equivalence_set.h"
18 #include "point_group.h"
19 #include "symmetrize.h"
20 #include "linalg.h"
21 #include "subspace.h"
22 
23 #include "debug.h"
24 
msymFindSymmetry(msym_context ctx)25 msym_error_t msymFindSymmetry(msym_context ctx){
26     msym_error_t ret = MSYM_SUCCESS;
27     int elementsl = 0, esl = 0;
28     msym_element_t *elements = NULL;
29     msym_thresholds_t *t = NULL;
30     msym_equivalence_set_t *es = NULL, *des = NULL;;
31     msym_point_group_t *pg = NULL;
32     int sopsl = 0;
33     msym_symmetry_operation_t *sops = NULL;
34     msym_equivalence_set_t *ses = NULL;
35     int sesl = 0;
36     msym_point_group_t *fpg = NULL;
37 
38     if(MSYM_SUCCESS != (ret = ctxGetElements(ctx, &elementsl, &elements))) goto err;
39 
40     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
41 
42     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &des))){
43         if(MSYM_SUCCESS != (ret = msymFindEquivalenceSets(ctx))) goto err;
44     }
45 
46     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))) goto err;
47     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))){
48         if(MSYM_SUCCESS != (ret = findSymmetryOperations(esl,es,t,&sopsl,&sops))) goto err;
49         if(MSYM_SUCCESS != (ret = findPointGroup(sopsl, sops, t, &fpg))) goto err;
50         pg = fpg;
51         if(MSYM_SUCCESS != (ret = ctxSetPointGroup(ctx, pg))) {
52             free(pg);
53             goto err;
54         }
55     }
56 
57     if(NULL != fpg || isLinearPointGroup(pg)){
58         // Reuild equivalence sets after determining poing group in case they are very similar
59         if(MSYM_SUCCESS != (ret = ctxReduceLinearPointGroup(ctx))) goto err;
60 
61         if(MSYM_SUCCESS != (ret = splitPointGroupEquivalenceSets(pg, esl, es, &sesl, &ses, t))) goto err;
62         if(MSYM_SUCCESS != (ret = ctxSetEquivalenceSets(ctx, sesl, ses))) goto err;
63         ses = NULL; sesl = 0;
64         if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))) goto err;
65     }
66 
67     if(MSYM_SUCCESS != (ret = msymFindEquivalenceSetPermutations(ctx))) goto err;
68 
69     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))) goto err; //This is only for printing, since permutation may regenerate sets
70 
71     free(sops);
72     return ret;
73 
74 err:
75     free(ses);
76     free(sops);
77     if(des == NULL) {
78         ctxDestroyEquivalcenceSets(ctx);
79     }
80     return ret;
81 }
82 
msymSetPointGroupByName(msym_context ctx,const char * name)83 msym_error_t msymSetPointGroupByName(msym_context ctx, const char *name){
84     msym_error_t ret = MSYM_SUCCESS;
85     msym_point_group_t *pg = NULL, *ppg = NULL;
86     msym_thresholds_t *t = NULL;
87 
88 
89     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
90 
91     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &ppg))){
92         double transform[3][3];
93         mleye(3,transform);
94         if(MSYM_SUCCESS != (ret = generatePointGroupFromName(name, transform, t, &pg))) goto err;
95     } else if(MSYM_SUCCESS != (ret = generatePointGroupFromName(name, ppg->transform, t, &pg))) goto err;
96 
97     if(MSYM_SUCCESS != (ret = ctxSetPointGroup(ctx, pg))) goto err;
98 
99     return ret;
100 
101 err:
102     free(pg);
103     return ret;
104 }
105 
msymSetPointGroupByType(msym_context ctx,msym_point_group_type_t type,int n)106 msym_error_t msymSetPointGroupByType(msym_context ctx, msym_point_group_type_t type, int n){
107     msym_error_t ret = MSYM_SUCCESS;
108     msym_point_group_t *pg = NULL, *ppg = NULL;
109     msym_thresholds_t *t = NULL;
110 
111 
112     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
113 
114     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &ppg))){
115         double transform[3][3];
116         mleye(3,transform);
117         if(MSYM_SUCCESS != (ret = generatePointGroupFromType(type, n, transform, t, &pg))) goto err;
118     } else if(MSYM_SUCCESS != (ret = generatePointGroupFromType(type, n, ppg->transform, t, &pg))) goto err;
119 
120     if(MSYM_SUCCESS != (ret = ctxSetPointGroup(ctx, pg))) goto err;
121 
122     return ret;
123 
124 err:
125     free(pg);
126     return ret;
127 }
128 
msymGenerateElements(msym_context ctx,int length,msym_element_t elements[length])129 msym_error_t msymGenerateElements(msym_context ctx, int length, msym_element_t elements[length]){
130     msym_error_t ret = MSYM_SUCCESS;
131     msym_point_group_t *pg = NULL;
132     msym_thresholds_t *t = NULL;
133     msym_element_t *gelements = NULL;
134     msym_equivalence_set_t *es = NULL;
135     msym_element_t **pelements = NULL;
136     double err = 0.0;
137     double cm[3];
138 
139     int glength = 0, plength = 0, esl = 0;
140     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
141     if(MSYM_SUCCESS != (ret = msymGetCenterOfMass(ctx, cm))) goto err;
142 
143     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
144     if(MSYM_SUCCESS != (ret = generateEquivalenceSet(pg, length, elements, cm, &glength, &gelements, &esl, &es,t))) goto err;
145     if(MSYM_SUCCESS != (ret = ctxSetElements(ctx, glength, gelements))) goto err;
146     if(MSYM_SUCCESS != (ret = ctxGetElementPtrs(ctx, &plength, &pelements))) goto err;
147     if(plength != glength){
148         ret = MSYM_INVALID_ELEMENTS;
149         msymSetErrorDetails("Inconsistency detected when setting elements");
150         goto err;
151     }
152     for(int i = 0;i < esl;i++){
153         for(int j = 0;j < es[i].length;j++){
154             long int index = es[i].elements[j] - gelements;
155             es[i].elements[j] = pelements[index];
156         }
157     }
158     if(MSYM_SUCCESS != (ret = ctxSetEquivalenceSets(ctx, esl, es))) goto err;
159     es = NULL; esl = 0;
160     if(MSYM_SUCCESS != (ret = msymFindEquivalenceSetPermutations(ctx))) goto err;
161     if(MSYM_SUCCESS != (ret = msymSymmetrizeElements(ctx, &err))) goto err;
162     if(MSYM_SUCCESS != (ret = msymSetCenterOfMass(ctx, cm))) goto err;
163     free(gelements);
164     return ret;
165 
166 err:
167     free(gelements);
168     free(es);
169     return ret;
170 }
171 
msymFindEquivalenceSets(msym_context ctx)172 msym_error_t msymFindEquivalenceSets(msym_context ctx){
173     msym_error_t ret = MSYM_SUCCESS;
174     int pelementsl = 0;
175     msym_element_t **pelements = NULL;
176     msym_thresholds_t *t = NULL;
177     msym_point_group_t *pg = NULL;
178     msym_geometry_t g = MSYM_GEOMETRY_UNKNOWN;
179     double eigvec[3][3];
180     double eigval[3];
181     int esl = 0;
182     msym_equivalence_set_t *es;
183 
184     if(MSYM_SUCCESS != (ret = ctxGetElementPtrs(ctx, &pelementsl, &pelements))) goto err;
185     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
186     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) {
187         if(MSYM_SUCCESS != (ret = ctxGetGeometry(ctx, &g, eigval, eigvec))) goto err;
188         if(MSYM_SUCCESS != (ret = findEquivalenceSets(pelementsl, pelements, g, &esl, &es, t))) goto err;
189     } else {
190         if(MSYM_SUCCESS != (ret = findPointGroupEquivalenceSets(pg, pelementsl, pelements, &esl, &es, t))) goto err;
191     }
192     if(MSYM_SUCCESS != (ret = ctxSetEquivalenceSets(ctx, esl, es))) goto err;
193 err:
194     return ret;
195 }
196 
197 
msymAlignAxes(msym_context ctx)198 msym_error_t msymAlignAxes(msym_context ctx){
199 
200     msym_error_t ret = MSYM_SUCCESS;
201     msym_element_t *elements = NULL;
202     msym_point_group_t *pg;
203     int elementsl = 0;
204     double zero[3] = {0,0,0};
205 
206     if(MSYM_SUCCESS != (ret = ctxGetElements(ctx, &elementsl, &elements))) goto err;
207     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
208 
209 
210     if(pg->sops == NULL || pg->order == 0){
211         msymSetErrorDetails("No symmetry operations in point group");
212         ret = MSYM_INVALID_POINT_GROUP;
213         goto err;
214     }
215 
216     if(MSYM_SUCCESS != (ret = msymSetCenterOfMass(ctx, zero))) goto err;
217 
218     for(int i = 0; i < elementsl;i++) mvmul(elements[i].v, pg->transform, elements[i].v);
219     for(int i = 0; i < pg->order;i++) mvmul(pg->sops[i].v, pg->transform, pg->sops[i].v);
220     mleye(3,pg->transform);
221     if(MSYM_SUCCESS != (ret = ctxUpdateExternalElementCoordinates(ctx))) goto err;
222 
223 err:
224     return ret;
225 }
226 
msymGetAlignmentAxes(msym_context ctx,double primary[3],double secondary[3])227 msym_error_t msymGetAlignmentAxes(msym_context ctx, double primary[3], double secondary[3]){
228     msym_error_t ret = MSYM_SUCCESS;
229     msym_point_group_t *pg;
230 
231     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
232 
233     double m[3][3], x[3] = {1,0,0}, z[3] = {0,0,1};
234     minv(pg->transform,m);
235     mvmul(z, m, primary);
236     mvmul(x, m, secondary);
237 
238 err:
239     return ret;
240 
241 }
242 
243 
msymGetAlignmentTransform(msym_context ctx,double transform[3][3])244 msym_error_t msymGetAlignmentTransform(msym_context ctx, double transform[3][3]){
245     msym_error_t ret = MSYM_SUCCESS;
246     msym_point_group_t *pg;
247 
248     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
249 
250     mcopy(pg->transform, transform);
251 
252 err:
253     return ret;
254 
255 }
256 
msymSetAlignmentTransform(msym_context ctx,double transform[3][3])257 msym_error_t msymSetAlignmentTransform(msym_context ctx, double transform[3][3]){
258     msym_error_t ret = MSYM_SUCCESS;
259     msym_point_group_t *pg;
260     msym_element_t *elements = NULL;
261     msym_thresholds_t *t = NULL;
262     msym_equivalence_set_t *es = NULL;
263     int elementsl = 0, esl = 0;
264     double m[3][3];
265 
266     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
267     if(MSYM_SUCCESS != (ret = ctxGetElements(ctx, &elementsl, &elements))){
268         elements = NULL;
269         elementsl = 0;
270     }
271 
272     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))){
273         es = NULL;
274         esl = 0;
275     }
276 
277     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
278 
279     if(pg->sops == NULL || pg->order == 0){
280         msymSetErrorDetails("No symmetry operations in point group for setting alignment transform");
281         ret = MSYM_INVALID_POINT_GROUP;
282         goto err;
283     }
284     /* Don't transform elements if we don't have an equivalence set the current pg is set manually */
285     if(NULL != es){
286         for(int i = 0; i < elementsl;i++) mvmul(elements[i].v, pg->transform, elements[i].v);
287     }
288     for(int i = 0; i < pg->order;i++) mvmul(pg->sops[i].v, pg->transform, pg->sops[i].v);
289 
290     minv(transform,m);
291     mcopy(transform, pg->transform);
292     if(NULL != es){
293         for(int i = 0; i < elementsl;i++) mvmul(elements[i].v, m, elements[i].v);
294     }
295     for(int i = 0; i < pg->order;i++) mvmul(pg->sops[i].v, m, pg->sops[i].v);
296 
297 err:
298     return ret;
299 }
300 
msymSetAlignmentAxes(msym_context ctx,double primary[3],double secondary[3])301 msym_error_t msymSetAlignmentAxes(msym_context ctx, double primary[3], double secondary[3]){
302 
303     msym_error_t ret = MSYM_SUCCESS;
304     msym_point_group_t *pg;
305     msym_element_t *elements = NULL;
306     msym_thresholds_t *t = NULL;
307     msym_equivalence_set_t *es = NULL;
308     int elementsl = 0, esl = 0;
309     double x[3] = {1,0,0}, z[3] = {0,0,1}, m[3][3], p[3], s[3];
310 
311     vnorm2(primary, p);
312     vnorm2(secondary,s);
313 
314     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
315     if(MSYM_SUCCESS != (ret = ctxGetElements(ctx, &elementsl, &elements))){
316         elements = NULL;
317         elementsl = 0;
318     }
319 
320     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))){
321         es = NULL;
322         esl = 0;
323     }
324 
325     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
326 
327     if(pg->sops == NULL || pg->order == 0){
328         msymSetErrorDetails("No symmetry operations in point group for setting alignment axes");
329         ret = MSYM_INVALID_POINT_GROUP;
330         goto err;
331     }
332 
333     if(!vperpendicular(primary, secondary, t->angle)) {
334         msymSetErrorDetails("Alignment axes are not orthogonal");
335         ret = MSYM_INVALID_AXES;
336         goto err;
337     }
338 
339     /* Don't transform elements if we don't have an equivalence set the current pg is set manually */
340     if(NULL != es){
341         for(int i = 0; i < elementsl;i++) mvmul(elements[i].v, pg->transform, elements[i].v);
342     }
343     for(int i = 0; i < pg->order;i++) mvmul(pg->sops[i].v, pg->transform, pg->sops[i].v);
344 
345     vproj_plane(s, p, s);
346     malign(p,z,pg->transform);
347     mvmul(s, pg->transform, s);
348     malign(s,x,m);
349     mmmul(m,pg->transform,pg->transform);
350     minv(pg->transform,m);
351 
352     if(NULL != es){
353         for(int i = 0; i < elementsl;i++) mvmul(elements[i].v, m, elements[i].v);
354     }
355     for(int i = 0; i < pg->order;i++) mvmul(pg->sops[i].v, m, pg->sops[i].v);
356 
357 
358 err:
359     return ret;
360 }
361 
msymSelectSubgroup(msym_context ctx,const msym_subgroup_t * sg)362 msym_error_t msymSelectSubgroup(msym_context ctx, const msym_subgroup_t *sg){
363     msym_error_t ret = MSYM_SUCCESS;
364     msym_subgroup_t *sgs;
365     msym_point_group_t *pg;
366     msym_thresholds_t *t = NULL;
367     int sgl = 0;
368 
369     if(MSYM_SUCCESS != (ret = ctxGetSubgroups(ctx, &sgl, &sgs))) goto err;
370     if(sg < sgs || sg >= sgs + sgl){
371         msymSetErrorDetails("Subgroup not available in current context");
372         ret = MSYM_INVALID_SUBGROUPS;
373         goto err;
374     }
375     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
376     if(MSYM_SUCCESS != (ret = pointGroupFromSubgroup(sg, t, &pg))) goto err;
377     if(MSYM_SUCCESS != (ret = ctxSetPointGroup(ctx, pg))) goto err;
378     if(MSYM_SUCCESS != (ret = msymFindEquivalenceSets(ctx))) goto err;
379     if(MSYM_SUCCESS != (ret = msymFindEquivalenceSetPermutations(ctx))) goto err;
380 
381 err:
382     return ret;
383 }
384 
msymSymmetrizeElements(msym_context ctx,double * oerr)385 msym_error_t msymSymmetrizeElements(msym_context ctx, double *oerr){
386     msym_error_t ret = MSYM_SUCCESS;
387 
388     msym_point_group_t *pg = NULL;
389     msym_equivalence_set_t *es = NULL;
390     msym_element_t *elements = NULL;
391 
392     msym_permutation_t **perm = NULL;
393     msym_thresholds_t *t = NULL;
394     double error = 0.0;
395     int perml = 0, esl = 0, elementsl = 0, sopsl = 0;
396 
397     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
398     if(MSYM_SUCCESS != (ret = ctxGetElements(ctx, &elementsl, &elements))) goto err;
399     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
400     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))){
401         if(MSYM_SUCCESS != (ret = msymFindEquivalenceSets(ctx))) goto err;
402         if(MSYM_SUCCESS != (ret = msymFindEquivalenceSetPermutations(ctx))) goto err;
403         if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))) goto err;
404     }
405     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSetPermutations(ctx, &perml, &sopsl, &perm))) goto err;
406     if(sopsl != pg->order || perml != esl) {
407         msymSetErrorDetails("Detected inconsistency between point group, equivalence sets and permutaions");
408         ret = MSYM_INVALID_PERMUTATION;
409         goto err;
410     }
411 
412     if(MSYM_SUCCESS != (ret = symmetrizeElements(pg, esl, es, perm, t, &error))) goto err;
413 
414     if(MSYM_SUCCESS != (ret = ctxUpdateGeometry(ctx))) goto err;
415 
416     if(MSYM_SUCCESS != (ret = ctxUpdateExternalElementCoordinates(ctx))) goto err;
417 
418     *oerr = error;
419 err:
420     return ret;
421 }
422 
msymApplyTranslation(msym_context ctx,msym_element_t * ext,double v[3])423 msym_error_t msymApplyTranslation(msym_context ctx, msym_element_t *ext, double v[3]){
424     msym_error_t ret = MSYM_SUCCESS;
425 
426     msym_point_group_t *pg;
427     msym_equivalence_set_t *es, *ees;
428     msym_element_t *eelements;
429     msym_equivalence_set_t **eesmap = NULL;
430     msym_permutation_t **perm;
431     msym_thresholds_t *t = NULL;
432     int perml = 0, esl = 0, eesl = 0, eelementsl = 0, sopsl = 0;
433 
434     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
435     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
436     if(MSYM_SUCCESS != (ret = ctxGetExternalElements(ctx, &eelementsl, &eelements))) goto err;
437     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))){
438         if(MSYM_SUCCESS != (ret = msymFindEquivalenceSets(ctx))) goto err;
439         if(MSYM_SUCCESS != (ret = msymFindEquivalenceSetPermutations(ctx))) goto err;
440         if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))) goto err;
441     }
442     if(MSYM_SUCCESS != (ret = ctxGetExternalEquivalenceSets(ctx, &eesl, &ees))) goto err;
443     if(MSYM_SUCCESS != (ret = ctxGetExternalElementEquivalenceSetMap(ctx, &eesmap))) goto err;
444 
445     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSetPermutations(ctx, &perml, &sopsl, &perm))) goto err;
446     if(sopsl != pg->order || perml != esl) {
447         msymSetErrorDetails("Detected inconsistency between point group, equivalence sets and permutaions");
448         ret = MSYM_INVALID_PERMUTATION;
449         goto err;
450     }
451 
452     int esmi = (int)(ext - eelements);
453 
454     if(esmi > eelementsl) {
455         msymSetErrorDetails("Element outside of memory block of external elements");
456         ret = MSYM_INVALID_ELEMENTS;
457         goto err;
458     }
459 
460     int fesi = (int)(eesmap[esmi] - ees);
461     msym_equivalence_set_t *fes = eesmap[esmi];
462     int fi = 0;
463     for(fi = 0;fi < fes->length;fi++){
464         if(fes->elements[fi] == ext) break;
465     }
466 
467     if(fi >= fes->length){
468         msymSetErrorDetails("Could not find index of element %s in equivalence set %d", ext->name, fesi);
469         ret = MSYM_INVALID_ELEMENTS;
470         goto err;
471     }
472 
473     if(MSYM_SUCCESS != (ret = symmetrizeTranslation(pg, &es[fesi], perm[fesi], fi, v))) goto err;
474 
475     if(MSYM_SUCCESS != (ret = ctxUpdateExternalElementCoordinates(ctx))) goto err;
476 
477     return ret;
478 err:
479     return ret;
480 }
481 
msymGenerateSubrepresentationSpaces(msym_context ctx)482 msym_error_t msymGenerateSubrepresentationSpaces(msym_context ctx){
483     msym_error_t ret = MSYM_SUCCESS;
484 
485     msym_point_group_t *pg = NULL;
486     msym_basis_function_t *basis = NULL;
487     msym_equivalence_set_t *es = NULL;
488     msym_equivalence_set_t **eesmap = NULL;
489     msym_permutation_t **perm = NULL;
490     msym_thresholds_t *t = NULL;
491     msym_subrepresentation_space_t *srs = NULL;
492     msym_basis_function_t **srsbf = NULL;
493     msym_element_t *elements = NULL;
494     const msym_subgroup_t *sg = NULL;
495     int *span = NULL;
496 
497     int basisl = 0, esl = 0, perml = 0, sopsl = 0, srsl = 0, elementsl = 0, sgl = 0;
498 
499     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
500     if(MSYM_SUCCESS != (ret = ctxGetExternalElements(ctx, &elementsl, &elements))) goto err;
501     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
502     if(pg->ct == NULL){
503         if(MSYM_SUCCESS != (ret = generateCharacterTable(pg->type, pg->n, pg->order, pg->sops, &pg->ct))) goto err;
504     }
505     if(MSYM_SUCCESS != (ret = ctxGetExternalEquivalenceSets(ctx, &esl, &es))) goto err;
506     if(MSYM_SUCCESS != (ret = ctxGetExternalElementEquivalenceSetMap(ctx, &eesmap))) goto err;
507     if(MSYM_SUCCESS != (ret = ctxGetBasisFunctions(ctx, &basisl, &basis))) goto err;
508     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSetPermutations(ctx, &perml, &sopsl, &perm))) goto err;
509     if(sopsl != pg->order || perml != esl) {ret = MSYM_INVALID_PERMUTATION; goto err;}
510 
511     if(MSYM_SUCCESS != (ret = msymGetSubgroups(ctx, &sgl, &sg))) goto err;
512 
513     if(MSYM_SUCCESS != (ret = generateSubrepresentationSpaces(pg, sgl, sg, esl, es, perm, basisl, basis, elements, eesmap, t, &srsl, &srs, &srsbf, &span))) goto err;
514 
515     if(MSYM_SUCCESS != (ret = ctxSetSubrepresentationSpaces(ctx,srsl,srs,srsbf,span))) goto err;
516 
517     return ret;
518 err:
519     freeSubrepresentationSpaces(srsl, srs);
520     free(srs);
521     free(span);
522     return ret;
523 }
524 
msymGetSALCs(msym_context ctx,int l,double c[l][l],int species[l],msym_partner_function_t pf[l])525 msym_error_t msymGetSALCs(msym_context ctx, int l, double c[l][l], int species[l], msym_partner_function_t pf[l]){
526     msym_error_t ret = MSYM_SUCCESS;
527 
528 
529     msym_subrepresentation_space_t *srs = NULL;
530     msym_basis_function_t *basis = NULL;
531 
532     int *span = NULL;
533 
534     int srsl = 0, basisl = 0;
535 
536     if(MSYM_SUCCESS != (ret = ctxGetBasisFunctions(ctx, &basisl, &basis))) goto err;
537 
538     if(MSYM_SUCCESS != (ret = ctxGetSubrepresentationSpaces(ctx, &srsl, &srs, &span))){
539         if(MSYM_SUCCESS != (ret = msymGenerateSubrepresentationSpaces(ctx))) goto err;
540         if(MSYM_SUCCESS != (ret = ctxGetSubrepresentationSpaces(ctx, &srsl, &srs, &span))) goto err;
541     }
542 
543     if(l != basisl){
544         ret = MSYM_INVALID_INPUT;
545         msymSetErrorDetails("Supplied SALC matrix size (%dx%d) does not match number of basis functions (%d)",l,l,basisl);
546         goto err;
547     }
548 
549     memset(c,0,sizeof(double[l][l]));
550     int wf = 0;
551     for(int i = 0;i < srsl;i++){
552         int s = srs[i].s;
553         for(int j = 0;j < srs[i].salcl;j++){
554             int pwf = wf;
555             double (*mpf)[srs[i].salc[j].fl] = srs[i].salc[j].pf;
556             for(int d = 0;d < srs[i].salc[j].d;d++){
557                 if(wf >= basisl){
558                     ret = MSYM_INVALID_SUBSPACE;
559                     msymSetErrorDetails("Generated more SALCs than the number of basis functions (%d)", basisl);
560                     goto err;
561                 }
562                 for(int f = 0;f < srs[i].salc[j].fl;f++){
563                     int index = (int)(srs[i].salc[j].f[f] - basis);
564                     c[wf][index] = mpf[d][f];
565                 }
566                 if(NULL != pf){
567                     pf[wf].i = pwf;
568                     pf[wf].d = d;
569                 }
570                 if(NULL != species) species[wf] = s;
571                 wf++;
572             }
573         }
574     }
575 
576     if(wf != basisl){
577         ret = MSYM_INVALID_BASIS_FUNCTIONS;
578         msymSetErrorDetails("Number of generated SALC wavefunctions (%d) does not match orbital basis (%d)",wf,basisl);
579         goto err;
580     }
581 
582 err:
583     return ret;
584 
585 }
586 
msymSymmetrySpeciesComponents(msym_context ctx,int wfl,double * wf,int sl,double * s)587 msym_error_t msymSymmetrySpeciesComponents(msym_context ctx, int wfl, double *wf, int sl, double *s){
588     msym_error_t ret = MSYM_SUCCESS;
589 
590     msym_point_group_t *pg = NULL;
591     msym_subrepresentation_space_t *srs = NULL;
592     msym_basis_function_t *basis = NULL;
593     int *span = NULL;
594 
595     int srsl = 0, basisl = 0;
596 
597     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
598     if(pg->ct == NULL){
599         if(MSYM_SUCCESS != (ret = generateCharacterTable(pg->type, pg->n, pg->order, pg->sops, &pg->ct))) goto err;
600     }
601 
602     if(MSYM_SUCCESS != (ret = ctxGetBasisFunctions(ctx, &basisl, &basis))) goto err;
603 
604     if(basisl != wfl) {
605         ret = MSYM_INVALID_INPUT;
606         msymSetErrorDetails("Supplied coefficient vector size (%d) does not match number of basis functions (%d)",wfl,basisl);
607         goto err;
608     }
609 
610     if(sl != pg->ct->d) {
611         ret = MSYM_INVALID_INPUT;
612         msymSetErrorDetails("Supplied symmetry species vector size (%d) does not match character table (%d)",sl,pg->ct->d);
613         goto err;
614     }
615 
616     if(MSYM_SUCCESS != (ret = ctxGetSubrepresentationSpaces(ctx, &srsl, &srs, &span))){
617         if(MSYM_SUCCESS != (ret = msymGenerateSubrepresentationSpaces(ctx))) goto err;
618         if(MSYM_SUCCESS != (ret = ctxGetSubrepresentationSpaces(ctx, &srsl, &srs, &span))) goto err;
619     }
620 
621     if(MSYM_SUCCESS != (ret = symmetrySpeciesComponents(pg, srsl, srs, basisl, basis, wf, s))) goto err;
622 
623 err:
624     return ret;
625 
626 }
627 
msymSymmetrizeWavefunctions(msym_context ctx,int l,double c[l][l],int species[l],msym_partner_function_t pf[l])628 msym_error_t msymSymmetrizeWavefunctions(msym_context ctx, int l, double c[l][l], int species[l], msym_partner_function_t pf[l]){
629     msym_error_t ret = MSYM_SUCCESS;
630     msym_point_group_t *pg = NULL;
631     msym_subrepresentation_space_t *srs = NULL;
632     msym_basis_function_t *basis = NULL;
633     int *span = NULL;
634 
635     int srsl = 0, basisl = 0;
636 
637     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
638     if(pg->ct == NULL){
639         if(MSYM_SUCCESS != (ret = generateCharacterTable(pg->type, pg->n, pg->order, pg->sops, &pg->ct))) goto err;
640     }
641 
642 
643     if(MSYM_SUCCESS != (ret = ctxGetBasisFunctions(ctx, &basisl, &basis))) goto err;
644 
645     if(basisl != l) {
646         ret = MSYM_INVALID_INPUT;
647         msymSetErrorDetails("Supplied wavefunction matrix size (%d) does not match number of basis functions (%d)",l,basisl);
648         goto err;
649     }
650 
651     if(MSYM_SUCCESS != (ret = ctxGetSubrepresentationSpaces(ctx, &srsl, &srs, &span))){
652         if(MSYM_SUCCESS != (ret = msymGenerateSubrepresentationSpaces(ctx))) goto err;
653         if(MSYM_SUCCESS != (ret = ctxGetSubrepresentationSpaces(ctx, &srsl, &srs, &span))) goto err;
654     }
655 
656     if(MSYM_SUCCESS != (ret = symmetrizeWavefunctions(pg, srsl, srs, span, basisl, basis, c , c, species, pf))) goto err;
657 
658 err:
659     return ret;
660 }
661 
msymFindEquivalenceSetPermutations(msym_context ctx)662 msym_error_t msymFindEquivalenceSetPermutations(msym_context ctx) {
663     msym_error_t ret = MSYM_SUCCESS;
664     //We can't allocate this as a double[][] unless we typecast it every time, since the compiler doesn't have the indexing information in the context
665     msym_permutation_t **perm = NULL;
666     msym_permutation_t *bperm = NULL;
667     msym_point_group_t *pg = NULL;
668     msym_equivalence_set_t *es = NULL;
669     msym_thresholds_t *t = NULL;
670     double (**esv)[3] = NULL;
671     int esl = 0;
672 
673     if(MSYM_SUCCESS != (ret = ctxGetThresholds(ctx, &t))) goto err;
674     if(MSYM_SUCCESS != (ret = ctxGetPointGroup(ctx, &pg))) goto err;
675     if(MSYM_SUCCESS != (ret = ctxGetEquivalenceSets(ctx, &esl, &es))) goto err;
676 
677     perm = (msym_permutation_t**)malloc(esl*sizeof(msym_permutation_t*) + esl*pg->order*sizeof(msym_permutation_t));
678     bperm = (msym_permutation_t*)(perm + esl);
679     memset(bperm,0,esl*pg->order*sizeof(msym_permutation_t));
680 
681 
682 
683     for(int i = 0; i < esl;i++){
684         perm[i] = bperm + i*pg->order;
685         if(es[i].length > pg->order){
686             msymSetErrorDetails("Equivalence set has more elements (%d) than the order of the point group %s (%d)",es[i].length,pg->name,pg->order);
687             ret = MSYM_INVALID_EQUIVALENCE_SET;
688             goto err;
689         }
690     }
691     /*
692     if(perm == NULL){
693         perm = (msym_permutation_t**)malloc(esl*sizeof(msym_permutation_t*) + esl*pg->sopsl*sizeof(msym_permutation_t));
694         bperm = (msym_permutation_t*)(perm + esl);
695         memset(bperm,0,esl*pg->sopsl*sizeof(msym_permutation_t));
696         for(int i = 0; i < esl;i++){ //This really shouldn't happen
697             perm[i] = bperm + i*pg->sopsl;
698             if(es[i].length > pg->order){
699                 msymSetErrorDetails("Equivalence set has more elements (%d) than the order of the point group %s (%d)",es[i].length,pg->name,pg->order);
700                 ret = MSYM_INVALID_EQUIVALENCE_SET;
701                 goto err;
702             }
703         }
704     }*/
705 
706     esv = malloc(sizeof(double (*[pg->order])[3]));
707     for(int i = 0; i < esl;i++){
708         for(int j = 0; j < es[i].length;j++){
709             esv[j] = &es[i].elements[j]->v;
710         }
711 
712         for(int j = 0; j < pg->order;j++){
713             if(MSYM_SUCCESS != (ret = findPermutation(&pg->sops[j], es[i].length, esv, t, &perm[i][j]))) goto err;
714         }
715     }
716 
717     if(MSYM_SUCCESS != (ret = ctxSetEquivalenceSetPermutations(ctx, esl, pg->order, perm))) goto err;
718 
719     free(esv);
720     return ret;
721 
722 err:
723     free(esv);
724     free(perm);
725     return ret;
726 }
727 
728