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