← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/input/find_symmetry.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/input/find_symmetry.cc
 
1
/*! \file
 
2
    \ingroup INPUT
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#define EXTERN
 
6
#include <cstdio>
 
7
#include <libciomr/libciomr.h>
 
8
#include <cstdlib>
 
9
#include <cstring>
 
10
#include <cmath>
 
11
#include "input.h"
 
12
#include "global.h"
 
13
#include "defines.h"
 
14
#include <physconst.h>
 
15
#include <symmetry.h>
 
16
 
 
17
namespace {
 
18
  void alloc_irr_char();
 
19
  void permute_axes(double **, int, int);
 
20
}
 
21
 
 
22
namespace psi { namespace input {
 
23
 
 
24
void find_symmetry()
 
25
{
 
26
  int i, j;
 
27
  double tmp;
 
28
  double *v1;
 
29
  int **redun_atom_orbit;     /*Redundant (i.e. for all 8 operations) atomic orbits*/
 
30
  int c2x_flag = 1;
 
31
  int c2y_flag = 1;
 
32
  int c2z_flag = 1;
 
33
  int inv_flag = 1;
 
34
  int sig_xy_flag = 1;
 
35
  int sig_xz_flag = 1;
 
36
  int sig_yz_flag = 1;
 
37
  int naxes, nplanes;
 
38
  int *row;
 
39
 
 
40
  if (num_atoms == 1) { /* Atomic calculation */
 
41
/*    symmetry = strdup("D2h ");
 
42
    nirreps = 8;
 
43
    irr_labels = (char **) malloc(sizeof(char *)*nirreps);
 
44
    sym_oper = init_int_array(nirreps);
 
45
    alloc_irr_char();
 
46
    irr_labels[0] = strdup("Ag  ");
 
47
    irr_labels[1] = strdup("B1g ");
 
48
    irr_labels[2] = strdup("B2g ");
 
49
    irr_labels[3] = strdup("B3g ");
 
50
    irr_labels[4] = strdup("Au  ");
 
51
    irr_labels[5] = strdup("B1u ");
 
52
    irr_labels[6] = strdup("B2u ");
 
53
    irr_labels[7] = strdup("B3u ");
 
54
    sym_oper[1] = C2ZFLAG;
 
55
    sym_oper[2] = C2YFLAG;
 
56
    sym_oper[3] = C2XFLAG;
 
57
    sym_oper[4] = IFLAG;
 
58
    sym_oper[5] = SIGXYFLAG;
 
59
    sym_oper[6] = SIGXZFLAG;
 
60
    sym_oper[7] = SIGYZFLAG;
 
61
    atom_orbit = init_int_matrix(1,nirreps);
 
62
    atom_orbit[0][0] = atom_orbit[0][1] = atom_orbit[0][2] = atom_orbit[0][3] = atom_orbit[0][4] =
 
63
    atom_orbit[0][5] = atom_orbit[0][6] = atom_orbit[0][7] = 0;
 
64
    return; */
 
65
  }
 
66
 
 
67
  v1 = init_array(3);
 
68
  redun_atom_orbit = init_int_matrix(num_atoms,8);
 
