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