1 //
2 //  symop.c
3 //  Symmetry
4 //
5 //  Created by Marcus Johansson on 30/10/14.
6 //  Copyright (c) 2014 Marcus Johansson.
7 //
8 //  Distributed under the MIT License ( See LICENSE file or copy at http://opensource.org/licenses/MIT )
9 //
10 
11 #include <math.h>
12 #include <stdlib.h>
13 #include "msym.h"
14 #include "context.h"
15 #include "symop.h"
16 #include "linalg.h"
17 
18 #include "debug.h"
19 
20 #ifndef M_PI
21 #define M_PI 3.14159265358979323846264338327950288419716939937510582
22 #endif
23 
24 int igcd(int a, int b);
25 
symopPow(msym_symmetry_operation_t * A,int pow,msym_symmetry_operation_t * O)26 void symopPow(msym_symmetry_operation_t *A, int pow, msym_symmetry_operation_t *O){
27     double origo[3] = {0.0,0.0,0.0};
28     int apow, gcd;
29     O->power = 1;
30     O->orientation = A->orientation;
31     switch (A->type){
32         case IDENTITY :
33             O->type = IDENTITY;
34             O->order = 0;
35             vcopy(origo,O->v);
36             break;
37         case REFLECTION :
38         case INVERSION :
39             if(pow % 2 == 0){
40                 O->type = IDENTITY;
41                 O->order = 0;
42                 vcopy(origo,O->v);
43             } else {
44                 O->type = A->type;
45                 O->order = 0;
46                 vcopy(A->v,O->v);
47             }
48             break;
49         case PROPER_ROTATION :
50             // A->power = 0 is assumed to be uninitialized and set to 1
51             apow = A->power == 0 ? pow % A->order : (A->power*pow) % A->order;
52             gcd = igcd(apow, A->order);
53             if(apow == 0){
54                 O->type = IDENTITY;
55                 O->order = 0;
56                 vcopy(origo,O->v);
57             } else {
58                 O->type = PROPER_ROTATION;
59                 O->order = A->order/gcd;
60                 O->power = apow/gcd;
61                 vcopy(A->v,O->v);
62             }
63             break;
64         case IMPROPER_ROTATION :
65             A->type = PROPER_ROTATION;
66             symopPow(A, pow, O);
67             A->type = IMPROPER_ROTATION;
68             apow = A->power == 0 ? pow % (2*A->order) : A->power*pow % (2*A->order);
69             if(O->type == IDENTITY && pow == 0){} // Do nothing
70             else if(O->type == IDENTITY && A->order % 2 == 1 && apow == A->order){
71                 O->type = REFLECTION;
72                 O->order = 0;
73                 vcopy(A->v,O->v);
74             } else {
75                 if (apow > A->order && A->order % 2 == 1 && apow % 2 == 1){
76                 //if (apow > O->order && O->order % 2 == 1 && apow % 2 == 1){
77                     O->power += A->order;
78                     O->power %= 2*O->order;
79                 }
80                 if(apow % 2 == 1){
81                     O->type = IMPROPER_ROTATION;
82                 }
83                 if(O->type == IMPROPER_ROTATION && O->order == 2){
84                     O->type = INVERSION;
85                     O->power = 1;
86                 }
87             }
88         default: break;
89     }
90 }
91 
92 #define PHI 1.618033988749894848204586834
93 
symmetryOperationCharacter(msym_symmetry_operation_t * sop,msym_basis_function_t * f)94 double symmetryOperationCharacter(msym_symmetry_operation_t *sop, msym_basis_function_t *f){
95     double c = 0.0;
96     switch (f->type) {
97         case MSYM_BASIS_TYPE_REAL_SPHERICAL_HARMONIC:
98             c = symmetryOperationYCharacter(sop,f->f.rsh.l);
99             break;
100         case MSYM_BASIS_TYPE_CARTESIAN :
101             c = symmetryOperationCartesianCharacter(sop);
102             break;
103         default:
104             break;
105     }
106     return c;
107 }
108 
symmetryOperationCartesianCharacter(msym_symmetry_operation_t * sop)109 double symmetryOperationCartesianCharacter(msym_symmetry_operation_t *sop){
110     double x = 0.0;
111     switch (sop->type) {
112         case IDENTITY          : x = 3; break;
113         case INVERSION         : x = -3; break;
114         case REFLECTION        : x = 1; break;
115         case PROPER_ROTATION   : x = 2*cos(sop->power*2*M_PI/sop->order) + 1; break;
116         case IMPROPER_ROTATION : x = 2*cos(sop->power*2*M_PI/sop->order) - 1; break;
117         default:
118             break;
119     }
120     return x;
121 }
122 
123 /* character of spherical harmonics basis */
symmetryOperationYCharacter(msym_symmetry_operation_t * sop,int l)124 double symmetryOperationYCharacter(msym_symmetry_operation_t *sop, int l){
125     double x = 0.0, a = sop->power*2*M_PI/sop->order;
126     switch (sop->type) {
127         case IDENTITY          : x = 2*l+1; break;
128         case INVERSION         : x = (2*l+1)*(1 - ((l & 1) << 1)); break;
129         case REFLECTION        : x = 1; break;
130         case PROPER_ROTATION   : x = sin((l+0.5)*a)/sin(a/2); break;
131         case IMPROPER_ROTATION : x = cos((l+0.5)*a)/cos(a/2); break;
132         default:
133             break;
134     }
135     return x;
136 
137 
138 }
139 
symmetryOperationName(msym_symmetry_operation_t * sop,int l,char buf[l])140 void symmetryOperationName(msym_symmetry_operation_t* sop, int l, char buf[l]){
141     char *rn = "";
142     char *cn = "";
143     switch(sop->orientation){
144         case HORIZONTAL : rn = "h"; break;
145         case VERTICAL   : rn = "v"; cn = "'"; break;
146         case DIHEDRAL   : rn = "d"; cn = "''"; break;
147         default: break;
148     }
149     switch (sop->type) {
150         case PROPER_ROTATION   :
151             if(sop->order == 2) snprintf(buf, l, "C%d%s",sop->order,cn);
152             else snprintf(buf, l, "C%d%s^%d",sop->order,cn,sop->power);
153             break;
154         case IMPROPER_ROTATION : snprintf(buf, l, "S%d^%d",sop->order,sop->power); break;
155         case REFLECTION        : snprintf(buf, l, "R%s",rn); break;
156         case INVERSION         : snprintf(buf, l, "i"); break;
157         case IDENTITY          : snprintf(buf, l, "E"); break;
158         default                : snprintf(buf, l, "?"); break;
159     }
160 }
161 
symmetryOperationShortName(msym_symmetry_operation_t * sop,int l,char buf[l])162 void symmetryOperationShortName(msym_symmetry_operation_t* sop, int l, char buf[l]){
163     switch (sop->type) {
164         case PROPER_ROTATION   : snprintf(buf, l, "C%d",sop->order); break;
165         case IMPROPER_ROTATION : snprintf(buf, l, "S%d",sop->order); break;
166         case REFLECTION        : snprintf(buf, l, "R"); break;
167         case INVERSION         : snprintf(buf, l, "i"); break;
168         case IDENTITY          : snprintf(buf, l, "E"); break;
169         default                : snprintf(buf, l, "?"); break;
170     }
171 }
172 
applySymmetryOperation(msym_symmetry_operation_t * sop,double iv[3],double ov[3])173 void applySymmetryOperation(msym_symmetry_operation_t *sop,double iv[3], double ov[3]){
174     switch (sop->type) {
175         case PROPER_ROTATION   : {
176             if(sop->order){
177                 double theta = sop->power*2*M_PI/sop->order;
178                 vrotate(theta, iv, sop->v, ov);
179             } else {
180                 // Treat C0 as a sum over all Cn operations i.e. a projection
181                 vproj(iv, sop->v, ov);
182             }
183             break;
184         }
185         case IMPROPER_ROTATION : {
186             double theta = sop->power*2*M_PI/sop->order;
187             vrotate(theta, iv, sop->v, ov);
188             vreflect(ov, sop->v, ov);
189             break;
190         }
191         case REFLECTION        :
192             vreflect(iv, sop->v, ov);
193             break;
194         case INVERSION         :
195             ov[0] = -iv[0];
196             ov[1] = -iv[1];
197             ov[2] = -iv[2];
198             break;
199         case IDENTITY          :
200             ov[0] = iv[0];
201             ov[1] = iv[1];
202             ov[2] = iv[2];
203             break;
204         default : fprintf(stderr,"UNKNOWN OPERATION\n"); //should never happen so will not handle this error atm
205     }
206 
207 }
208 
invertSymmetryOperation(msym_symmetry_operation_t * sop,msym_symmetry_operation_t * isop)209 void invertSymmetryOperation(msym_symmetry_operation_t *sop, msym_symmetry_operation_t *isop){
210     copySymmetryOperation(isop, sop);
211     switch (sop->type) {
212         case PROPER_ROTATION   : {
213             isop->power = sop->order - sop->power;
214             break;
215         }
216         case IMPROPER_ROTATION : {
217             if(sop->order % 2 == 0) isop->power = sop->order - sop->power;
218             else isop->power = 2*sop->order - sop->power;
219             break;
220         }
221         case REFLECTION        :
222         case INVERSION         :
223         case IDENTITY          :
224             break;
225         default : fprintf(stderr,"UNKNOWN OPERATION\n"); //should never happen so will not handle this error atm
226     }
227 }
228 
symmetryOperationMatrix(msym_symmetry_operation_t * sop,double m[3][3])229 void symmetryOperationMatrix(msym_symmetry_operation_t *sop, double m[3][3]){
230     switch (sop->type) {
231         case PROPER_ROTATION   : {
232             if(sop->order){
233                 double theta = sop->power*2*M_PI/sop->order;
234                 mrotate(theta, sop->v, m);
235             } else {
236                 // Treat C0 as a sum over all Cn operations i.e. a projection
237                 kron2(3, 1, &sop->v, 1, 3, &sop->v, m);
238             }
239             break;
240         }
241         case IMPROPER_ROTATION : {
242             double theta = sop->power*2*M_PI/sop->order;
243             double m1[3][3], m2[3][3];
244             mrotate(theta, sop->v, m1);
245             mreflect(sop->v, m2);
246             mmmul(m2, m1, m);
247             break;
248         }
249         case REFLECTION        :
250             mreflect(sop->v, m);
251             break;
252         case INVERSION         : //TODO: fix this
253             for(int i = 0; i < 3; i++){
254                 for(int j = 0; j < 3; j++){
255                     m[i][j] = -(i == j);
256                 }
257             }
258             break;
259         case IDENTITY          :
260             for(int i = 0; i < 3; i++){
261                 for(int j = 0; j < 3; j++){
262                     m[i][j] = (i == j);
263                 }
264             }
265             break;
266         default : fprintf(stderr,"UNKNOWN OPERATION\n"); //should never happen so will not handle this error atm
267     }
268 }
269 
copySymmetryOperation(msym_symmetry_operation_t * dst,msym_symmetry_operation_t * src)270 void copySymmetryOperation(msym_symmetry_operation_t *dst, msym_symmetry_operation_t *src){
271     dst->type = src->type;
272     dst->order = src->order;
273     dst->power = src->power;
274     dst->cla = src->cla;
275     dst->orientation = src->orientation;
276     vcopy(src->v, dst->v);
277 }
278 
findSymmetryOperation(msym_symmetry_operation_t * sop,msym_symmetry_operation_t * sops,int l,msym_thresholds_t * thresholds)279 msym_symmetry_operation_t *findSymmetryOperation(msym_symmetry_operation_t *sop, msym_symmetry_operation_t *sops, int l, msym_thresholds_t *thresholds){
280     msym_symmetry_operation_t *fsop = NULL;
281     for(msym_symmetry_operation_t* s = sops; s < (sops + l);s++){
282         if((s->type == INVERSION && sop->type == INVERSION) || (s->type == IDENTITY && sop->type == IDENTITY)){
283             fsop = s;
284             break;
285         }
286 
287         else if(s->type == sop->type && ((sop->type != PROPER_ROTATION && sop->type != IMPROPER_ROTATION) || (s->order == sop->order && s->power == sop->power)) && vparallel(s->v, sop->v,thresholds->angle)){
288             fsop = s;
289             break;
290         }
291     }
292 
293     return fsop;
294 }
295 
printSymmetryOperation(msym_symmetry_operation_t * sop)296 void printSymmetryOperation(msym_symmetry_operation_t *sop){
297     char buf[12];
298     symmetryOperationName(sop, 12, buf);
299     if(sop->type == INVERSION || sop->type == IDENTITY)
300         clean_debug_printf("%s(%d)\n",buf,sop->cla);
301     else
302         clean_debug_printf("%s(%d) [%lf;%lf;%lf]\n",buf,sop->cla, sop->v[0],sop->v[1],sop->v[2]);
303 }
304 
igcd(int a,int b)305 int igcd(int a, int b) {
306     int c;
307     while (a) {
308         c = a;
309         a = b % a;
310         b = c;
311     }
312     return abs(b);
313 }
314