69
 
 
70
  for(i=0;i<num_atoms;i++) {
 
71
    if (c2x_flag && !redun_atom_orbit[i][C2XFLAG]) {
 
72
      c2x(geometry[i],v1);
 
73
      c2x_flag = 0;
 
74
      for(j=0;j<num_atoms;j++)
 
75
        if ((int)nuclear_charges[i] == (int)nuclear_charges[j] && vect_equiv(geometry[j],v1)){
 
76
          redun_atom_orbit[i][C2XFLAG] = j;
 
77
          redun_atom_orbit[j][C2XFLAG] = i;
 
78
          c2x_flag = 1;
 
79
          break;
 
80
        }
 
81
    }
 
82
    if (c2y_flag && !redun_atom_orbit[i][C2YFLAG]) {
 
83
      c2y(geometry[i],v1);
 
84
      c2y_flag = 0;
 
85
      for(j=0;j<num_atoms;j++)
 
86
        if ((int)nuclear_charges[i] == (int)nuclear_charges[j] && vect_equiv(geometry[j],v1)){
 
87
          redun_atom_orbit[i][C2YFLAG] = j;
 
88
          redun_atom_orbit[j][C2YFLAG] = i;
 
89
          c2y_flag = 1;
 
90
          break;
 
91
        }
 
92
    }
 
93
    if (c2z_flag && !redun_atom_orbit[i][C2ZFLAG]) {
 
94
      c2z(geometry[i],v1);
 
95
      c2z_flag = 0;
 
96
      for(j=0;j<num_atoms;j++)
 
97
        if ((int)nuclear_charges[i] == (int)nuclear_charges[j] && vect_equiv(geometry[j],v1)){
 
98
          redun_atom_orbit[i][C2ZFLAG] = j;
 
99
          redun_atom_orbit[j][C2ZFLAG] = i;
 
100
          c2z_flag = 1;
 
101
          break;
 
102
        }
 
103
    }
 
104
    if (inv_flag && !redun_atom_orbit[i][IFLAG]) {
 
105
      inversion(geometry[i],v1);
 
106
      inv_flag = 0;
 
107
      for(j=0;j<num_atoms;j++)
 
108
        if ((int)nuclear_charges[i] == (int)nuclear_charges[j] && vect_equiv(geometry[j],v1)){
 
109
          redun_atom_orbit[i][IFLAG] = j;
 
110
          redun_atom_orbit[j][IFLAG] = i;
 
111
          inv_flag = 1;
 
112
          break;
 
113
        }
 
114
    }
 
115
    if (sig_xy_flag && !redun_atom_orbit[i][SIGXYFLAG]) {
 
116
      sig_xy(geometry[i],v1);
 
117
      sig_xy_flag = 0;
 
118
      for(j=0;j<num_atoms;j++)
 
119
        if ((int)nuclear_charges[i] == (int)nuclear_charges[j] && vect_equiv(geometry[j],v1)){
 
120
          redun_atom_orbit[i][SIGXYFLAG] = j;
 
121
          redun_atom_orbit[j][SIGXYFLAG] = i;
 
122
          sig_xy_flag = 1;
 
123
          break;
 
124
        }
 
125
    }
 
126
    if (sig_xz_flag && !redun_atom_orbit[i][SIGXZFLAG]) {
 
127
      sig_xz(geometry[i],v1);
 
128
      sig_xz_flag = 0;
 
129
      for(j=0;j<num_atoms;j++)
 
130
        if ((int)nuclear_charges[i] == (int)nuclear_charges[j] && vect_equiv(geometry[j],v1)){
 
131
          redun_atom_orbit[i][SIGXZFLAG] = j;
 
132
          redun_atom_orbit[j][SIGXZFLAG] = i;
 
133
          sig_xz_flag = 1;
 
134
          break;
 
135
        }
 
136
    }
 
137
    if (sig_yz_flag && !redun_atom_orbit[i][SIGYZFLAG]) {
 
138
      sig_yz(geometry[i],v1);
 
139
      sig_yz_flag = 0;
 
140
      for(j=0;j<num_atoms;j++)
 
141
        if ((int)nuclear_charges[i] == (int)nuclear_charges[j] && vect_equiv(geometry[j],v1)){
 
142
          redun_atom_orbit[i][SIGYZFLAG] = j;
 
143
          redun_atom_orbit[j][SIGYZFLAG] = i;
 
144
          sig_yz_flag = 1;
 
145
          break;
 
146
        }
 
147
    }
 
148
  }
 
149
 
 
150
  free(v1);
 
151
 
 
152
  if (print_lvl >= DEBUGPRINT) {
 
153
    fprintf(outfile,"\n  -Symmetry flags:\n");
 
154
    fprintf(outfile,"    C2x = %d\n",c2x_flag);
 
155
    fprintf(outfile,"    C2y = %d\n",c2y_flag);
 
156
    fprintf(outfile,"    C2z = %d\n",c2z_flag);
 
157
    fprintf(outfile,"    I   = %d\n",inv_flag);
 
158
    fprintf(outfile,"    sxy = %d\n",sig_xy_flag);
 
159
    fprintf(outfile,"    sxz = %d\n",sig_xz_flag);
 
160
    fprintf(outfile,"    syz = %d\n",sig_yz_flag);
 
161
  }
 
