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