source: ThirdParty/mpqc_open/src/lib/math/symmetry/maketab.cc@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 25.1 KB
Line 
1//
2// maketab.cc
3//
4// Modifications are
5// Copyright (C) 1996 Limit Point Systems, Inc.
6//
7// Author: Edward Seidl <seidl@janed.com>
8// Maintainer: LPS
9//
10// This file is part of the SC Toolkit.
11//
12// The SC Toolkit is free software; you can redistribute it and/or modify
13// it under the terms of the GNU Library General Public License as published by
14// the Free Software Foundation; either version 2, or (at your option)
15// any later version.
16//
17// The SC Toolkit is distributed in the hope that it will be useful,
18// but WITHOUT ANY WARRANTY; without even the implied warranty of
19// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20// GNU Library General Public License for more details.
21//
22// You should have received a copy of the GNU Library General Public License
23// along with the SC Toolkit; see the file COPYING.LIB. If not, write to
24// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
25//
26// The U.S. Government is granted a limited license as per AL 91-7.
27//
28
29/* maketab.cc
30 *
31 * THIS SOFTWARE FITS THE DESCRIPTION IN THE U.S. COPYRIGHT ACT OF A
32 * "UNITED STATES GOVERNMENT WORK". IT WAS WRITTEN AS A PART OF THE
33 * AUTHOR'S OFFICIAL DUTIES AS A GOVERNMENT EMPLOYEE. THIS MEANS IT
34 * CANNOT BE COPYRIGHTED. THIS SOFTWARE IS FREELY AVAILABLE TO THE
35 * PUBLIC FOR USE WITHOUT A COPYRIGHT NOTICE, AND THERE ARE NO
36 * RESTRICTIONS ON ITS USE, NOW OR SUBSEQUENTLY.
37 *
38 * Author:
39 * E. T. Seidl
40 * Bldg. 12A, Rm. 2033
41 * Computer Systems Laboratory
42 * Division of Computer Research and Technology
43 * National Institutes of Health
44 * Bethesda, Maryland 20892
45 * Internet: seidl@alw.nih.gov
46 * June, 1993
47 */
48
49#include <util/misc/math.h>
50#include <stdio.h>
51#include <string.h>
52
53#include <math/symmetry/pointgrp.h>
54#include <util/misc/formio.h>
55
56using namespace std;
57using namespace sc;
58
59/*
60 * This function will generate a character table for the point group.
61 * This character table is in the order that symmetry operations are
62 * generated, not in Cotton order. If this is a problem, tough.
63 * Also generate the transformation matrices.
64 */
65
66int CharacterTable::make_table()
67{
68 int i,j,ei,gi;
69 char label[4];
70
71 if (!g) return 0;
72
73 gamma_ = new IrreducibleRepresentation[nirrep_];
74
75 symop = new SymmetryOperation[g];
76 SymmetryOperation so;
77
78 _inv = new int[g];
79
80 // this array forms a reducible representation for rotations about x,y,z
81 double *rot = new double[g];
82 memset(rot,0,sizeof(double)*g);
83
84 // this array forms a reducible representation for translations along x,y,z
85 double *trans = new double[g];
86 memset(trans,0,sizeof(double)*g);
87
88 // the angle to rotate about the principal axis
89 double theta = (nt) ? 2.0*M_PI/nt : 2.0*M_PI;
90
91 switch (pg) {
92
93 case C1:
94 // no symmetry case
95 gamma_[0].init(1,1,"A");
96 gamma_[0].nrot_ = 3;
97 gamma_[0].ntrans_ = 3;
98 gamma_[0].rep[0][0][0] = 1.0;
99
100 symop[0].E();
101
102 break;
103
104 case CI:
105 // equivalent to S2 about the z axis
106 gamma_[0].init(2,1,"Ag");
107 gamma_[0].rep[0][0][0] = 1.0;
108 gamma_[0].rep[1][0][0] = 1.0;
109 gamma_[0].nrot_=3;
110
111 gamma_[1].init(2,1,"Au");
112 gamma_[1].rep[0][0][0] = 1.0;
113 gamma_[1].rep[1][0][0] = -1.0;
114 gamma_[1].ntrans_=3;
115
116 symop[0].E();
117 symop[1].i();
118
119 break;
120
121 case CS: // reflection through the xy plane
122 gamma_[0].init(2,1,"A'","Ap");
123 gamma_[0].rep[0][0][0] = 1.0;
124 gamma_[0].rep[1][0][0] = 1.0;
125 gamma_[0].nrot_=1;
126 gamma_[0].ntrans_=2;
127
128 gamma_[1].init(2,1,"A\"","App");
129 gamma_[1].rep[0][0][0] = 1.0;
130 gamma_[1].rep[1][0][0] = -1.0;
131 gamma_[1].nrot_=2;
132 gamma_[1].ntrans_=1;
133
134 symop[0].E();
135 symop[1].sigma_h();
136
137 break;
138
139 case CN:
140 // clockwise rotation about z axis by theta*i radians
141 //
142 // for odd n, the irreps are A and E1...E(nir-1)
143 // for even n, the irreps are A, B, and E1...E(nir-2)
144 //
145 gamma_[0].init(g,1,"A");
146 for (gi=0; gi < g; gi++)
147 gamma_[0].rep[gi][0][0] = 1.0;
148
149 i=1;
150
151 if (!(nt%2)) {
152 gamma_[1].init(g,1,"B");
153 for (gi=0; gi < g; gi++)
154 gamma_[1].rep[gi][0][0] = (gi%2) ? -1.0 : 1.0;
155
156 i++;
157 }
158
159 ei=1;
160 for (; i < nirrep_; i++, ei++) {
161 IrreducibleRepresentation& ir = gamma_[i];
162
163 if (nt==3 || nt==4)
164 sprintf(label,"E");
165 else
166 sprintf(label,"E%d",ei);
167
168 ir.init(g,2,label);
169 ir.complex_=1;
170
171 // identity
172 ir.rep[0].E();
173
174 // Cn
175 ir.rep[1].rotation(ei*theta);
176
177 // the other n-1 Cn's
178 for (j=2; j < g; j++)
179 ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
180 }
181
182 // identity
183 symop[0].E();
184
185 // Cn
186 symop[1].rotation(theta);
187
188 // the other n-2 Cn's
189 for (i=2; i < nt; i++)
190 symop[i] = symop[i-1].operate(symop[1]);
191
192 for (i=0; i < nt ; i++)
193 rot[i] = trans[i] = symop[i].trace();
194
195 break;
196
197 case CNV:
198 // clockwise rotation about z axis by theta*i radians, then
199 // reflect through the xz plane
200 //
201 // for odd n, the irreps are A1, A2, and E1...E(nir-2)
202 // for even n, the irreps are A1, A2, B1, B2, and E1...E(nir-4)
203 //
204
205 gamma_[0].init(g,1,"A1");
206 gamma_[1].init(g,1,"A2");
207
208 for (gi=0; gi < nt; gi++) {
209 // Cn's
210 gamma_[0].rep[gi][0][0] = 1.0;
211 gamma_[1].rep[gi][0][0] = 1.0;
212
213 // sigma's
214 gamma_[0].rep[gi+nt][0][0] = 1.0;
215 gamma_[1].rep[gi+nt][0][0] = -1.0;
216 }
217
218 if (!(nt%2)) {
219 gamma_[2].init(g,1,"B1");
220 gamma_[3].init(g,1,"B2");
221
222 for (gi=0; gi < nt ; gi++) {
223 double ci = (gi%2) ? -1.0 : 1.0;
224
225 // Cn's
226 gamma_[2].rep[gi][0][0] = ci;
227 gamma_[3].rep[gi][0][0] = ci;
228
229 // sigma's
230 gamma_[2].rep[gi+nt][0][0] = ci;
231 gamma_[3].rep[gi+nt][0][0] = -ci;
232 }
233 }
234
235 ei=1;
236 for (i = (nt%2) ? 2 : 4; i < nirrep_; i++, ei++) {
237 IrreducibleRepresentation& ir = gamma_[i];
238
239 char lab[4];
240 if (nt==3 || nt==4)
241 sprintf(lab,"E");
242 else
243 sprintf(lab,"E%d",ei);
244
245 ir.init(g,2,lab);
246
247 // identity
248 ir.rep[0].E();
249
250 // Cn
251 ir.rep[1].rotation(ei*theta);
252
253 // the other n-2 Cn's
254 for (j=2; j < nt; j++)
255 ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
256
257 // sigma xz
258 ir.rep[nt].sigma_xz();
259
260 SymRep sr(2);
261 sr.rotation(ei*theta/2.0);
262
263 // the other n-1 sigma's
264 for (j=nt+1; j < g; j++)
265 ir.rep[j] = ir.rep[j-1].transform(sr);
266 }
267
268 // identity
269 symop[0].E();
270
271 // Cn
272 symop[1].rotation(theta);
273
274 // the other n-2 Cn's
275 for (i=2; i < nt; i++)
276 symop[i] = symop[i-1].operate(symop[1]);
277
278 // sigma xz
279 symop[nt].sigma_xz();
280
281 so.rotation(theta/2.0);
282
283 // the other n-1 sigma's
284 for (j=nt+1; j < g; j++)
285 symop[j] = symop[j-1].transform(so);
286
287 for (i=0; i < nt ; i++) {
288 rot[i] = trans[i] = symop[i].trace();
289
290 rot[i+nt] = -symop[i+nt].trace();
291 trans[i+nt] = symop[i+nt].trace();
292 }
293
294 break;
295
296 case CNH:
297 // lockwise rotation about z axis by theta*i radians,
298 // as well as rotation-reflection about same axis
299
300 //
301 // for odd n, the irreps are A', A", E1'...E(nir/2-1)', E1"...E(nir/2-1)''
302 // for even n, the irreps are Ag, Bg, Au, Bu,
303 // E1g...E(nir/2-1)g, E1u...E(nir/2-1)u
304 //
305 gamma_[0].init(g,1, (nt%2) ? "A'" : "Ag", (nt%2) ? "Ap" : 0);
306 gamma_[nirrep_/2].init(g,1, (nt%2) ? "A\"" : "Au", (nt%2) ? "Ap" : 0);
307
308 for (gi=0; gi < nt; gi++) {
309 // Cn's
310 gamma_[0].rep[gi][0][0] = 1.0;
311 gamma_[nirrep_/2].rep[gi][0][0] = 1.0;
312
313 // Sn's
314 gamma_[0].rep[gi+nt][0][0] = 1.0;
315 gamma_[nirrep_/2].rep[gi+nt][0][0] = -1.0;
316 }
317
318 if (!(nt%2)) {
319 gamma_[1].init(g,1,"Bg");
320 gamma_[1+nirrep_/2].init(g,1,"Bu");
321
322 for (gi=0; gi < nt; gi++) {
323 double ci = (gi%2) ? -1.0 : 1.0;
324
325 // Cn's
326 gamma_[1].rep[gi][0][0] = ci;
327 gamma_[1+nirrep_/2].rep[gi][0][0] = ci;
328
329 // Sn's
330 gamma_[1].rep[gi+nt][0][0] = ci;
331 gamma_[1+nirrep_/2].rep[gi+nt][0][0] = -ci;
332 }
333 }
334
335 ei=1;
336 for (i = (nt%2) ? 1 : 2; i < nirrep_/2 ; i++, ei++) {
337 IrreducibleRepresentation& ir1 = gamma_[i];
338 IrreducibleRepresentation& ir2 = gamma_[i+nirrep_/2];
339
340 if (nt==3 || nt==4)
341 sprintf(label,(nt%2) ? "E'" : "Eg");
342 else
343 sprintf(label,"E%d%s", ei, (nt%2) ? "'" : "g");
344
345 ir1.init(g,2,label);
346
347 if (nt==3 || nt==4)
348 sprintf(label,(nt%2) ? "E\"" : "Eu");
349 else
350 sprintf(label,"E%d%s", ei, (nt%2) ? "\"" : "u");
351
352 ir2.init(g,2,label);
353
354 ir1.complex_=1;
355 ir2.complex_=1;
356
357 // identity
358 ir1.rep[0].E();
359 ir2.rep[0].E();
360
361 // Cn
362 ir1.rep[1].rotation(ei*theta);
363 ir2.rep[1].rotation(ei*theta);
364
365 for (j=2; j < nt; j++) {
366 ir1.rep[j] = ir1.rep[j-1].operate(ir1.rep[1]);
367 ir2.rep[j] = ir2.rep[j-1].operate(ir2.rep[1]);
368 }
369
370 // Sn's
371 SymRep sr(2);
372 sr.i();
373
374 for (j=nt; j < g; j++) {
375 ir1.rep[j] = ir1.rep[j-nt];
376 ir2.rep[j] = ir2.rep[j-nt].operate(sr);
377 }
378 }
379
380 // identity
381 symop[0].E();
382
383 // Cn
384 symop[1].rotation(theta);
385
386 // the other n-2 Cn's
387 for (i=2; i < nt; i++)
388 symop[i] = symop[i-1].operate(symop[1]);
389
390 // Sn's, for odd nt, operate on Cn's with sigma_h, for even nt,
391 // operate Cn's with i
392 if (nt%2)
393 so.sigma_h();
394 else
395 so.i();
396
397 for (i=0; i < nt ; i++) {
398 symop[i+nt] = symop[i].operate(so);
399
400 rot[i] = trans[i] = symop[i].trace();
401 trans[i+nt] = symop[i+nt].trace();
402 rot[i+nt] = -trans[i+nt];
403 }
404
405 break;
406
407 case SN:
408 // clockwise rotation-reflection by theta*i radians about z axis
409 //
410 // for odd n/2, the irreps are Ag, Au, E1g...E(nir/2-1)g,E1u...E(nir/2-1)u
411 // for even n/2, the irreps are A, B, E1...E(nir-2)
412 //
413 if ((nt/2)%2) {
414 gamma_[0].init(g, 1, "Ag");
415 gamma_[nirrep_/2].init(g, 1, "Au");
416
417 for (gi=0; gi < nt/2; gi++) {
418 gamma_[0].rep[gi][0][0] = 1.0;
419 gamma_[nirrep_/2].rep[gi][0][0] = 1.0;
420
421 gamma_[0].rep[gi+nt/2][0][0] = 1.0;
422 gamma_[nirrep_/2].rep[gi+nt/2][0][0] = -1.0;
423 }
424
425 ei=1;
426 for (i=1; i < nirrep_/2 ; i++, ei++) {
427 IrreducibleRepresentation& ir1 = gamma_[i];
428 IrreducibleRepresentation& ir2 = gamma_[i+nirrep_/2];
429
430 if (nt==6)
431 sprintf(label,"Eg");
432 else
433 sprintf(label,"E%dg",ei);
434
435 ir1.init(g,2,label);
436 ir1.complex_=1;
437
438 if (nt==6)
439 sprintf(label,"Eu");
440 else
441 sprintf(label,"E%du", ei);
442
443 ir2.init(g,2,label);
444 ir2.complex_=1;
445
446 // identity
447 ir1.rep[0].E();
448 ir2.rep[0].E();
449
450 // C(n/2)
451 ir1.rep[1].rotation(ei*theta*2.0);
452 ir2.rep[1].rotation(ei*theta*2.0);
453
454 for (j=2; j < nt/2; j++) {
455 ir1.rep[j] = ir1.rep[j-1].operate(ir1.rep[1]);
456 ir2.rep[j] = ir2.rep[j-1].operate(ir2.rep[1]);
457 }
458
459 SymRep sr(2);
460 sr.i();
461
462 // Sn
463 for (j=nt/2; j < nt; j++) {
464 ir1.rep[j] = ir1.rep[j-nt/2];
465 ir2.rep[j] = ir2.rep[j-nt/2].operate(sr);
466 }
467 }
468
469 // identity
470 symop[0].E();
471
472 // Cn
473 symop[1].rotation(2.0*theta);
474
475 for (i=2; i < nt/2 ; i++)
476 symop[i] = symop[i-1].operate(symop[1]);
477
478 so.i();
479
480 // Sn
481 for (i=nt/2; i < nt; i++)
482 symop[i] = symop[i-nt/2].operate(so);
483
484 for (i=0; i < nt/2 ; i++) {
485 rot[i] = trans[i] = symop[i].trace();
486
487 trans[i+nt/2] = symop[i+nt/2].trace();
488 rot[i+nt/2] = -trans[i+nt/2];
489 }
490
491 } else {
492 gamma_[0].init(g, 1, "A");
493 gamma_[1].init(g, 1, "B");
494
495 for (gi=0; gi < nt; gi++) {
496 gamma_[0].rep[gi][0][0] = 1.0;
497 gamma_[1].rep[gi][0][0] = (gi%2) ? -1.0 : 1.0;
498 }
499
500 ei=1;
501 for (i=2; i < nirrep_; i++, ei++) {
502 IrreducibleRepresentation& ir = gamma_[i];
503
504 if (nt==4)
505 sprintf(label,"E");
506 else
507 sprintf(label,"E%d",ei);
508
509 ir.init(g,2,label);
510 ir.complex_ = 1;
511
512 // identity
513 ir.rep[0].E();
514
515 // Sn
516 ir.rep[1].rotation(ei*theta);
517
518 for (j=2; j < nt; j++)
519 ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
520 }
521
522 // identity
523 symop[0].E();
524
525 // Sn
526 symop[1].rotation(theta);
527 symop[1][2][2] = -1.0;
528
529 for (i=2; i < nt ; i++)
530 symop[i] = symop[i-1].operate(symop[1]);
531
532 for (i=0; i < nt ; i++) {
533 trans[i] = symop[i].trace();
534 rot[i] = (i%2) ? -trans[i] : trans[i];
535 }
536 }
537
538 break;
539
540 case DN:
541 // clockwise rotation about z axis, followed by C2 about x axis
542
543 // D2 is a special case
544 if (nt==2) {
545 gamma_[0].init(g,1,"A");
546 gamma_[1].init(g,1,"B1");
547 gamma_[2].init(g,1,"B2");
548 gamma_[3].init(g,1,"B3");
549
550 for (i=0; i < g; i++) {
551 gamma_[0].rep[i][0][0] = 1.0;
552 gamma_[1].rep[i][0][0] = (i < 2) ? 1.0 : -1.0;
553 gamma_[2].rep[i][0][0] = (i % 2) ? -1.0 : 1.0;
554 gamma_[3].rep[i][0][0] = (i < 2) ?
555 ((i % 2) ? -1.0 : 1.0) : ((i%2) ? 1.0 : -1.0);
556 }
557 } else {
558 // Dn is isomorphic with Cnv
559 //
560 // for odd n, the irreps are A1, A2, and E1...E(nir-2)
561 // for even n, the irreps are A1, A2, B1, B2, and E1...E(nir-4)
562 //
563 gamma_[0].init(g,1,"A1");
564 gamma_[1].init(g,1,"A2");
565
566 for (gi=0; gi < nt; gi++) {
567 // Cn's
568 gamma_[0].rep[gi][0][0] = 1.0;
569 gamma_[1].rep[gi][0][0] = 1.0;
570
571 // C2's
572 gamma_[0].rep[gi+nt][0][0] = 1.0;
573 gamma_[1].rep[gi+nt][0][0] = -1.0;
574 }
575
576 i=2;
577
578 if (!(nt%2)) {
579 gamma_[2].init(g,1,"B1");
580 gamma_[3].init(g,1,"B2");
581
582 for (gi=0; gi < nt ; gi++) {
583 double ci = (gi%2) ? -1.0 : 1.0;
584
585 // Cn's
586 gamma_[2].rep[gi][0][0] = ci;
587 gamma_[3].rep[gi][0][0] = ci;
588
589 // sigma's
590 gamma_[2].rep[gi+nt][0][0] = ci;
591 gamma_[3].rep[gi+nt][0][0] = -ci;
592 }
593
594 i = 4;
595 }
596
597 ei=1;
598 for (; i < nirrep_; i++, ei++) {
599 IrreducibleRepresentation& ir = gamma_[i];
600
601 char lab[4];
602 if (nt==3 || nt==4)
603 sprintf(lab,"E");
604 else
605 sprintf(lab,"E%d",ei);
606
607 ir.init(g,2,lab);
608
609 // identity
610 ir.rep[0].E();
611
612 // Cn
613 ir.rep[1].rotation(ei*theta);
614
615 for (j=2; j < nt; j++)
616 ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
617
618 // C2(x)
619 ir.rep[nt].c2_y();
620
621 SymRep sr(2);
622 sr.rotation(ei*theta/2.0);
623
624 for (j=nt+1; j < 2*nt; j++)
625 ir.rep[j] = ir.rep[j-1].transform(sr);
626 }
627 }
628
629 // identity
630 symop[0].E();
631
632 // Cn
633 symop[1].rotation(theta);
634
635 for (i=2; i < nt; i++)
636 symop[i] = symop[i-1].operate(symop[1]);
637
638 // C2(x)
639 symop[nt].c2_y();
640
641 so.rotation(theta/2.0);
642
643 for (i=nt+1; i < 2*nt; i++)
644 symop[i] = symop[i-1].transform(so);
645
646 for (i=0; i < 2*nt ; i++)
647 rot[i] = trans[i] = symop[i].trace();
648
649 break;
650
651 case DND:
652 // rotation reflection about z axis by theta/2 radians, followed
653 // by c2 about x axis, then reflection through yz plane
654 //
655 // for odd n, the irreps are A1g, A2g, A1u, A2u, E1g...E(nir/2-2)g,
656 // E1u...E(nir/2-2)u
657 // for even n, the irreps are A1, A2, B1, B2, E1...E(nir-4)
658 //
659
660 if (nt%2) {
661 gamma_[0].init(g,1,"A1g");
662 gamma_[1].init(g,1,"A2g");
663
664 for (gi=0; gi < g; gi++) {
665 gamma_[0].rep[gi][0][0] = 1.0;
666 gamma_[1].rep[gi][0][0] = (gi/nt==0 || gi/nt==2) ? 1.0 : -1.0;
667 }
668
669 i=nirrep_/2;
670 j=i+1;
671 gamma_[i].init(g,1,"A1u");
672 gamma_[j].init(g,1,"A2u");
673
674 for (gi=0; gi < g/2; gi++) {
675 gamma_[i].rep[gi][0][0] = gamma_[0].rep[gi][0][0];
676 gamma_[j].rep[gi][0][0] = gamma_[1].rep[gi][0][0];
677
678 gamma_[i].rep[gi+g/2][0][0] = -gamma_[0].rep[gi][0][0];
679 gamma_[j].rep[gi+g/2][0][0] = -gamma_[1].rep[gi][0][0];
680 }
681
682 ei=1;
683
684 for (i=2; i < nirrep_/2 ; i++, ei++) {
685 IrreducibleRepresentation& ir1 = gamma_[i];
686 IrreducibleRepresentation& ir2 = gamma_[i+nirrep_/2];
687
688 if (nt==3) {
689 ir1.init(g,2,"Eg");
690 ir2.init(g,2,"Eu");
691 } else {
692 sprintf(label,"E%dg",ei);
693 ir1.init(g,2,label);
694 sprintf(label,"E%du",ei);
695 ir2.init(g,2,label);
696 }
697
698 // identity
699 ir1.rep[0].E();
700
701 // Cn
702 ir1.rep[1].rotation(ei*theta);
703
704 for (j=2; j < nt; j++)
705 ir1.rep[j] = ir1.rep[j-1].operate(ir1.rep[1]);
706
707 // C2(x)
708 ir1.rep[nt].c2_y();
709
710 for (j=nt+1; j < 2*nt; j++)
711 ir1.rep[j] = ir1.rep[j-1].transform(ir1.rep[1]);
712
713 for (j=0; j < 2*nt; j++)
714 ir2.rep[j] = ir1.rep[j];
715
716 // Sn and sigma d
717 SymRep sr(2);
718 sr.i();
719
720 for (j=2*nt; j < g; j++) {
721 ir1.rep[j] = ir1.rep[j-2*nt];
722 ir2.rep[j] = ir2.rep[j-2*nt].operate(sr);
723 }
724 }
725
726 // identity
727 symop[0].E();
728
729 // Cn
730 symop[1].rotation(theta);
731
732 for (i=2; i < nt; i++)
733 symop[i] = symop[i-1].operate(symop[1]);
734
735 // C2(x)
736 symop[nt].c2_y();
737
738 for (i=nt+1; i < 2*nt; i++)
739 symop[i] = symop[i-1].transform(symop[1]);
740
741 // i + n-1 S2n + n sigma
742 so.i();
743 for (i=2*nt; i < g; i++)
744 symop[i] = symop[i-2*nt].operate(so);
745
746 for (i=0; i < g; i++) {
747 trans[i] = symop[i].trace();
748 rot[i] = (i < g/2) ? trans[i] : -trans[i];
749 }
750
751 } else { // even nt
752
753 gamma_[0].init(g,1,"A1");
754 gamma_[1].init(g,1,"A2");
755 gamma_[2].init(g,1,"B1");
756 gamma_[3].init(g,1,"B2");
757
758 for (gi=0; gi < 2*nt; gi++) {
759 // Sn
760 gamma_[0].rep[gi][0][0] = 1.0;
761 gamma_[1].rep[gi][0][0] = 1.0;
762 gamma_[2].rep[gi][0][0] = (gi%2) ? -1.0 : 1.0;
763 gamma_[3].rep[gi][0][0] = (gi%2) ? -1.0 : 1.0;
764
765 // n C2's and n sigma's
766 gamma_[0].rep[gi+2*nt][0][0] = 1.0;
767 gamma_[1].rep[gi+2*nt][0][0] = -1.0;
768 gamma_[2].rep[gi+2*nt][0][0] = (gi%2) ? -1.0 : 1.0;
769 gamma_[3].rep[gi+2*nt][0][0] = (gi%2) ? 1.0 : -1.0;
770 }
771
772 ei=1;
773 for (i=4; i < nirrep_; i++, ei++) {
774 IrreducibleRepresentation& ir = gamma_[i];
775
776 if (nt==2)
777 sprintf(label,"E");
778 else
779 sprintf(label,"E%d",ei);
780
781 ir.init(g,2,label);
782
783 // identity
784 ir.rep[0].E();
785
786 // S2n
787 ir.rep[1].rotation(ei*theta/2.0);
788
789 for (j=2; j < 2*nt; j++)
790 ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
791
792 // C2(x) + sigma_d
793 ir.rep[2*nt].c2_y();
794
795 for (j=2*nt+1; j < g; j++)
796 ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
797 }
798
799 // identity
800 symop[0].E();
801
802 // Sn's
803 symop[1].rotation(theta/2.0);
804 symop[1][2][2] = -1.0;
805
806 for (i=2; i < 2*nt; i++)
807 symop[i] = symop[i-1].operate(symop[1]);
808
809 // C2(x)
810 symop[2*nt].c2_y();
811
812 for (i=2*nt+1; i < g; i++)
813 symop[i] = symop[i-1].operate(symop[1]);
814
815 for (i=0; i < g; i++) {
816 trans[i] = symop[i].trace();
817 rot[i] = (i%2) ? -trans[i] : trans[i];
818 }
819 }
820
821 break;
822
823 case DNH:
824 // clockwise rotation and rotation-reflection about z axis,
825 // followed by c2 about x axis and then reflection
826 // through xz
827
828 i=nirrep_/2; j=i+1;
829
830 if (nt%2) {
831 gamma_[0].init(g,1,"A1'");
832 gamma_[1].init(g,1,"A2'");
833 gamma_[i].init(g,1,"A1\"");
834 gamma_[j].init(g,1,"A2\"");
835 } else {
836 if (nt==2) {
837 gamma_[0].init(g,1,"Ag");
838 gamma_[1].init(g,1,"B1g");
839 gamma_[4].init(g,1,"Au");
840 gamma_[5].init(g,1,"B1u");
841 } else {
842 gamma_[0].init(g,1,"A1g");
843 gamma_[1].init(g,1,"A2g");
844 gamma_[i].init(g,1,"A1u");
845 gamma_[j].init(g,1,"A2u");
846 }
847 }
848
849 for (gi=0; gi < nt; gi++) {
850 // E + n-1 Cn's
851 gamma_[0].rep[gi][0][0] = gamma_[1].rep[gi][0][0] =
852 gamma_[i].rep[gi][0][0] = gamma_[j].rep[gi][0][0] = 1.0;
853
854 // n C2's
855 gamma_[0].rep[gi+nt][0][0] = gamma_[i].rep[gi+nt][0][0] = 1.0;
856 gamma_[1].rep[gi+nt][0][0] = gamma_[j].rep[gi+nt][0][0] = -1.0;
857
858 // i + n-1 S2n's
859 gamma_[0].rep[gi+2*nt][0][0] = gamma_[1].rep[gi+2*nt][0][0] = 1.0;
860 gamma_[i].rep[gi+2*nt][0][0] = gamma_[j].rep[gi+2*nt][0][0] = -1.0;
861
862 // n sigma's
863 gamma_[0].rep[gi+3*nt][0][0] = gamma_[j].rep[gi+3*nt][0][0] = 1.0;
864 gamma_[i].rep[gi+3*nt][0][0] = gamma_[1].rep[gi+3*nt][0][0] = -1.0;
865 }
866
867 if (!(nt%2)) {
868 if (nt==2) {
869 gamma_[2].init(g,1,"B2g");
870 gamma_[3].init(g,1,"B3g");
871 gamma_[6].init(g,1,"B2u");
872 gamma_[7].init(g,1,"B3u");
873 } else {
874 gamma_[2].init(g,1,"B1g");
875 gamma_[3].init(g,1,"B2g");
876 gamma_[i+2].init(g,1,"B1u");
877 gamma_[j+2].init(g,1,"B2u");
878 }
879
880 for (gi=0; gi < nt; gi++) {
881 // E + n-1 Cn's
882 gamma_[2].rep[gi][0][0] = gamma_[3].rep[gi][0][0] =
883 gamma_[i+2].rep[gi][0][0] = gamma_[j+2].rep[gi][0][0] =
884 (gi%2) ? -1.0 : 1.0;
885
886 // n C2's
887 gamma_[2].rep[gi+nt][0][0] = gamma_[i+2].rep[gi+nt][0][0] =
888 (gi%2) ? -1.0 : 1.0;
889 gamma_[3].rep[gi+nt][0][0] = gamma_[j+2].rep[gi+nt][0][0] =
890 (gi%2) ? 1.0 : -1.0;
891
892 // i + n-1 S2n's
893 gamma_[2].rep[gi+2*nt][0][0] = gamma_[3].rep[gi+2*nt][0][0] =
894 (gi%2) ? -1.0 : 1.0;
895 gamma_[i+2].rep[gi+2*nt][0][0] = gamma_[j+2].rep[gi+2*nt][0][0] =
896 (gi%2) ? 1.0 : -1.0;
897
898 // n sigma's
899 gamma_[2].rep[gi+3*nt][0][0] = gamma_[j+2].rep[gi+3*nt][0][0] =
900 (gi%2) ? -1.0 : 1.0;
901 gamma_[i+2].rep[gi+3*nt][0][0] = gamma_[3].rep[gi+3*nt][0][0] =
902 (gi%2) ? 1.0 : -1.0;
903 }
904 }
905
906 ei=1;
907 for (i = (nt%2) ? 2 : 4; i < nirrep_/2 ; i++, ei++) {
908 IrreducibleRepresentation& ir1 = gamma_[i];
909 IrreducibleRepresentation& ir2 = gamma_[i+nirrep_/2];
910
911 if (nt==3) {
912 ir1.init(g,2,"E'");
913 ir2.init(g,2,"E\"");
914 } else if (nt==4) {
915 ir1.init(g,2,"Eg");
916 ir2.init(g,2,"Eu");
917 } else {
918 sprintf(label,"E%d%s", ei, (nt%2) ? "'" : "g");
919 ir1.init(g,2,label);
920
921 sprintf(label,"E%d%s", ei, (nt%2) ? "\"" : "u");
922 ir2.init(g,2,label);
923 }
924
925 // identity
926 ir1.rep[0].E();
927
928 // n-1 Cn's
929 ir1.rep[1].rotation(ei*theta);
930
931 for (j=2; j < nt; j++)
932 ir1.rep[j] = ir1.rep[j-1].operate(ir1.rep[1]);
933
934 // n C2's
935 ir1.rep[nt].c2_y();
936
937 SymRep sr(2);
938 sr.rotation(ei*theta/2.0);
939
940 for (j=nt+1; j < 2*nt; j++)
941 ir1.rep[j] = ir1.rep[j-1].transform(sr);
942
943 sr.i();
944 for (j=0; j < 2*nt; j++) {
945 ir1.rep[j+2*nt] = ir1.rep[j];
946 ir2.rep[j] = ir1.rep[j];
947 ir2.rep[j+2*nt] = ir1.rep[j].operate(sr);
948 }
949 }
950
951 // identity
952 symop[0].E();
953
954 // n-1 Cn's
955 symop[1].rotation(theta);
956
957 for (i=2; i < nt; i++)
958 symop[i] = symop[i-1].operate(symop[1]);
959
960 // n C2's
961 symop[nt].c2_y();
962
963 so.rotation(theta/2.0);
964 for (i=nt+1; i < 2*nt; i++)
965 symop[i] = symop[i-1].transform(so);
966
967 if (nt%2)
968 so.sigma_h();
969 else
970 so.i();
971
972 for (i=2*nt; i < g; i++)
973 symop[i] = symop[i-2*nt].operate(so);
974
975 for (i=0,j=2*nt; i < 2*nt ; i++,j++) {
976 rot[i] = trans[i] = symop[i].trace();
977 trans[j] = symop[j].trace();
978 rot[j] = -trans[j];
979 }
980
981 break;
982
983 case T:
984 t();
985 break;
986
987 case TH:
988 th();
989 break;
990
991 case TD:
992 td();
993 break;
994
995 case O:
996 o();
997 break;
998
999 case OH:
1000 oh();
1001 break;
1002
1003 case I:
1004 this->i();
1005 break;
1006
1007 case IH:
1008 ih();
1009 break;
1010
1011 default:
1012 return -1;
1013
1014 }
1015
1016 /* ok, we have the reducible representation of the rotations and
1017 * translations, now let's use projection operators to find out how many
1018 * rotations and translations there are in each irrep
1019 */
1020
1021 if (pg != C1 && pg != CI && pg != CS && pg != T && pg != TD && pg != TH &&
1022 pg != O && pg != OH && pg != I && pg != IH) {
1023
1024 for (i=0; i < nirrep_; i++) {
1025 double nr=0; double nt=0;
1026
1027 for (j=0; j < gamma_[i].g; j++) {
1028 nr += rot[j]*gamma_[i].character(j);
1029 nt += trans[j]*gamma_[i].character(j);
1030 }
1031
1032 gamma_[i].nrot_ = (int) ((nr+0.5)/gamma_[i].g);
1033 gamma_[i].ntrans_ = (int) ((nt+0.5)/gamma_[i].g);
1034 }
1035 }
1036
1037 delete[] rot;
1038 delete[] trans;
1039
1040 // now find the inverse of each symop
1041 for (gi=0; gi < g; gi++) {
1042 int gj;
1043 for (gj=0; gj < g; gj++) {
1044 so = symop[gi].operate(symop[gj]);
1045
1046 // is so a unit matrix?
1047 if (fabs(1.0-so[0][0]) < 1.0e-8 &&
1048 fabs(1.0-so[1][1]) < 1.0e-8 &&
1049 fabs(1.0-so[2][2]) < 1.0e-8) break;
1050 }
1051
1052 if (gj==g) {
1053 ExEnv::err0() << indent
1054 << "make_table: uh oh, can't find inverse of " << gi << endl;
1055 abort();
1056 }
1057
1058 _inv[gi] = gj;
1059 }
1060
1061 return 0;
1062}
1063
1064/////////////////////////////////////////////////////////////////////////////
1065
1066// Local Variables:
1067// mode: c++
1068// c-file-style: "ETS"
1069// End:
Note: See TracBrowser for help on using the repository browser.