162
  naxes = c2x_flag + c2y_flag + c2z_flag;
 
163
  nplanes = sig_xy_flag + sig_xz_flag + sig_yz_flag;
 
164
  nirreps = 1 + naxes + inv_flag + nplanes;
 
165
 
 
166
 
 
167
  /*-------------------
 
168
    Handling subgroups
 
169
   --------------------*/
 
170
 
 
171
  if (subgroup != NULL) {
 
172
    switch (nirreps) {
 
173
      case 8: /*Can be C2v, C2h, D2, etc.*/
 
174
              if (!strcmp(subgroup,"D2")) {
 
175
                nirreps = 4; nplanes = 0; inv_flag = 0;
 
176
                sig_xy_flag = 0; sig_xz_flag = 0; sig_yz_flag = 0;
 
177
                break;
 
178
              }
 
179
              if (!strcmp(subgroup,"C2V")) {
 
180
                if (unique_axis == NULL) {
 
181
                  printf("  UNIQUE_AXIS must be specified for this symmetry.\n  The largest Abelian point group will be used.\n\n");
 
182
                  break;
 
183
                }
 
184
                nirreps = 4; naxes = 1; nplanes = 2; inv_flag = 0;
 
185
                if (!strcmp(unique_axis,"X")) {
 
186
                  c2y_flag = 0; c2z_flag = 0; sig_yz_flag = 0;
 
187
                }
 
188
                else if (!strcmp(unique_axis,"Y")) {
 
189
                  c2x_flag = 0; c2z_flag = 0; sig_xz_flag = 0;
 
190
                }
 
191
                else if (!strcmp(unique_axis,"Z")) {
 
192
                  c2x_flag = 0; c2y_flag = 0; sig_xy_flag = 0;
 
193
                }
 
194
                break;
 
195
              }
 
196
              if (!strcmp(subgroup,"C2H")) {
 
197
                if (unique_axis == NULL) {
 
198
                  printf("  UNIQUE_AXIS must be specified for this symmetry.\n  The largest Abelian point group will be used.\n\n");
 
199
                  break;
 
200
                }
 
201
                nirreps = 4; naxes = 1; nplanes = 1; inv_flag = 1;
 
202
                if (!strcmp(unique_axis,"X")) {
 
203
                  c2y_flag = 0; c2z_flag = 0; sig_xy_flag = 0; sig_xz_flag = 0;
 
204
                }
 
205
                else if (!strcmp(unique_axis,"Y")) {
 
206
                  c2x_flag = 0; c2z_flag = 0; sig_xy_flag = 0; sig_yz_flag = 0;
 
207
                }
 
208
                else if (!strcmp(unique_axis,"Z")) {
 
209
                  c2x_flag = 0; c2y_flag = 0; sig_xz_flag = 0; sig_yz_flag = 0;
 
210
                }
 
211
                break;
 
212
              }
 
213
      case 4: /*Can be C2, Cs, Ci, etc.*/
 
214
              if (inv_flag && !strcmp(subgroup,"CI")) {
 
215
                nirreps = 2; naxes = 0; nplanes = 0;
 
216
                c2x_flag = 0; c2y_flag = 0; c2z_flag = 0;
 
217
                sig_xy_flag = 0; sig_xz_flag = 0; sig_yz_flag = 0;
 
218
                break;
 
219
              }
 
220
              if ((naxes > 0) && !strcmp(subgroup,"C2")) {
 
221
                if (naxes == 3) {
 
222
                  if (unique_axis == NULL) {
 
223
                    printf("  UNIQUE_AXIS must be specified for this symmetry.\n  The largest Abelian point group will be used.\n\n");
 
224
                    break;
 
225
                  }
 
226
                  if (!strcmp(unique_axis,"X")) {
 
227
                    c2y_flag = 0; c2z_flag = 0;
 
228
                  }
 
229
                  else if (!strcmp(unique_axis,"Y")) {
 
230
                    c2x_flag = 0; c2z_flag = 0;
 
231
                  }
 
232
                  else if (!strcmp(unique_axis,"Z")) {
 
233
                    c2x_flag = 0; c2y_flag = 0;
 
234
                  }
 
235
                  naxes = 1;
 
236
                }
 
237
                nirreps = 2; nplanes = 0; inv_flag = 0;
 
238
                sig_xy_flag = 0; sig_xz_flag = 0; sig_yz_flag = 0;
 
239
                break;
 
240
              }
 
241
              if ((nplanes > 0) && !strcmp(subgroup,"CS")) {
 
242
                if (nplanes > 1) {
 
243
                  if (unique_axis == NULL) {
 
244
                    printf("  UNIQUE_AXIS must be specified for this symmetry.\n  The largest Abelian point group will be used.\n\n");
 
245
                    break;
 
246
                  }
 
247
                  if (!strcmp(unique_axis,"Z"))
 
248
                    if (sig_xy_flag) {
 
249
                      sig_xz_flag = 0; sig_yz_flag = 0;
 
250
                    }
 
251
                    else
 
252
                      punt("  UNIQUE_AXIS defined incorrectly.");
 
253
                  else if (!strcmp(unique_axis,"Y"))
 
254
                    if (sig_xz_flag) {
 
255
                      sig_xy_flag = 0; sig_yz_flag = 0;
 
256
                    }
 
257
                    else
 
258
                      punt("  UNIQUE_AXIS defined incorrectly.");
 
259
                  else if (!strcmp(unique_axis,"X"))
 
260
                    if (sig_yz_flag) {
 
261
                      sig_xy_flag = 0; sig_xz_flag = 0;
 
262
                    }
 
263
                    else
 
264
                      punt("  UNIQUE_AXIS defined incorrectly.");
 
265
                  nplanes = 1;
 
266
                }
 
267
                nirreps = 2; naxes = 0; inv_flag = 0;
 
268
                c2x_flag = 0; c2y_flag = 0; c2z_flag = 0;
 
269
                break;
 
270
              }
 
271
      case 2: /*Can only be C1*/
 
272
              if (!strcmp(subgroup,"C1")) {
 
273
                nirreps = 1; naxes = 0; nplanes = 0;
 
274
                c2x_flag = 0; c2y_flag = 0; c2z_flag = 0; inv_flag = 0;
 
275
                sig_xy_flag = 0; sig_xz_flag = 0; sig_yz_flag = 0;
 
276
                break;
 
277
              }
 
278
      case 1: /*What subgroup???*/
 
279
              if (strcmp(subgroup,"C1")) {
 
280
                printf("  Subgroup specification doesn't comply with the actual symmetry.\n");
 
281
                printf("  The largest Abelian point group will be used.\n\n");
 
282
              }
 
283
              break;
 
284
    }
 
285
    if (print_lvl >= DEBUGPRINT) {
 
286
      fprintf(outfile,"\n  -Flags within the subgroup:\n");
 
287
      fprintf(outfile,"    C2x = %d\n",c2x_flag);
 
288
      fprintf(outfile,"    C2y = %d\n",c2y_flag);
 
289
      fprintf(outfile,"    C2z = %d\n",c2z_flag);
 
290
      fprintf(outfile,"    I   = %d\n",inv_flag);
 
291
      fprintf(outfile,"    sxy = %d\n",sig_xy_flag);
 
292
      fprintf(outfile,"    sxz = %d\n",sig_xz_flag);
 
293
      fprintf(outfile,"    syz = %d\n",sig_yz_flag);
 
294
    }
 
295
  }
 
