1 //
2 // pointgrp.cc
3 //
4 // Copyright (C) 1997 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <janssen@limitpt.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 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <util/misc/formio.h>
33 #include <math/symmetry/corrtab.h>
34
35 using namespace std;
36 using namespace sc;
37
38 ////////////////////////////////////////////////////////////////////////
39
CorrelationTable()40 CorrelationTable::CorrelationTable():
41 n_(0),
42 ngamma_(0),
43 gamma_(0)
44 {
45 }
46
CorrelationTable(const Ref<PointGroup> & group,const Ref<PointGroup> & subgroup)47 CorrelationTable::CorrelationTable(const Ref<PointGroup>& group,
48 const Ref<PointGroup>& subgroup):
49 n_(0),
50 ngamma_(0),
51 gamma_(0)
52 {
53 int rc = initialize_table(group,subgroup);
54 if (rc != 0) {
55 ExEnv::err0()
56 << "ERROR: CorrelationTable: " << error(rc) << endl;
57 abort();
58 }
59 }
60
~CorrelationTable()61 CorrelationTable::~CorrelationTable()
62 {
63 clear();
64 }
65
66 int
initialize_table(const Ref<PointGroup> & group,const Ref<PointGroup> & subgroup)67 CorrelationTable::initialize_table(const Ref<PointGroup>& group,
68 const Ref<PointGroup>& subgroup)
69 {
70 clear();
71
72 group_ = group;
73 subgroup_ = subgroup;
74
75 int i, j, k, l;
76
77 CharacterTable ct = group_->char_table();
78 CharacterTable subct = subgroup_->char_table();
79
80 n_ = ct.nirrep();
81 subn_ = subct.nirrep();
82 ngamma_ = new int[n_];
83 gamma_ = new int*[n_];
84
85 // CAN ONLY HANDLE NONDEGENERATE POINT GROUPS
86
87 for (i=0; i<n_; i++) {
88 ngamma_[i] = 0;
89 gamma_[i] = 0;
90 }
91
92 // map the ops in the high order to low order groups
93 int *so_to_subso = new int[ct.order()];
94 int *subso_to_so = new int[subct.order()];
95 for (i=0; i<subct.order(); i++) subso_to_so[i] = -1;
96 for (i=0; i<ct.order(); i++) {
97 SymmetryOperation so = ct.symm_operation(i);
98 int found = 0;
99 so_to_subso[i] = -1;
100 for (j=0; j<subct.order(); j++) {
101 SymmetryOperation subso = subct.symm_operation(j);
102 double sumsquare = 0.0;
103 for (k=0; k<3; k++) {
104 for (l=0; l<3; l++) {
105 double diff = so(k,l)-subso(k,l);
106 sumsquare += diff*diff;
107 }
108 }
109 if (sumsquare < 1.0e-12) {
110 found++;
111 //ExEnv::outn() << scprintf("symmop %d in %s is %d in %s",
112 // i,ct.symbol(),j,subct.symbol()) << endl;
113 so_to_subso[i] = j;
114 subso_to_so[j] = i;
115 }
116 }
117 if (found > 1) {
118 delete[] so_to_subso;
119 delete[] subso_to_so;
120 return -1;
121 }
122 }
123 for (i=0; i<subct.order(); i++) {
124 if (subso_to_so[i] == -1) {
125 delete[] so_to_subso;
126 delete[] subso_to_so;
127 return -2;
128 }
129 }
130
131 // compute the correlation table
132 for (i=0; i<ct.nirrep(); i++) {
133 for (j=0; j<subct.nirrep(); j++) {
134 double nmatch = 0.0;
135 for (k=0; k<ct.order(); k++) {
136 double chr = ct.gamma(i).character(k);
137 if (so_to_subso[k] >= 0) {
138 double subchr = subct.gamma(j).character(so_to_subso[k]);
139 nmatch += subchr*chr;
140 }
141 }
142 nmatch /= subct.order();
143 int inmatch = (int)(nmatch+0.5);
144 if (fabs(nmatch-inmatch)>1.0e-6) {
145 delete[] so_to_subso;
146 delete[] subso_to_so;
147 return -4;
148 }
149 if (inmatch > 0) {
150 int *newgamma = new int[ngamma_[i] + inmatch];
151 memcpy(newgamma,gamma_[i],ngamma_[i]*sizeof(int));
152 for (k=0; k<inmatch; k++) newgamma[ngamma_[i]+k] = j;
153 ngamma_[i] += inmatch;
154 delete[] gamma_[i];
155 gamma_[i] = newgamma;
156 }
157 }
158 }
159
160 delete[] so_to_subso;
161 delete[] subso_to_so;
162
163 for (i=0; i<n(); i++) {
164 int degen = ct.gamma(i).degeneracy();
165 int subdegen = 0;
166 for (j=0; j<ngamma(i); j++) {
167 subdegen += subct.gamma(gamma(i,j)).degeneracy();
168 }
169 if (degen != subdegen) {
170 return -3;
171 }
172 }
173
174 return 0;
175 }
176
177 void
clear()178 CorrelationTable::clear()
179 {
180 for (int i=0; i<n_; i++) {
181 delete[] gamma_[i];
182 }
183 delete[] ngamma_;
184 delete[] gamma_;
185 }
186
187 const char *
error(int rc)188 CorrelationTable::error(int rc)
189 {
190 if (rc == -1) {
191 return "too many symop matches";
192 }
193 else if (rc == -2) {
194 return "not a subgroup or wrong ref frame";
195 }
196 else if (rc == -3) {
197 return "only groups with non-complex characters supported (degen sum)";
198 }
199 else if (rc == -4) {
200 return "only groups with non-complex characters supported (nonint nirr)";
201 }
202 else if (rc != 0) {
203 return "unknown error";
204 }
205 return "no error";
206 }
207
208 int
degen(int i) const209 CorrelationTable::degen(int i) const
210 {
211 return group_->char_table().gamma(i).degeneracy();
212 }
213
214 int
subdegen(int i) const215 CorrelationTable::subdegen(int i) const
216 {
217 return subgroup_->char_table().gamma(i).degeneracy();
218 }
219
220 void
print(ostream & o) const221 CorrelationTable::print(ostream &o) const
222 {
223 o << indent
224 << "Correlation Table from "
225 << group_->symbol() << " to " << subgroup_->symbol() << ":" << endl;
226
227 CharacterTable ct = group_->char_table();
228 CharacterTable subct = subgroup_->char_table();
229
230 o << incindent;
231 for (int i=0; i<n(); i++) {
232 o << indent
233 << ct.gamma(i).symbol() << ":";
234 for (int j=0; j<ngamma(i); j++) {
235 o << indent
236 << " " << subct.gamma(gamma(i,j)).symbol();
237 }
238 o << endl;
239 }
240 o << decindent;
241 }
242
243 /////////////////////////////////////////////////////////////////////////////
244
245 // Local Variables:
246 // mode: c++
247 // c-file-style: "CLJ-CONDENSED"
248 // End:
249