1 //
2 //  rsh.c
3 //  libmsym
4 //
5 //  Created by Marcus Johansson on 21/10/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 // Ref:
12 //
13 // J. Ivanic and K. Ruedenberg, "Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion",
14 // J. Phys. Chem., vol. 100, no. 15, pp. 6342-6347, 1996. http://pubs.acs.org/doi/pdf/10.1021/jp953350u
15 
16 #include <math.h>
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <float.h>
20 #include <string.h>
21 #include "rsh.h"
22 #include "linalg.h"
23 #include "symop.h"
24 
25 #define SQR(x) ((x)*(x))
26 
27 #ifndef M_SQRT2
28 #define M_SQRT2 1.41421356237309504880
29 #endif
30 
31 void rshSymmetryOperationRepresentation(msym_symmetry_operation_t *sops, int index, int l, rsh_representations_t *lts);
32 void rshCalculateUVWCoefficients(int l, int m1, int m2, double* u, double* v, double* w);
33 void rshRotationRepresentation(int index, int l, rsh_representations_t *lrs);
34 double rshRotationP(int index, int l, int i, int m1, int m2, rsh_representations_t *lrs);
35 double rshRotationU(int index, int l, int m1, int m2, rsh_representations_t *lrs);
36 double rshRotationV(int index, int l, int m1, int m2, rsh_representations_t *lrs);
37 double rshRotationW(int index, int l, int m1, int m2, rsh_representations_t *lrs);
38 
39 
generateRSHRepresentations(int sopsl,msym_symmetry_operation_t sops[sopsl],int lmax,rsh_representations_t * lrs)40 msym_error_t generateRSHRepresentations(int sopsl, msym_symmetry_operation_t sops[sopsl], int lmax, rsh_representations_t *lrs){
41     msym_error_t ret = MSYM_SUCCESS;
42     for(int l = 0;l <= lmax;l++){
43         if(lrs[l].d != 2*l+1){
44             ret = MSYM_INVALID_BASIS_FUNCTIONS;
45             msymSetErrorDetails("Invalid dimension of real spherical harmonic (expected %d, got %d)",2*l+1, lrs[l].d);
46             goto err;
47         }
48         for(int i = 0; i < sopsl;i++){
49             rshSymmetryOperationRepresentation(sops,i,l,lrs);
50         }
51     }
52 
53 err:
54     return ret;
55 }
56 
57 
58 
rshSymmetryOperationRepresentation(msym_symmetry_operation_t * sops,int index,int l,rsh_representations_t * lrs)59 void rshSymmetryOperationRepresentation(msym_symmetry_operation_t *sops, int index, int l, rsh_representations_t *lrs){
60     if(0 == l){
61         int d = lrs[0].d;
62         double (*st)[d][d] = lrs[0].t;
63         st[index][0][0] = 1;
64     } else if(1 == l){
65         int d = lrs[1].d;
66         double (*st)[d][d] = lrs[1].t;
67         double (*r)[d] = st[index];
68         double t[3][3];
69         symmetryOperationMatrix(&sops[index], t);
70         // x,y,z -> py, pz, px
71         r[0][0] = t[1][1];
72         r[0][1] = t[1][2];
73         r[0][2] = t[1][0];
74         r[1][0] = t[2][1];
75         r[1][1] = t[2][2];
76         r[1][2] = t[2][0];
77         r[2][0] = t[0][1];
78         r[2][1] = t[0][2];
79         r[2][2] = t[0][0];
80     } else {
81         int d = lrs[l].d;
82         double (*st)[d][d] = lrs[l].t;
83         double (*r)[d] = st[index];
84         switch (sops[index].type) {
85             case INVERSION:
86                 if(l & 1){
87                     memset(r,0,sizeof(double[d][d]));
88                     for(int i = 0;i < d;i++) r[i][i] = -1;
89                     break;
90                 } //fallthrough
91             case IDENTITY:
92                 mleye(d,r);
93                 break;
94             case REFLECTION:
95             case IMPROPER_ROTATION:
96             case PROPER_ROTATION:
97             default:
98                 rshRotationRepresentation(index, l, lrs);
99                 break;
100 
101         }
102     }
103 }
104 
rshRotationP(int index,int l,int i,int m1,int m2,rsh_representations_t * lrs)105 double rshRotationP(int index, int l, int i, int m1, int m2, rsh_representations_t *lrs){
106     int d1 = lrs[1].d, dl = lrs[l-1].d, ol = (dl - 1)/2;
107     double (*st1)[d1][d1] = lrs[1].t, (*stl)[dl][dl] = lrs[l-1].t;
108     double (*r1)[d1] = st1[index], (*rl)[dl] = stl[index];
109     double ret = 0;
110 
111     if (m2 == l) {
112         ret = r1[i + 1][2]*rl[m1 + ol][l - 1 + ol] - r1[i + 1][0]*rl[m1 + ol][1 - l + ol];
113     } else if (m2 == -l) {
114         ret = r1[i + 1][2]*rl[m1 + ol][1 - l + ol] + r1[i + 1][0]*rl[m1 + ol][l - 1 + ol];
115     } else {
116         ret = r1[i + 1][1]*rl[m1 + ol][m2 + ol];
117     }
118 
119     return ret;
120 }
121 
rshRotationU(int index,int l,int m1,int m2,rsh_representations_t * lrs)122 double rshRotationU(int index, int l, int m1, int m2, rsh_representations_t *lrs){
123     return rshRotationP(index, l, 0, m1, m2, lrs);
124 }
125 
rshRotationV(int index,int l,int m1,int m2,rsh_representations_t * lrs)126 double rshRotationV(int index, int l, int m1, int m2, rsh_representations_t *lrs){
127     double ret = 0;
128     if (m1 == 0) {
129         ret = rshRotationP(index, l, 1, 1, m2, lrs) + rshRotationP(index, l, -1, -1, m2, lrs);
130     } else if (m1 == 1){
131         ret = M_SQRT2*rshRotationP(index, l, 1, 0, m2, lrs);
132     } else if (m1 == -1){
133         ret = M_SQRT2*rshRotationP(index, l, -1, 0, m2, lrs);
134     } else if (m1 > 0){
135         ret = rshRotationP(index, l, 1, m1 - 1, m2, lrs) - rshRotationP(index,l, -1, -m1 + 1, m2, lrs);
136     } else {
137         ret = rshRotationP(index, l, 1, m1 + 1, m2, lrs) + rshRotationP(index,l, -1, -m1 - 1, m2, lrs);
138     }
139 
140     return ret;
141 }
142 
rshRotationW(int index,int l,int m1,int m2,rsh_representations_t * lrs)143 double rshRotationW(int index, int l, int m1, int m2, rsh_representations_t *lrs){
144     double ret = 0;
145     if (m1 > 0) {
146         ret = rshRotationP(index, l, 1, m1 + 1, m2, lrs) + rshRotationP(index, l, -1, -m1 - 1, m2, lrs);
147     } else {
148         ret = rshRotationP(index, l, 1, m1 - 1, m2, lrs) - rshRotationP(index, l, -1, -m1 + 1, m2, lrs);
149     }
150 
151     return ret;
152 }
153 
rshRotationRepresentation(int index,int l,rsh_representations_t * lrs)154 void rshRotationRepresentation(int index, int l, rsh_representations_t *lrs){
155     int d = lrs[l].d;
156     double (*st)[d][d] = lrs[l].t;
157     double (*r)[d] = st[index];
158     for (int m1 = -l; m1 <= l; m1++) {
159         for (int m2 = -l; m2 <= l; m2++) {
160             double u, v, w;
161             rshCalculateUVWCoefficients(l, m1, m2, &u, &v, &w);
162             if(u != 0){
163                 u *= rshRotationU(index, l, m1, m2, lrs);
164             }
165             if(v != 0){
166                 v *= rshRotationV(index, l, m1, m2, lrs);
167             }
168             if(w != 0){
169                 w *= rshRotationW(index, l, m1, m2, lrs);
170             }
171             r[m1 + l][m2 + l] = u + v + w;
172         }
173     }
174 }
175 
rshCalculateUVWCoefficients(int l,int m1,int m2,double * u,double * v,double * w)176 void rshCalculateUVWCoefficients(int l, int m1, int m2, double *u, double *v, double *w) {
177     if(m1 == 0){ // d = 1
178         if(abs(m2) == l){
179             *u =  sqrt(l/(4*l - 2.0));
180             *v = -0.5*sqrt((l - 1.0)/(2*l - 1.0));
181         } else {
182             double m22 = SQR(m2), l2 = SQR(l), l2m22 = l2 - m22;
183             *u =  sqrt(l2/l2m22);
184             *v = -0.5*sqrt(2*(l2 - l)/l2m22);
185         }
186         *w = 0;
187     } else { // d = 0
188         double am1 = abs(m1);
189         double div = (abs(m2) == l ? 2.0*l*(2.0*l - 1) : (l + m2)*(l - m2));
190         *u =  sqrt((l + m1)*(l - m1)/div);
191         *v =  0.5*sqrt((l + am1 - 1)*(l + am1)/div);
192         *w = -0.5*sqrt((l - am1 - 1)*(l - am1)/div);
193     }
194 }
195