296
 
 
297
  
 
298
  /*-------------------------
 
299
    Allocate character table
 
300
   -------------------------*/
 
301
 
 
302
  alloc_irr_char();
 
303
  irr_labels = (char **) malloc(sizeof(char *)*nirreps);
 
304
 
 
305
  
 
306
  /*---------------------------------------------------------------------
 
307
    Figure out point group, reorient, if necessary, and get rid of zeros
 
308
    in the array of atomic orbits
 
309
   ---------------------------------------------------------------------*/
 
310
  
 
311
  atom_orbit = init_int_matrix(num_atoms,nirreps);
 
312
  for(i=0;i<num_atoms;i++)
 
313
    atom_orbit[i][EFLAG] = i;
 
314
  sym_oper = init_int_array(nirreps);
 
315
  sym_oper[0] = EFLAG;
 
316
  switch (nirreps) {
 
317
    case 1: /* There's one irreducible representation - it is C1 */
 
318
            symmetry = strdup("C1  ");
 
319
            irr_labels[0] = strdup("A   ");
 
320
            break;
 
321
    case 2: /* Could be Ci, C2, or Cs */
 
322
            if (inv_flag) {
 
323
              symmetry = strdup("Ci  ");
 
324
              irr_labels[0] = strdup("Ag  ");
 
325
              irr_labels[1] = strdup("Au  ");
 
326
              sym_oper[1] = IFLAG;
 
327
              for(i=0;i<num_atoms;i++)
 
328
                atom_orbit[i][1] = redun_atom_orbit[i][IFLAG]; /* i */
 
329
            }
 
330
            else if (naxes) {
 
331
              symmetry = strdup("C2  ");
 
332
              irr_labels[0] = strdup("A   ");
 
333
              irr_labels[1] = strdup("B   ");
 
334
              sym_oper[1] = C2ZFLAG;
 
335
              /* Making sure that the symmetry axis is along Z */
 
336
              if (c2z_flag)
 
337
                for(i=0;i<num_atoms;i++)
 
338
                  atom_orbit[i][1] = redun_atom_orbit[i][C2ZFLAG]; /* C2 = C2z */
 
339
              else
 
340
                if (c2x_flag) { /* Unique axis is along X */
 
341
                  permute_axes(geometry,XAXIS,ZAXIS);
 
342
                  for(i=0;i<num_atoms;i++) {
 
343
                    atom_orbit[i][1] = redun_atom_orbit[i][C2XFLAG]; /* C2 = C2x */
 
344
                  }
 
345
                }
 
346
                else { /* Unique axis is along Y */
 
347
                  permute_axes(geometry,YAXIS,ZAXIS);
 
348
                  for(i=0;i<num_atoms;i++) {
 
349
                    atom_orbit[i][1] = redun_atom_orbit[i][C2YFLAG]; /* C2 = C2y */
 
350
                  }
 
351
                }
 
352
            }
 
353
            else if (nplanes) { /* Cs symmetry */
 
354
              symmetry = strdup("Cs  ");
 
355
              irr_labels[0] = strdup("Ap  ");
 
356
              irr_labels[1] = strdup("App ");
 
357
              sym_oper[1] = SIGXYFLAG;
 
358
              if (sig_xy_flag)
 
359
                for(i=0;i<num_atoms;i++)
 
360
                  atom_orbit[i][1] = redun_atom_orbit[i][SIGXYFLAG]; /* sig_xy */
 
361
              else
 
362
                if (sig_yz_flag) { /* Unique axis is along X */
 
363
                  permute_axes(geometry,XAXIS,ZAXIS);
 
364
                  for(i=0;i<num_atoms;i++) {
 
365
                    atom_orbit[i][1] = redun_atom_orbit[i][SIGYZFLAG]; /* sig_xy = sig_yz */
 
366
                  }
 
367
                }
 
368
                else { /* Unique axis is along Y */
 
369
                  permute_axes(geometry,YAXIS,ZAXIS);
 
370
                  for(i=0;i<num_atoms;i++) {
 
371
                    atom_orbit[i][1] = redun_atom_orbit[i][SIGXZFLAG]; /* sig_xy = sig_xz */
 
372
                  }
 
373
                }
 
374
            }
 
375
            else
 
376
              punt("Unidentified symmetry.");
 
377
            break;
 
378
            
 
379
    case 4: /* It could be C2v, C2h, or D2 */
 
380
            if (inv_flag) {
 
381
              symmetry = strdup("C2h ");
 
382
              /* Swap rows 1 and 2 in irr_char since ordering of irreps in C2h is different from that in C2v and D2*/
 
383
              row = irr_char[1];
 
384
              irr_char[1] = irr_char[2];
 
385
              irr_char[2] = row;
 
386
              irr_labels[0] = strdup("Ag  ");
 
387
              irr_labels[1] = strdup("Bg  ");
 
388
              irr_labels[2] = strdup("Au  ");
 
389
              irr_labels[3] = strdup("Bu  ");
 
390
              sym_oper[1] = C2ZFLAG;
 
391
              sym_oper[2] = IFLAG;
 
392
              sym_oper[3] = SIGXYFLAG;
 
393
              if (c2z_flag)
 
394
                for(i=0;i<num_atoms;i++) {
 
395
                  atom_orbit[i][1] = redun_atom_orbit[i][C2ZFLAG]; /* C2 */
 
396
                  atom_orbit[i][2] = redun_atom_orbit[i][IFLAG]; /* i */
 
397
                  atom_orbit[i][3] = redun_atom_orbit[i][SIGXYFLAG]; /* sig_xy */
 
398
                }
 
399
              else
 
400
                if (c2x_flag) { /* Unique axis is along X */
 
401
                  permute_axes(geometry,XAXIS,ZAXIS);
 
402
                  for(i=0;i<num_atoms;i++) {
 
403
                    atom_orbit[i][1] = redun_atom_orbit[i][C2XFLAG]; /* C2 = C2x */
 
404
                    atom_orbit[i][2] = redun_atom_orbit[i][IFLAG]; /* i */
 
405
                    atom_orbit[i][3] = redun_atom_orbit[i][SIGYZFLAG]; /* sig_xy = sig_yz */
 
406
                  }
 
407
                }
 
408
                else { /* Unique axis is along Y */
 
409
                  permute_axes(geometry,YAXIS,ZAXIS);
 
410
                  for(i=0;i<num_atoms;i++) {
 
411
                    atom_orbit[i][1] = redun_atom_orbit[i][C2YFLAG]; /* C2 = C2y */
 
412
                    atom_orbit[i][2] = redun_atom_orbit[i][IFLAG]; /* i */
 
413
                    atom_orbit[i][3] = redun_atom_orbit[i][SIGXZFLAG]; /* sig_xy = sig_xz */
 
414
                  }
 
415
                }
 
416
            }
 
417
            else if (nplanes == 2) {
 
418
              symmetry = strdup("C2v ");
 
419
              irr_labels[0] = strdup("A1  ");
 
420
              irr_labels[1] = strdup("A2  ");
 
421
              irr_labels[2] = strdup("B1  ");
 
422
              irr_labels[3] = strdup("B2  ");
 
423
              sym_oper[1] = C2ZFLAG;
 
424
              sym_oper[2] = SIGXZFLAG;
 
425
              sym_oper[3] = SIGYZFLAG;
 
426
              if (c2z_flag)
 
427
                for(i=0;i<num_atoms;i++) {
 
428
                  atom_orbit[i][1] = redun_atom_orbit[i][C2ZFLAG]; /* C2 */
 
429
                  atom_orbit[i][2] = redun_atom_orbit[i][SIGXZFLAG]; /* sig_xz */
 
430
                  atom_orbit[i][3] = redun_atom_orbit[i][SIGYZFLAG]; /* sig_yz */
 
431
                }
 
432
              else
 
433
                if (c2x_flag) { /* Unique axis is along X */
 
434
                  permute_axes(geometry,XAXIS,ZAXIS);
 
435
                  for(i=0;i<num_atoms;i++) {
 
436
                    atom_orbit[i][1] = redun_atom_orbit[i][C2XFLAG]; /* C2 = C2x */
 
437
                    atom_orbit[i][2] = redun_atom_orbit[i][SIGXYFLAG]; /* sig_xz = sig_yx */
 
438
                    atom_orbit[i][3] = redun_atom_orbit[i][SIGXZFLAG]; /* sig_yz = sig_zx */
 
439
                  }
 
440
                }
 
441
                else { /* Unique axis is along Y */
 
442
                  permute_axes(geometry,YAXIS,ZAXIS);
 
443
                  for(i=0;i<num_atoms;i++) {
 
444
                    atom_orbit[i][1] = redun_atom_orbit[i][C2YFLAG]; /* C2 = C2y */
 
445
                    atom_orbit[i][2] = redun_atom_orbit[i][SIGYZFLAG]; /* sig_xz = sig_zy */
 
446
                    atom_orbit[i][3] = redun_atom_orbit[i][SIGXYFLAG]; /* sig_yz = sig_xy */
 
447
                  }
 
448
                }
 
449
            }
 
450
            else if (naxes == 3) { /* D2 symmetry */
 
451
              symmetry = strdup("D2  ");
 
452
              irr_labels[0] = strdup("A   ");
 
453
              irr_labels[1] = strdup("B1  ");
 
454
              irr_labels[2] = strdup("B2  ");
 
455
              irr_labels[3] = strdup("B3  ");
 
456
              sym_oper[1] = C2ZFLAG;
 
457
              sym_oper[2] = C2YFLAG;
 
458
              sym_oper[3] = C2XFLAG;
 
459
              for(i=0;i<num_atoms;i++) {
 
460
                atom_orbit[i][1] = redun_atom_orbit[i][C2ZFLAG]; /* C2z */
 
461
                atom_orbit[i][2] = redun_atom_orbit[i][C2YFLAG]; /* C2y */
 
462
                atom_orbit[i][3] = redun_atom_orbit[i][C2XFLAG]; /* C2x */
 
463
              }
 
464
            }
 
465
            else
 
466
              punt("Unidentified symmetry.");
 
467
            break;
 
468
            
 
469
    case 8: /* D2h symmetry */
 
470
            if (nplanes == 3 && naxes == 3 && inv_flag) {
 
471
              symmetry = strdup("D2h ");
 
472
              irr_labels[0] = strdup("Ag  ");
 
473
              irr_labels[1] = strdup("B1g ");
 
474
              irr_labels[2] = strdup("B2g ");
 
475
              irr_labels[3] = strdup("B3g ");
 
476
              irr_labels[4] = strdup("Au  ");
 
477
              irr_labels[5] = strdup("B1u ");
 
478
              irr_labels[6] = strdup("B2u ");
 
479
              irr_labels[7] = strdup("B3u ");
 
480
              sym_oper[1] = C2ZFLAG;
 
481
              sym_oper[2] = C2YFLAG;
 
482
              sym_oper[3] = C2XFLAG;
 
483
              sym_oper[4] = IFLAG;
 
484
              sym_oper[5] = SIGXYFLAG;
 
485
              sym_oper[6] = SIGXZFLAG;
 
486
              sym_oper[7] = SIGYZFLAG;
 
487
              for(i=0;i<num_atoms;i++) {
 
488
                atom_orbit[i][1] = redun_atom_orbit[i][C2ZFLAG]; /* C2z */
 
489
                atom_orbit[i][2] = redun_atom_orbit[i][C2YFLAG]; /* C2y */
 
490
                atom_orbit[i][3] = redun_atom_orbit[i][C2XFLAG]; /* C2x */
 
491
                atom_orbit[i][4] = redun_atom_orbit[i][IFLAG]; /* i */
 
492
                atom_orbit[i][5] = redun_atom_orbit[i][SIGXYFLAG]; /* sig_xy */
 
493
                atom_orbit[i][6] = redun_atom_orbit[i][SIGXZFLAG]; /* sig_xz */
 
494
                atom_orbit[i][7] = redun_atom_orbit[i][SIGYZFLAG]; /* sig_yz */
 
495
              }
 
496
            }
 
497
            else
 
498
              punt("Unidentified symmetry.");
 
499
            break;
 
500
            
 
501
    default: /* Error */
 
502
            punt("Unidentified symmetry.");
 
503
  } /* end of switch(nirreps) */
 
