1 //
2 // ico.cc --- implementation of icosahedral operations
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Edward Seidl <seidl@janed.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #include <util/misc/math.h>
29 #include <string.h>
30
31 #include <math/symmetry/pointgrp.h>
32
33 using namespace sc;
34
35 // these are the operations which make up T
36 static void
i_ops(SymRep * t1rep,SymRep * t2rep,SymRep * grep,SymRep * hrep)37 i_ops(SymRep *t1rep, SymRep *t2rep, SymRep *grep, SymRep *hrep)
38 {
39 int i;
40
41 // identity
42 t1rep[0].E();
43 t2rep[0].E();
44 grep[0].E();
45 hrep[0].E();
46
47 //
48 // 12 C5's
49 //
50 // first the 2 C5's about the z axis
51 t1rep[1].rotation(2.0*(double)M_PI/5.0);
52 t1rep[2].rotation(8.0*(double)M_PI/5.0);
53
54 t2rep[1] = t1rep[1].operate(t1rep[1]);
55 t2rep[2] = t1rep[2].operate(t1rep[2]);
56
57 grep[1].rotation(2.0*(double)M_PI/5.0);
58 grep[2].rotation(8.0*(double)M_PI/5.0);
59
60 hrep[1].rotation(2.0*(double)M_PI/5.0);
61 hrep[2].rotation(8.0*(double)M_PI/5.0);
62
63 // form rotation matrices for the C3 axis about the zx axis (these were
64 // taken from turbomole version 2, which claims they were sort of inherited
65 // from hondo
66 SymRep t1so(3);
67 SymRep gso(4);
68 SymRep hso(5);
69
70 double c2p5 = cos(2.0*(double)M_PI/5.0);
71 double s2p5 = sin(2.0*(double)M_PI/5.0);
72 double cosd = s2p5/((1.0-c2p5)*sqrt(3.0));
73 double cosd2 = cosd*cosd;
74 double sind2 = 1.0 - cosd2;
75 double sind = sqrt(sind2);
76
77 t1so[0][0] = 1.0 - 1.5*cosd2;
78 t1so[1][0] = 0.5*sqrt(3.0)*cosd;
79 t1so[2][0] = 1.5*cosd*sind;
80 t1so[0][1] = -0.5*sqrt(3.0)*cosd;
81 t1so[1][1] = -0.5;
82 t1so[2][1] = 0.5*sqrt(3.0)*sind;
83 t1so[0][2] = 1.5*cosd*sind;
84 t1so[1][2] = -0.5*sqrt(3.0)*sind;
85 t1so[2][2] = 1.0 - 1.5*sind2;
86
87 gso[0][0] = (3.0*sqrt(5.0)+5.0)/20.0;
88 gso[0][1] = cosd*sqrt(3.0)*(sqrt(5.0)-1.0)/4.0;
89 gso[0][2] = 3.0*sqrt(5.0)/10.0;
90 gso[0][3] = -sqrt(5.0-2.0*sqrt(5.0))*sqrt(5.0)/10.0;
91 gso[1][0] = -gso[0][1];
92 gso[1][1] = (1-sqrt(5.0))/4.0;
93 gso[1][2] = cosd*sqrt(3.0)/2.0;
94 gso[1][3] = cosd*sqrt(5-2*sqrt(5.0))*sqrt(3.0)/2.0;
95 gso[2][0] = gso[0][2];
96 gso[2][1] = -gso[1][2];
97 gso[2][2] = (5-3*sqrt(5.0))/20.0;
98 gso[2][3] = sqrt(5.0-2*sqrt(5.0))*(sqrt(5.0)+5)/20;
99 gso[3][0] = -gso[0][3];
100 gso[3][1] = gso[1][3];
101 gso[3][2] = -gso[2][3];
102 gso[3][3] = (sqrt(5.0)+1)/4.0;
103
104 hso[0][0] = -1.0/5.0;
105 hso[0][4] = sqrt(3.0)*(sqrt(5.0)+1)/10.0;
106 hso[0][3] = 3.0*cosd*(3.0*sqrt(5.0)-5.0)/10.0;
107 hso[0][2] = 3.0*cosd*(5.0-sqrt(5.0))/10.0;
108 hso[0][1] = sqrt(3.0)*(sqrt(5.0)-1.0)/10.0;
109 hso[4][0] = hso[0][4];
110 hso[4][4] = (2.0*sqrt(5.0)+1.0)/10.0;
111 hso[4][3] = sqrt(3.0)*cosd*(5.0-2.0*sqrt(5.0))/10.0;
112 hso[4][2] = sqrt(3.0)*cosd*(5.0-3.0*sqrt(5.0))/5.0;
113 hso[4][1] = 2.0/5.0;
114 hso[3][0] = -hso[0][3];
115 hso[3][4] = -hso[4][3];
116 hso[3][3] = -1.0/2.0;
117 hso[3][2] = 0.0;
118 hso[3][1] = sqrt(3.0)*cosd*(5.0-sqrt(5.0))/5.0;
119 hso[2][0] = -hso[0][2];
120 hso[2][4] = -hso[4][2];
121 hso[2][3] = 0.0;
122 hso[2][2] = -1.0/2.0;
123 hso[2][1] = -sqrt(3.0)*sqrt(5.0)*cosd/10.0;
124 hso[1][0] = hso[0][1];
125 hso[1][4] = hso[4][1];
126 hso[1][3] = -hso[3][1];
127 hso[1][2] = -hso[2][1];
128 hso[1][1] = (1.0-2.0*sqrt(5.0))/10.0;
129
130 // now rotate the first C5's by 2pi/3 degrees about the zx axis (sort of)
131 t1rep[3] = t1rep[1].transform(t1so);
132 t1rep[4] = t1rep[2].transform(t1so);
133
134 grep[3] = grep[1].transform(gso);
135 grep[4] = grep[2].transform(gso);
136
137 hrep[3] = hrep[1].transform(hso);
138 hrep[4] = hrep[2].transform(hso);
139
140 // rotate twice to get the first one aligned along the x axis
141 t1rep[3] = t1rep[3].transform(t1rep[1]).transform(t1rep[1]);
142 t1rep[4] = t1rep[4].transform(t1rep[1]).transform(t1rep[1]);
143
144 grep[3] = grep[3].transform(grep[1]).transform(grep[1]);
145 grep[4] = grep[4].transform(grep[1]).transform(grep[1]);
146
147 hrep[3] = hrep[3].transform(hrep[1]).transform(hrep[1]);
148 hrep[4] = hrep[4].transform(hrep[1]).transform(hrep[1]);
149
150 t2rep[3] = t1rep[4].operate(t1rep[4]);
151 t2rep[4] = t1rep[3].operate(t1rep[3]);
152
153 t2rep[13] = t1rep[2];
154 t2rep[14] = t1rep[1];
155
156 t2rep[15] = t1rep[3];
157 t2rep[16] = t1rep[4];
158
159 // and then rotate those by 2pi/5 about the z axis 4 times
160 for (i=5; i < 13; i++) {
161 t1rep[i] = t1rep[i-2].transform(t1rep[1]);
162 grep[i] = grep[i-2].transform(grep[1]);
163 hrep[i] = hrep[i-2].transform(hrep[1]);
164
165 t2rep[i] = t2rep[i-2].transform(t2rep[1]);
166 t2rep[i+12] = t2rep[i+10].transform(t2rep[1]);
167 }
168
169 //
170 // 12 C5^2's
171 //
172 // get these from operating on each of the C5's with itself
173 for (i=13; i < 25; i++) {
174 t1rep[i] = t1rep[i-12].operate(t1rep[i-12]);
175 grep[i] = grep[i-12].operate(grep[i-12]);
176 hrep[i] = hrep[i-12].operate(hrep[i-12]);
177 }
178
179 //
180 // 20 C3's
181 //
182 // first we have 2 C3's about the zx axis
183 t1rep[25] = t1so;
184 t1rep[26] = t1so.operate(t1so);
185
186 grep[25] = gso;
187 grep[26] = gso.operate(gso);
188
189 hrep[25] = hso;
190 hrep[26] = hso.operate(hso);
191
192 // and then rotate those by 2pi/5 about the z axis 4 times
193 for (i=27; i < 35; i++) {
194 t1rep[i] = t1rep[i-2].transform(t1rep[1]);
195 grep[i] = grep[i-2].transform(grep[1]);
196 hrep[i] = hrep[i-2].transform(hrep[1]);
197 }
198
199 // now rotate one of the above C3's by 2pi/3 about the zx axis
200 t1rep[35] = t1rep[27].transform(t1so);
201 t1rep[36] = t1rep[28].transform(t1so);
202
203 grep[35] = grep[27].transform(gso);
204 grep[36] = grep[28].transform(gso);
205
206 hrep[35] = hrep[27].transform(hso);
207 hrep[36] = hrep[28].transform(hso);
208
209 // and then rotate those by 2pi/5 about the z axis 4 times
210 for (i=37; i < 45; i++) {
211 t1rep[i] = t1rep[i-2].transform(t1rep[1]);
212 grep[i] = grep[i-2].transform(grep[1]);
213 hrep[i] = hrep[i-2].transform(hrep[1]);
214 }
215
216 t2rep[25] = t1rep[35];
217 t2rep[26] = t1rep[36];
218
219 for (i=27; i < 35; i++)
220 t2rep[i] = t2rep[i-2].transform(t2rep[1]);
221
222 t2rep[35] = t1rep[26];
223 t2rep[36] = t1rep[25];
224
225 for (i=37; i < 45; i++)
226 t2rep[i] = t2rep[i-2].transform(t2rep[1]);
227
228 //
229 // 15 C2's
230 //
231 // first we have a C2 about the y axis
232 t1rep[45][0][0] = -1.0;
233 t1rep[45][1][1] = 1.0;
234 t1rep[45][2][2] = -1.0;
235
236 t2rep[45] = t1rep[45];
237
238 grep[45][0][0] = -1.0;
239 grep[45][1][1] = 1.0;
240 grep[45][2][2] = -1.0;
241 grep[45][3][3] = 1.0;
242
243 hrep[45][0][0] = 1.0;
244 hrep[45][1][1] = 1.0;
245 hrep[45][2][2] = -1.0;
246 hrep[45][3][3] = -1.0;
247 hrep[45][4][4] = 1.0;
248
249 // and rotate that by 2pi/5 about the z axis 4 times
250 for (i=46; i < 50; i++) {
251 t1rep[i] = t1rep[i-1].transform(t1rep[1]);
252 t2rep[i] = t2rep[i-1].transform(t2rep[1]);
253 grep[i] = grep[i-1].transform(grep[1]);
254 hrep[i] = hrep[i-1].transform(hrep[1]);
255 }
256
257 // now take the C2 about the y axis and rotate it by 2pi/3 about the zx axis
258 t1rep[50] = t1rep[45].transform(t1so);
259 grep[50] = grep[45].transform(gso);
260 hrep[50] = hrep[45].transform(hso);
261
262 // align this c2 along the x axis
263 t1rep[50] = t1rep[50].transform(t1rep[2]).transform(t1rep[2]);
264 grep[50] = grep[50].transform(grep[2]).transform(grep[2]);
265 hrep[50] = hrep[50].transform(hrep[2]).transform(hrep[2]);
266
267 // and rotate that by 2pi/5 about the z axis 4 times
268 for (i=51; i < 55; i++) {
269 t1rep[i] = t1rep[i-1].transform(t1rep[1]);
270 grep[i] = grep[i-1].transform(grep[1]);
271 hrep[i] = hrep[i-1].transform(hrep[1]);
272 }
273
274 // finally, take a C2 about the y axis, and rotate it by 2pi/3 about the
275 // xz axis, and align it along the x axis
276 t1rep[55] = t1rep[45].transform(t1rep[35]).transform(t1rep[1]);
277 grep[55] = grep[45].transform(grep[35]).transform(grep[1]);
278 hrep[55] = hrep[45].transform(hrep[35]).transform(hrep[1]);
279
280 // and then rotate that by 2pi/5 about the z axis 4 times
281 for (i=56; i < 60; i++) {
282 t1rep[i] = t1rep[i-1].transform(t1rep[1]);
283 grep[i] = grep[i-1].transform(grep[1]);
284 hrep[i] = hrep[i-1].transform(hrep[1]);
285 }
286
287 t2rep[50] = t1rep[55];
288 t2rep[55] = t1rep[50];
289
290 for (i=51; i < 55; i++) {
291 t2rep[i] = t2rep[i-1].transform(t2rep[1]);
292 t2rep[i+5] = t2rep[i+4].transform(t2rep[1]);
293 }
294 }
295
296 void
i()297 CharacterTable::i()
298 {
299 int i;
300
301 IrreducibleRepresentation& ira = gamma_[0];
302 IrreducibleRepresentation& ir1 = gamma_[1];
303 IrreducibleRepresentation& ir2 = gamma_[2];
304 IrreducibleRepresentation& irg = gamma_[3];
305 IrreducibleRepresentation& irh = gamma_[4];
306
307 ira.init(g,1,"A");
308 ir1.init(g,3,"T1");
309 ir2.init(g,3,"T2");
310 irg.init(g,4,"G");
311 irh.init(g,5,"H");
312
313 // i_ops gives us all the symmetry operations we need
314 i_ops(ir1.rep, ir2.rep, irg.rep, irh.rep);
315
316 ir1.nrot_ = 1;
317 ir1.ntrans_ = 1;
318
319 for (i=0; i < g; i++) {
320 ira.rep[i][0][0] = 1.0;
321 symop[i] = ir1.rep[i];
322 }
323 }
324
325
ih()326 void CharacterTable::ih()
327 {
328 int i;
329
330 IrreducibleRepresentation& irag = gamma_[0];
331 IrreducibleRepresentation& ir1g = gamma_[1];
332 IrreducibleRepresentation& ir2g = gamma_[2];
333 IrreducibleRepresentation& irgg = gamma_[3];
334 IrreducibleRepresentation& irhg = gamma_[4];
335
336 IrreducibleRepresentation& irau = gamma_[5];
337 IrreducibleRepresentation& ir1u = gamma_[6];
338 IrreducibleRepresentation& ir2u = gamma_[7];
339 IrreducibleRepresentation& irgu = gamma_[8];
340 IrreducibleRepresentation& irhu = gamma_[9];
341
342 irag.init(g,1,"Ag");
343 ir1g.init(g,3,"T1g");
344 ir2g.init(g,3,"T2g");
345 irgg.init(g,4,"Gg");
346 irhg.init(g,5,"Hg");
347
348 irau.init(g,1,"Au");
349 ir1u.init(g,3,"T1u");
350 ir2u.init(g,3,"T2u");
351 irgu.init(g,4,"Gu");
352 irhu.init(g,5,"Hu");
353
354 // i_ops gives us all the symmetry operations we need
355 i_ops(ir1g.rep, ir2g.rep, irgg.rep, irhg.rep);
356
357 ir1g.nrot_ = 1;
358 ir1u.ntrans_ = 1;
359
360 SymRep ti(3), gi(4), hi(5);
361 ti.i();
362 gi.i();
363 hi.i();
364
365 for (i=0; i < g/2; i++) {
366 irag.rep[i][0][0] = 1.0;
367 irau.rep[i][0][0] = 1.0;
368
369 irag.rep[i+60][0][0] = 1.0;
370 irau.rep[i+60][0][0] = -1.0;
371
372 ir1g.rep[i+60] = ir1g.rep[i];
373 ir2g.rep[i+60] = ir2g.rep[i];
374 irgg.rep[i+60] = irgg.rep[i];
375 irhg.rep[i+60] = irhg.rep[i];
376
377 ir1u.rep[i] = ir1g.rep[i];
378 ir2u.rep[i] = ir2g.rep[i];
379 irgu.rep[i] = irgg.rep[i];
380 irhu.rep[i] = irhg.rep[i];
381
382 ir1u.rep[i+60] = ir1g.rep[i].operate(ti);
383 ir2u.rep[i+60] = ir2g.rep[i].operate(ti);
384 irgu.rep[i+60] = irgg.rep[i].operate(gi);
385 irhu.rep[i+60] = irhg.rep[i].operate(hi);
386
387 symop[i] = ir1u.rep[i];
388 symop[i+60] = ir1u.rep[i+60];
389 }
390 }
391
392 /////////////////////////////////////////////////////////////////////////////
393
394 // Local Variables:
395 // mode: c++
396 // c-file-style: "ETS"
397 // End:
398