1 //
2 // equivalence_set.c
3 // libmsym
4 //
5 // Created by Marcus Johansson on 08/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 #include <math.h>
14 #include <time.h>
15 #include "equivalence_set.h"
16 #include "linalg.h"
17 #include "context.h"
18 #include "elements.h"
19
20 #include "debug.h"
21
22 #define SQR(x) ((x)*(x))
23
24 msym_error_t partitionEquivalenceSets(int length, msym_element_t *elements[length], msym_element_t *pelements[length], msym_geometry_t g, int *esl, msym_equivalence_set_t **es, msym_thresholds_t *thresholds);
25 msym_error_t partitionPointGroupEquivalenceSets(msym_point_group_t *pg, int length, msym_element_t *elements[length], msym_element_t *pelements[length], int *esl, msym_equivalence_set_t **es, msym_thresholds_t *thresholds);
26
27
28
copyEquivalenceSets(int length,msym_equivalence_set_t es[length],msym_equivalence_set_t ** ces)29 msym_error_t copyEquivalenceSets(int length, msym_equivalence_set_t es[length], msym_equivalence_set_t **ces){
30 msym_error_t ret = MSYM_SUCCESS;
31 int el = 0;
32
33 for(int i = 0;i < length;i++) el += es[i].length;
34 msym_equivalence_set_t *nes = malloc(sizeof(msym_equivalence_set_t[length]) + sizeof(msym_element_t *[el]));
35 msym_element_t **ep = (msym_element_t **) &es[length];
36 msym_element_t **nep = (msym_element_t **) &nes[length];
37 memcpy(nes, es, sizeof(msym_equivalence_set_t[length]) + sizeof(msym_element_t *[el]));
38 for(int i = 0;i < length;i++) nes[i].elements = nes[i].elements - ep + nep;
39 *ces = nes;
40 //err:
41 return ret;
42 }
43
44 //TODO: Use a preallocated pointer array instead of multiple mallocs
generateEquivalenceSet(msym_point_group_t * pg,int length,msym_element_t elements[length],double cm[3],int * glength,msym_element_t ** gelements,int * esl,msym_equivalence_set_t ** es,msym_thresholds_t * thresholds)45 msym_error_t generateEquivalenceSet(msym_point_group_t *pg, int length, msym_element_t elements[length], double cm[3], int *glength, msym_element_t **gelements, int *esl, msym_equivalence_set_t **es,msym_thresholds_t *thresholds){
46 msym_error_t ret = MSYM_SUCCESS;
47 msym_element_t *ge = calloc(length,sizeof(msym_element_t[pg->order]));
48 msym_equivalence_set_t *ges = calloc(length,sizeof(msym_equivalence_set_t));
49 int gel = 0;
50 int gesl = 0;
51
52 if(pg->order <= 0){
53 ret = MSYM_INVALID_POINT_GROUP;
54 msymSetErrorDetails("Point group of zero order when determining equivalence set");
55 goto err;
56 }
57
58 for(int i = 0;i < length;i++){
59 double ev[3];
60 vsub(elements[i].v, cm, ev);
61 msym_equivalence_set_t *aes = NULL;
62 int f;
63 for(f = 0;f < gel;f++){
64 if(ge[f].n == elements[i].n && ge[f].m == elements[i].m && 0 == strncmp(ge[f].name, elements[i].name, sizeof(ge[f].name)) && vequal(ge[f].v, ev, thresholds->permutation)){
65 break;
66 }
67 }
68 if(f == gel){
69 aes = &ges[gesl++];
70 aes->elements = calloc(pg->order,sizeof(msym_element_t*));
71 aes->length = 0;
72 } else {
73 continue;
74 }
75
76 for(msym_symmetry_operation_t *s = pg->sops;s < (pg->sops + pg->order);s++){
77 double v[3];
78 applySymmetryOperation(s, ev, v);
79
80 for(f = 0;f < gel;f++){
81 if(ge[f].n == elements[i].n && ge[f].m == elements[i].m && 0 == strncmp(ge[f].name, elements[i].name, sizeof(ge[f].name)) && vequal(ge[f].v, v, thresholds->permutation)){
82 break;
83 }
84 }
85 if(f == gel){
86 memcpy(&ge[gel],&elements[i],sizeof(msym_element_t));
87 ge[gel].id = NULL;
88 vcopy(v, ge[gel].v);
89 aes->elements[aes->length++] = &ge[gel++];
90 }
91 }
92
93 if(!aes->length || (pg->order % aes->length != 0)){
94 msymSetErrorDetails("Equivalence set length (%d) not a divisor of point group order (%d)",pg->order);
95 ret = MSYM_INVALID_EQUIVALENCE_SET;
96 goto err;
97 }
98
99 aes->elements = realloc(aes->elements,sizeof(msym_element_t*[aes->length]));
100 }
101
102 msym_element_t *geo = ge;
103 ge = realloc(ge,sizeof(msym_element_t[gel]));
104 ges = realloc(ges,sizeof(msym_equivalence_set_t[gesl]) + sizeof(msym_element_t *[gel]));
105
106 msym_element_t **ep = (msym_element_t **) &ges[gesl];
107 for(int i = 0;i < gesl;i++){
108 msym_element_t **tep = ep;
109 for(int j = 0;j < ges[i].length;j++){
110 *ep = ges[i].elements[j] - geo + ge;
111 ep++;
112 }
113 free(ges[i].elements);
114 ges[i].elements = tep;
115 }
116
117 *glength = gel;
118 *gelements = ge;
119 *es = ges;
120 *esl = gesl;
121 return ret;
122
123 err:
124 free(ge);
125 for(int i = 0; i < gesl;i++) free(ges[i].elements);
126 free(ges);
127 return ret;
128 }
129
splitPointGroupEquivalenceSets(msym_point_group_t * pg,int esl,msym_equivalence_set_t es[esl],int * sesl,msym_equivalence_set_t ** ses,msym_thresholds_t * thresholds)130 msym_error_t splitPointGroupEquivalenceSets(msym_point_group_t *pg, int esl, msym_equivalence_set_t es[esl], int *sesl, msym_equivalence_set_t **ses, msym_thresholds_t *thresholds){
131 msym_error_t ret = MSYM_SUCCESS;
132 int length = 0, gesl = 0;
133 for(int i = 0;i < esl;i++) length += es[i].length;
134 msym_equivalence_set_t *ges = NULL;
135 msym_element_t **pelements = calloc(length,sizeof(msym_element_t*));
136 msym_element_t **ep = (msym_element_t **) &es[esl];
137
138 for(int i = 0; i < esl;i++){
139 msym_equivalence_set_t *pes = NULL;
140 int pesl = 0;
141 if(MSYM_SUCCESS != (ret = partitionPointGroupEquivalenceSets(pg, es[i].length, es[i].elements, es[i].elements - ep + pelements, &pesl, &pes, thresholds))) goto err;
142 ges = realloc(ges, sizeof(msym_equivalence_set_t[gesl+pesl]));
143 memcpy(&ges[gesl], pes, sizeof(msym_equivalence_set_t[pesl]));
144 free(pes);
145 gesl += pesl;
146 }
147
148 ges = realloc(ges, sizeof(msym_equivalence_set_t[gesl]) + sizeof(msym_element_t *[length]));
149 ep = (msym_element_t **) &ges[gesl];
150 memcpy(ep, pelements, sizeof(msym_element_t *[length]));
151
152 for(int i = 0;i < gesl;i++){
153 ges[i].elements = ep;
154 ep += ges[i].length;
155 }
156
157 *sesl = gesl;
158 *ses = ges;
159
160 free(pelements);
161 return ret;
162 err:
163 free(ges);
164 free(pelements);
165 return ret;
166 }
167
findPointGroupEquivalenceSets(msym_point_group_t * pg,int length,msym_element_t * elements[length],int * esl,msym_equivalence_set_t ** es,msym_thresholds_t * thresholds)168 msym_error_t findPointGroupEquivalenceSets(msym_point_group_t *pg, int length, msym_element_t *elements[length], int *esl, msym_equivalence_set_t **es, msym_thresholds_t *thresholds){
169 msym_error_t ret = MSYM_SUCCESS;
170 msym_equivalence_set_t *ges = NULL;
171 msym_element_t **pelements = calloc(length,sizeof(msym_element_t*));
172 int gesl = 0;
173 if(MSYM_SUCCESS != (ret = partitionPointGroupEquivalenceSets(pg, length, elements, pelements, &gesl, &ges, thresholds))) goto err;
174
175 ges = realloc(ges,sizeof(msym_equivalence_set_t[gesl]) + sizeof(msym_element_t *[length]));
176 msym_element_t **ep = (msym_element_t **) &ges[gesl];
177 msym_element_t **epo = ep;
178 memcpy(ep, pelements, sizeof(msym_element_t *[length]));
179 for(int i = 0;i < gesl;i++){
180 if(ep > epo + length){
181 msymSetErrorDetails("Equivalence set pointer (%ld) extends beyond number of elements (%d)",ep-epo,length);
182 ret = MSYM_INVALID_EQUIVALENCE_SET;
183 goto err;
184 }
185 ges[i].elements = ep;
186 ep += ges[i].length;
187 }
188
189 *es = ges;
190 *esl = gesl;
191
192 free(pelements);
193 return ret;
194 err:
195 free(ges);
196 free(pelements);
197 return ret;
198
199 }
200
partitionPointGroupEquivalenceSets(msym_point_group_t * pg,int length,msym_element_t * elements[length],msym_element_t * pelements[length],int * esl,msym_equivalence_set_t ** es,msym_thresholds_t * thresholds)201 msym_error_t partitionPointGroupEquivalenceSets(msym_point_group_t *pg, int length, msym_element_t *elements[length], msym_element_t *pelements[length], int *esl, msym_equivalence_set_t **es, msym_thresholds_t *thresholds){
202 msym_error_t ret = MSYM_SUCCESS;
203 msym_equivalence_set_t *ges = calloc(length,sizeof(msym_equivalence_set_t));
204 int *eqi = malloc(sizeof(int[length]));
205 memset(eqi,-1,sizeof(int[length]));
206 int gesl = 0, pelementsl = 0;
207 for(int i = 0;i < length;i++){
208 if(eqi[i] >= 0) continue;
209 if(pelementsl >= length){
210 msymSetErrorDetails("Size of equivalence sets (%d) larger than number of elements (%d)",pelementsl,length);
211 ret = MSYM_INVALID_EQUIVALENCE_SET;
212 goto err;
213 }
214
215 msym_equivalence_set_t *aes = &ges[gesl++];
216 aes->elements = &pelements[pelementsl];
217 for(msym_symmetry_operation_t *s = pg->sops;s < (pg->sops + pg->order);s++){
218 double v[3];
219 int f;
220 applySymmetryOperation(s, elements[i]->v, v);
221 for(f = 0;f < length;f++){
222 if(elements[f]->n == elements[i]->n && elements[f]->m == elements[i]->m && 0 == strncmp(elements[f]->name, elements[i]->name, sizeof(elements[f]->name)) && vequal(elements[f]->v, v, thresholds->permutation)){
223 break;
224 }
225 }
226
227 if(f < length && eqi[f] >= 0 && eqi[f] != gesl-1){
228 char buf[64];
229 symmetryOperationName(s, 64, buf);
230 msymSetErrorDetails("Symmetry operation %s on element %d yeilded element (%d) in two diffenrent equivalence sets (%d and %d)",buf,i,f,eqi[f],gesl-1);
231 ret = MSYM_INVALID_EQUIVALENCE_SET;
232 goto err;
233 } else if(f < length && eqi[f] == gesl-1){
234 //clean_debug_printf("element[%d] %s belongs to equivalence set %d, but already added\n",f,elements[f]->name, eqi[f]);
235 } else if(f < length){
236 eqi[f] = gesl - 1;
237 aes->elements[aes->length++] = elements[f];
238 //clean_debug_printf("element[%d] %s belongs to equivalence set %d, adding\n",f,elements[f]->name, eqi[f]);
239 } else {
240 char buf[64];
241 symmetryOperationName(s, 64, buf);
242 msymSetErrorDetails("Cannot find permutation for %s when determining equivalence set from point group %s",buf,pg->name);
243 ret = MSYM_INVALID_EQUIVALENCE_SET;
244 goto err;
245 }
246 }
247 //printf("generated equivalance set %d of length %d\n",gesl-1,aes->length);
248 pelementsl += aes->length;
249 }
250
251 if(pelementsl != length){
252 msymSetErrorDetails("Size of equivalence sets (%d) is not equal to number of elements (%d)",pelementsl,length);
253 ret = MSYM_INVALID_EQUIVALENCE_SET;
254 goto err;
255 }
256
257 *es = ges;
258 *esl = gesl;
259
260 free(eqi);
261 return ret;
262 err:
263 free(eqi);
264 free(ges);
265 return ret;
266
267 }
268
findEquivalenceSets(int length,msym_element_t * elements[length],msym_geometry_t g,int * esl,msym_equivalence_set_t ** es,msym_thresholds_t * thresholds)269 msym_error_t findEquivalenceSets(int length, msym_element_t *elements[length], msym_geometry_t g, int *esl, msym_equivalence_set_t **es, msym_thresholds_t *thresholds) {
270 msym_error_t ret = MSYM_SUCCESS;
271 int sesl = 0;
272 msym_equivalence_set_t *ses = NULL;
273 msym_element_t **pelements = calloc(length,sizeof(msym_element_t *));
274
275 if(MSYM_SUCCESS != (ret = partitionEquivalenceSets(length, elements,pelements,g,&sesl,&ses,thresholds))) goto err;
276
277 if(sesl > 1){
278 for(int i = 0; i < sesl;i++){
279 int rsesl = 0;
280 msym_equivalence_set_t *rses = NULL;
281 if(MSYM_SUCCESS != (ret = partitionEquivalenceSets(ses[i].length, ses[i].elements,ses[i].elements,g, &rsesl,&rses,thresholds))) goto err;
282
283 if(rsesl > 1){
284 ses[i].elements = rses[0].elements;
285 ses[i].length = rses[0].length;
286 ses = realloc(ses, sizeof(msym_equivalence_set_t[sesl+rsesl-1]));
287 memcpy(&ses[sesl], &rses[1], sizeof(msym_equivalence_set_t[rsesl-1]));
288 sesl += rsesl-1;
289 i--;
290 }
291 free(rses);
292 }
293 }
294
295 ses = realloc(ses, sizeof(msym_equivalence_set_t[sesl]) + sizeof(msym_element_t *[length]));
296 msym_element_t **ep = (msym_element_t **) &ses[sesl];
297
298 for(int i = 0;i < sesl;i++){
299 memcpy(ep, ses[i].elements, sizeof(msym_element_t *[ses[i].length]));
300 ses[i].elements = ep;
301 ep += ses[i].length;
302 }
303
304 *esl = sesl;
305 *es = ses;
306 free(pelements);
307 return ret;
308 err:
309 free(pelements);
310 free(ses);
311 return ret;
312
313 }
314
315
partitionEquivalenceSets(int length,msym_element_t * elements[length],msym_element_t * pelements[length],msym_geometry_t g,int * esl,msym_equivalence_set_t ** es,msym_thresholds_t * thresholds)316 msym_error_t partitionEquivalenceSets(int length, msym_element_t *elements[length], msym_element_t *pelements[length], msym_geometry_t g, int *esl, msym_equivalence_set_t **es, msym_thresholds_t *thresholds) {
317
318 int ns = 0, gd = geometryDegenerate(g);
319 double *e = calloc(length,sizeof(double));
320 double *s = calloc(length,sizeof(double));
321
322 int *sp = calloc(length,sizeof(int)); //set partition
323 int *ss = calloc(length,sizeof(int)); //set size
324
325 double (*ev)[3] = calloc(length,sizeof(double[3]));
326 double (*ep)[3] = calloc(length,sizeof(double[3]));
327
328 double (*vec)[3] = calloc(length, sizeof(double[3]));
329 double *m = calloc(length, sizeof(double));
330
331 for(int i = 0;i < length;i++){
332 vcopy(elements[i]->v, vec[i]);
333 m[i] = elements[i]->m;
334 }
335
336 for(int i=0; i < length; i++){
337 for(int j = i+1; j < length;j++){
338 double w = m[i]*m[j]/(m[i]+m[j]);
339 double dist;
340 double v[3];
341 double proji[3], projj[3];
342
343 vnorm2(vec[i],v);
344 vproj_plane(vec[j], v, proji);
345 vscale(w, proji, proji);
346 vadd(proji,ep[i],ep[i]);
347
348 vnorm2(vec[j],v);
349 vproj_plane(vec[i], v, projj);
350 vscale(w, projj, projj);
351 vadd(projj,ep[j],ep[j]);
352
353 vsub(vec[j],vec[i],v);
354
355 dist = vabs(v);
356
357 vscale(w/dist,v,v);
358
359 vadd(v,ev[i],ev[i]);
360 vsub(ev[j],v,ev[j]);
361
362 double dij = w*dist; //This is sqrt(I) for a diatomic molecule along an axis perpendicular to the bond with O at center of mass.
363 e[i] += dij;
364 e[j] += dij;
365
366 s[i] += SQR(dij);
367 s[j] += SQR(dij);
368 }
369 vsub(vec[i],ev[i],ev[i]);
370
371 }
372
373 for(int i = 0; i < length; i++){
374
375 double v[3];
376 double w = m[i]/2.0;
377 double dist = vabs(elements[i]->v);
378 double dii = w*dist;
379 vscale(w,elements[i]->v,v);
380 vsub(ev[i],v,ev[i]);
381
382 // Plane projection can't really differentiate certain types of structures when we add the initial vector,
383 // but not doing so will result in huge cancellation errors on degenerate point groups,
384 // also large masses will mess up the eq check when this is 0.
385 if(gd) vadd(ep[i],v,ep[i]);
386
387 e[i] += dii;
388 s[i] += SQR(dii);
389 }
390 for(int i = 0; i < length; i++){
391 if(e[i] >= 0.0){
392 sp[i] = i;
393 for(int j = i+1; j < length;j++){
394 if(e[j] >= 0.0){
395 double vabsevi = vabs(ev[i]), vabsevj = vabs(ev[j]), vabsepi = vabs(ep[i]), vabsepj = vabs(ep[j]);
396 double eep = 0.0, eev = fabs((vabsevi)-(vabsevj))/((vabsevi)+(vabsevj)), ee = fabs((e[i])-(e[j]))/((e[i])+(e[j])), es = fabs((s[i])-(s[j]))/((s[i])+(s[j]));
397
398 if(!(vabsepi < thresholds->zero && vabsepj < thresholds->zero)){
399 eep = fabs((vabsepi)-(vabsepj))/((vabsepi)+(vabsepj));
400 }
401
402 double max = fmax(eev,fmax(eep,fmax(ee, es)));
403
404 if(max < thresholds->equivalence && elements[i]->n == elements[j]->n){
405 e[j] = max > 0.0 ? -max : -1.0;
406 sp[j] = i;
407 }
408 }
409 }
410 e[i] = -1.0;
411 }
412 }
413
414 for(int i = 0; i < length;i++){
415 int j = sp[i];
416 ns += (ss[j] == 0);
417 ss[j]++;
418 }
419
420 msym_equivalence_set_t *eqs = calloc(ns,sizeof(msym_equivalence_set_t));
421 msym_element_t **lelements = elements;
422 msym_element_t **pe = pelements;
423
424 if(elements == pelements){
425 lelements = malloc(sizeof(msym_element_t *[length]));
426 memcpy(lelements, elements, sizeof(msym_element_t *[length]));
427 }
428
429 for(int i = 0, ni = 0; i < length;i++){
430 if(ss[i] > 0){
431 int ei = 0;
432 eqs[ni].elements = pe;
433 eqs[ni].length = ss[i];
434 for(int j = 0; j < length;j++){
435 if(sp[j] == i){
436 double err = (e[j] > -1.0) ? fabs(e[j]) : 0.0;
437 eqs[ni].err = fmax(eqs[ni].err,err);
438 eqs[ni].elements[ei++] = lelements[j];
439 }
440 }
441 pe += ss[i];
442 ni++;
443 }
444 }
445
446 if(elements == pelements){
447 free(lelements);
448 }
449 free(m);
450 free(vec);
451 free(s);
452 free(e);
453 free(sp);
454 free(ss);
455 free(ev);
456 free(ep);
457 *es = eqs;
458 *esl = ns;
459 return MSYM_SUCCESS;
460 }
461