504
 
 
505
  free_int_matrix(redun_atom_orbit);
 
506
  
 
507
  return;
 
508
}
 
509
  
 
510
}} // namespace psi::input
 
511
 
 
512
namespace {
 
513
 
 
514
using namespace psi::input;
 
515
 
 
516
void alloc_irr_char()
 
517
{
 
518
  int i;
 
519
 
 
520
  /* Allocate global arrays */
 
521
  irr_char = init_int_matrix(nirreps,nirreps);
 
522
 
 
523
  for(i=0;i<nirreps;i++)
 
524
    irr_char[0][i] = 1;
 
525
  switch (nirreps) {
 
526
    case 1:
 
527
      break;
 
528
    case 2:
 
529
      irr_char[1][0] =  1; irr_char[1][1] = -1;
 
530
      break;
 
531
    case 4:
 
532
      irr_char[1][0] =  1; irr_char[1][1] =  1; irr_char[1][2] = -1; irr_char[1][3] = -1;
 
533
      irr_char[2][0] =  1; irr_char[2][1] = -1; irr_char[2][2] =  1; irr_char[2][3] = -1;
 
534
      irr_char[3][0] =  1; irr_char[3][1] = -1; irr_char[3][2] = -1; irr_char[3][3] =  1;
 
535
      break;
 
536
    case 8:
 
537
      irr_char[1][0] =  1; irr_char[1][1] =  1; irr_char[1][2] = -1; irr_char[1][3] = -1;
 
538
      irr_char[1][4] =  1; irr_char[1][5] =  1; irr_char[1][6] = -1; irr_char[1][7] = -1;
 
539
      irr_char[2][0] =  1; irr_char[2][1] = -1; irr_char[2][2] =  1; irr_char[2][3] = -1;
 
540
      irr_char[2][4] =  1; irr_char[2][5] = -1; irr_char[2][6] =  1; irr_char[2][7] = -1;
 
541
      irr_char[3][0] =  1; irr_char[3][1] = -1; irr_char[3][2] = -1; irr_char[3][3] =  1;
 
542
      irr_char[3][4] =  1; irr_char[3][5] = -1; irr_char[3][6] = -1; irr_char[3][7] =  1;
 
543
      irr_char[4][0] =  1; irr_char[4][1] =  1; irr_char[4][2] =  1; irr_char[4][3] =  1;
 
544
      irr_char[4][4] = -1; irr_char[4][5] = -1; irr_char[4][6] = -1; irr_char[4][7] = -1;
 
545
      irr_char[5][0] =  1; irr_char[5][1] =  1; irr_char[5][2] = -1; irr_char[5][3] = -1;
 
546
      irr_char[5][4] = -1; irr_char[5][5] = -1; irr_char[5][6] =  1; irr_char[5][7] =  1;
 
547
      irr_char[6][0] =  1; irr_char[6][1] = -1; irr_char[6][2] =  1; irr_char[6][3] = -1;
 
548
      irr_char[6][4] = -1; irr_char[6][5] =  1; irr_char[6][6] = -1; irr_char[6][7] =  1;
 
549
      irr_char[7][0] =  1; irr_char[7][1] = -1; irr_char[7][2] = -1; irr_char[7][3] =  1;
 
550
      irr_char[7][4] = -1; irr_char[7][5] =  1; irr_char[7][6] =  1; irr_char[7][7] = -1;
 
551
      break;
 
552
    default:
 
553
      punt("Number of irreps in find_symmetry is invalid.");
 
554
  }
 
555
 
 
556
  return;
 
557
}
 
558
 
 
559
/* This function permutes axes so that axis1 becomes axis2
 
560
   This is done by a rotation about C3 axis relating x, y, and z
 
561
   unit vectors */
 
562
void permute_axes(double **geom, int axis1, int axis2)
 
563
{
 
564
  int atom,i,j;
 
565
  int axis3;
 
566
  double tmp;
 
567
  double **R;    /* The matrix representation of the forward rotation */
 
568
 
 
569
  if (axis1 == axis2)
 
570
      punt("Called permute_axes with axis1=axis2");
 
571
 
 
572
  for(i=0;i<3;i++)
 
573
    if ( (i != axis1) && (i != axis2)) {
 
574
      axis3 = i;
 
575
      break;
 
576
    }
 
577
 
 
578
  R = block_matrix(3,3);
 
579
  R[axis1][axis2] = 1.0;
 
580
  R[axis2][axis3] = 1.0;
 
581
  R[axis3][axis1] = 1.0;
 
582
  rotate_full_geom(R);
 
583
  free_block(R);
 
584
 
 
585
  return;
 
586
}
 
587
 
 
588
} // namespace