1 //
2 // symmetrize.c
3 // Symmetry
4 //
5 // Created by Marcus Johansson on 04/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 <float.h>
15
16 #include "symmetrize.h"
17 #include "linalg.h"
18
19 #include "debug.h"
20
21 #define SQR(x) ((x)*(x))
22
23 /* This is a projection into the fully symmetric space.
24 * A little more computation than if we just recreate it from one atom,
25 * but it is independant of the chosen atom and we can get the size
26 * of the fully symmetric component.
27 * The sizes of the individual equivalence sets are rather small anyways.
28 */
29
symmetrizeElements(msym_point_group_t * pg,int esl,msym_equivalence_set_t * es,msym_permutation_t ** perm,msym_thresholds_t * thresholds,double * err)30 msym_error_t symmetrizeElements(msym_point_group_t *pg, int esl, msym_equivalence_set_t *es, msym_permutation_t **perm, msym_thresholds_t *thresholds, double *err){
31 msym_error_t ret = MSYM_SUCCESS;
32 double e = 0.0;
33 double (*v)[3] = malloc(sizeof(double[pg->order][3]));
34 for(int i = 0; i < esl;i++){
35 if(es[i].length > pg->order){
36 ret = MSYM_SYMMETRIZATION_ERROR;
37 msymSetErrorDetails("Equivalence set (%d elements) larger than order of point group (%d)",es[i].length,pg->order);
38 goto err;
39 }
40 memset(v, 0, sizeof(double[pg->order][3]));
41 for(int j = 0; j < pg->order;j++){
42 for(int k = 0; k < es[i].length;k++){
43 int p = perm[i][j].p[k];
44 double sv[3];
45 applySymmetryOperation(&pg->sops[j], es[i].elements[k]->v, sv);
46 vadd(sv, v[p], v[p]);
47 }
48 }
49 double sl = 0.0, ol = 0.0;
50 for(int j = 0; j < es[i].length;j++){
51 ol += vdot(es[i].elements[j]->v,es[i].elements[j]->v);
52 sl += vdot(v[j],v[j]);
53 vscale(1.0/((double)pg->order), v[j], es[i].elements[j]->v);
54 }
55 sl /= SQR((double)pg->order);
56 if(!(es[i].length == 1 && ol <= thresholds->zero)) e += (ol-sl)/ol; //e = fmax(e,(ol-sl)/ol);
57 }
58
59 *err = sqrt(fmax(e,0.0)); //should never be < 0, but it's a dumb way to die
60 err:
61 free(v);
62 return ret;
63 }
64
symmetrizeWavefunctions(msym_point_group_t * pg,int srsl,msym_subrepresentation_space_t * srs,int * span,int basisl,msym_basis_function_t basis[basisl],double wf[basisl][basisl],double symwf[basisl][basisl],int specieso[basisl],msym_partner_function_t pfo[basisl])65 msym_error_t symmetrizeWavefunctions(msym_point_group_t *pg, int srsl, msym_subrepresentation_space_t *srs, int *span, int basisl, msym_basis_function_t basis[basisl], double wf[basisl][basisl], double symwf[basisl][basisl], int specieso[basisl], msym_partner_function_t pfo[basisl]){
66 msym_error_t ret = MSYM_SUCCESS;
67
68 if(srsl != pg->ct->d){
69 ret = MSYM_SYMMETRIZATION_ERROR;
70 msymSetErrorDetails("Unexpected subspace length (expected %d got %d)",pg->ct->d, srsl);
71 return ret;
72 }
73
74 int *ispan = calloc(pg->ct->d,sizeof(*ispan));
75 int *species = (NULL == specieso ? malloc(sizeof(int[basisl])) : specieso);
76
77
78 memset(species,0,sizeof(int[basisl]));
79 if(NULL != pfo) memset(pfo,0,sizeof(msym_partner_function_t[basisl]));
80
81 int md = 1;
82 //could deduce from pg type but can't be bothered
83 for(int k = 0;k < pg->ct->d;k++) md = (md > pg->ct->s[k].d ? md : pg->ct->s[k].d);
84 int (*pf)[md] = calloc(basisl+1,sizeof(*pf));
85
86 int psalcl = 0;
87
88 int *psalck = calloc(pg->ct->d, sizeof(*psalck));
89
90 for(int i = 0;i < srsl;i++){
91 psalck[i] = psalcl;
92 psalcl += srs[i].salcl;
93 }
94 /* This is a bit of overkill since we only need the sign of the component
95 * We could just as well have used a basis change, but the sign in degenerate
96 * representations can't be used for partner function determination,
97 * and this is a sparce matrix
98 */
99 double (*psalc)[basisl][psalcl] = calloc(md+1,sizeof(*psalc));
100 double (*bfd)[md] = calloc(basisl+1, sizeof(*bfd));
101 double *dmpf = bfd[basisl];
102
103 /* Determine salc components, and build information vectors (e.g. indexing/offsets/irreps) */
104 for(int o = 0;o < basisl;o++){
105 double mcomp = -1.0;
106 int psalci = 0;
107 for(int k = 0;k < srsl;k++){
108 double mabs = 0.0;
109 for(int s = 0;s < srs[k].salcl;s++){
110 msym_salc_t *salc = &srs[k].salc[s];
111 double (*space)[salc->fl] = (double (*)[salc->fl]) salc->pf;
112 double psalcabs = 0.0;
113 for(int d = 0;d < salc->d;d++){
114 double c = 0.0;
115 for(int j = 0; j < salc->fl;j++){
116 c += wf[o][salc->f[j] - basis]*space[d][j];
117 }
118 double c2 = SQR(c);
119 psalc[d][o][psalci] = c; //Won't thrash the cache too bad since d is small
120 psalcabs += c2;
121 bfd[o][d] += c2;
122 }
123 mabs += psalcabs;
124 psalc[md][o][psalci] = sqrt(psalcabs);
125 psalci++;
126 }
127 if(mabs > mcomp){
128 species[o] = k;
129 mcomp = mabs;
130 }
131 }
132 ispan[species[o]]++;
133 }
134
135 for(int k = 0;k < pg->ct->d;k++){
136 if(ispan[k] != span[k]*pg->ct->s[k].d){
137 msymSetErrorDetails("Projected orbitals do not span the expected irredicible representations. Expected %d%s, got %d",span[k],pg->ct->s[k].name,ispan[k]);
138 ret = MSYM_SYMMETRIZATION_ERROR;
139 goto err;
140 }
141 }
142
143 /* Find parner functions */
144 for(int o = 0;o < basisl;o++){
145 int ko = species[o], dim = pg->ct->s[ko].d;
146
147 struct _fpf {int i; int j;} fpf = {.i = 0, .j = 0};
148
149 for(int i = 1;i < md;i++){
150 pf[o][i] = -1;
151 pf[basisl][i] = -1;
152 };
153
154 if(dim <= 1) continue;
155
156 /* check if this functions has alredy been assigned partners */
157 for(int i = 0;i < o && !fpf.j;i++){
158 for(int j = 1;j < md;j++){
159 if(pf[i][j] == o){
160 fpf.i = i;
161 fpf.j = j;
162 break;
163 }
164 }
165 }
166
167 if(fpf.j){
168 for(int i = 1; i < fpf.j;i++){
169 pf[pf[fpf.i][i]][0]--;
170 }
171 pf[o][0] -= fpf.j;
172 continue;
173 }
174
175 for(int i = 0;i < md;i++){dmpf[i] = DBL_MAX;}
176
177 for(int po = 0; po < basisl;po++){
178 if(species[po] != ko || o == po || abs(pf[po][0]) == dim - 1) continue;
179 double c = 0, mc = 0;
180 /* length of v1-v2 */
181 for(int i = 0;i < psalcl;i++){
182 double sub = psalc[md][o][i] - psalc[md][po][i];
183 c += SQR(sub);
184 }
185 c = sqrt(c);
186
187 /* find the <dim> smallest diffs */
188 int mic = 0;
189 for(int i = 1;i < dim;i++){
190 double diff = fabs(dmpf[i] - c);
191 if(c < dmpf[i] && (diff > mc)){
192 mic = i;
193 mc = diff;
194 }
195 }
196 if(mic > 0){
197 dmpf[mic] = c;
198 pf[o][mic] = po;
199 pf[basisl][mic] = po;
200 }
201 }
202
203 for(int i = 1;i < dim;i++){
204 pf[o][0] += pf[basisl][i] > 0;
205 }
206 }
207
208 /* verify that we have partners for everything */
209 for(int o = 0;o < basisl;o++){
210 int dim = pg->ct->s[species[o]].d;
211 if(abs(pf[o][0])+1 != dim){
212 for(int i = 0;i < md;i++) clean_debug_printf("%d = %d\n",i,pf[o][i]);
213 msymSetErrorDetails("Unexpected number of partner functions for wave function %d in %s (expected %d got %d)", o, pg->ct->s[species[o]].name, dim, abs(pf[o][0])+1);
214 ret = MSYM_SYMMETRIZATION_ERROR;
215 goto err;
216 }
217
218 for(int i = 0;pf[o][0] >= 0 && i < dim;i++){
219 if(pf[o][i] == -1){
220 msymSetErrorDetails("Could not determine partner function %d of wave function %d",i, o);
221 ret = MSYM_SYMMETRIZATION_ERROR;
222 goto err;
223 }
224 }
225 }
226
227
228 memset(symwf,0,sizeof(double[basisl][basisl]));
229
230 for(int o = 0;o < basisl;o++){
231 int k = species[o];
232 int dim = pg->ct->s[k].d;
233
234 if(pf[o][0] < 0) continue;
235
236 pf[o][0] = o;
237 for(int i = 0;i < dim;i++) pf[basisl][i] = -1;
238
239 /* Get the unique dimensions for each partner function in which they have the largest component.
240 * This is only needed when the symmetry is really broken but the degenerate functions can be averaged,
241 * but it also keeps the ordering intact.
242 * This is could be improved with a best match algorithm */
243 for(int i = 0;i < dim;i++){
244 double cmax = 0.0;
245 for(int d = 0;d < dim;d++){
246 double c = bfd[pf[o][i]][d]; //component of i:th partner in dimension d
247 if(c > cmax){
248 int found = 0;
249 for(int j = 0;j < i;j++){
250 if(pf[basisl][j] == d){
251 found = 1;
252 break;
253 }
254 }
255 if(!found){
256 pf[basisl][i] = d;
257 cmax = c;
258 }
259 }
260 }
261 }
262
263 /*for(int i = 0;i < dim;i++){
264 clean_debug_printf("partner function %d has maximum component in dimension %d\n",i,pf[basisl][i]);
265 }*/
266
267 /* calculate average component in each salc subspace and rotate onto the partner functions with largest component */
268 for(int s = 0;s < srs[k].salcl;s++){
269 int psalci = psalck[k] + s, pfmin = 0;
270 double avg = 0;
271
272
273 for(int d = 0;d < dim;d++){
274 int wfi = pf[o][d];
275 avg += psalc[md][wfi][psalci];
276 if(pf[basisl][d] == 0) pfmin = wfi;
277 }
278
279 avg /= dim;
280
281 //printf("average component in salc %d(%s) for wf %d = %lf\n",psalci,pg->ct->s[k].name,o,avg);
282
283 msym_salc_t *salc = &srs[k].salc[s];
284 double (*space)[salc->fl] = (double (*)[salc->fl]) salc->pf;
285
286 for(int d = 0; d < dim;d++){
287 int wfi = pf[o][d], di = pf[basisl][d];
288 /* use the sign of the projection onto the largest component */
289 double c = copysign(avg,psalc[di][wfi][psalci]);
290 for(int j = 0; j < salc->fl;j++){
291 symwf[wfi][salc->f[j] - basis] += c*space[di][j];
292 }
293 if(NULL != pfo){
294 pfo[wfi].d = di;
295 pfo[wfi].i = pfmin;
296 }
297 }
298 }
299 }
300
301 err:
302
303 if(species != specieso) free(species);
304 free(ispan);
305 free(pf);
306 free(psalc);
307 free(bfd);
308 free(psalck);
309
310 return ret;
311 }
312
symmetrizeTranslation(msym_point_group_t * pg,msym_equivalence_set_t * es,msym_permutation_t * perm,int pi,double translation[3])313 msym_error_t symmetrizeTranslation(msym_point_group_t *pg, msym_equivalence_set_t *es, msym_permutation_t *perm, int pi, double translation[3]){
314 msym_error_t ret = MSYM_SUCCESS;
315 double (*v)[3] = calloc(es->length,sizeof(double[3]));
316
317 for(int j = 0; j < pg->order;j++){
318 int p = perm[j].p[pi];
319 double stranslation[3];
320 applySymmetryOperation(&pg->sops[j], translation, stranslation);
321 vadd(stranslation, v[p], v[p]);
322 }
323
324 double scale = ((double)es->length)/pg->order;
325
326 for(int i = 0;i < es->length;i++){
327 vscale(scale, v[i], v[i]);
328 vadd(es->elements[i]->v,v[i],es->elements[i]->v);
329 }
330
331 //err:
332 free(v);
333 return ret;
334 }
335
336
337
338