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