MADNESS  0.10.1
pointgroup.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31  $Id$
32 */
33 #ifndef MAD_POINTGROUP_H
34 #define MAD_POINTGROUP_H
35 
36 /// \file pointgroup.h
37 /// \brief Implements basic functionality for Abelian point groups (D2h and subgroups)
38 
39 #include <madness/tensor/tensor.h>
40 #include <madness/world/vector.h>
41 #include <string>
42 #include <ostream>
43 #include <sstream>
44 #include <iomanip>
45 #include <algorithm>
46 
47 class PointGroup;
48 std::ostream& operator<<(std::ostream& s, const PointGroup& g);
49 
50 class PointGroup {
51  std::string name; ///< group name
52  int order; ///< group order
53  std::string irs[8]; ///< Names of the irreps
54  std::string ops[8]; ///< Names of the operators
55  int c[8][8]; ///< Character table
56 
57 public:
59 
60  /// Constructs point group by name (D2h and subgroups only)
61  PointGroup(const std::string name) {
62  this->name = name;
63  if (name == "C1") {
64  order = 1;
65  irs[0] = "a";
66  ops[0] = "e";
67  c[0][0] = 1;
68  }
69  else if (name == "C2") {
70  order = 2;
71  irs[0] = "a"; irs[1] = "b";
72  ops[0] = "e"; ops[1] = "c2z";
73  c[0][0] = 1; c[0][1] = 1;
74  c[1][0] = 1; c[1][1] =-1;
75  }
76  else if (name == "Ci") {
77  order = 2;
78  irs[0] = "ag"; irs[1] = "au";
79  ops[0] = "e"; ops[1] = "i";
80  c[0][0] = 1; c[0][1] = 1;
81  c[1][0] = 1; c[1][1] =-1;
82  }
83  else if (name == "Cs") {
84  order = 2;
85  irs[0] = "a"; irs[1] = "a'";
86  ops[0] = "e"; ops[1] = "sxy";
87  c[0][0] = 1; c[0][1] = 1;
88  c[1][0] = 1; c[1][1] =-1;
89  }
90  else if (name == "C2h") {
91  order = 4;
92  irs[0] = "ag"; irs[1] = "au"; irs[2] = "bg"; irs[3] = "bu";
93  ops[0] = "e"; ops[1] = "c2z";ops[2] = "sxy";ops[3] = "i";
94  c[0][0] = 1; c[0][1] = 1; c[0][2] = 1; c[0][3] = 1;
95  c[1][0] = 1; c[1][1] = 1; c[1][2] =-1; c[1][3] =-1;
96  c[2][0] = 1; c[2][1] =-1; c[2][2] =-1; c[2][3] = 1;
97  c[3][0] = 1; c[3][1] =-1; c[3][2] = 1; c[3][3] =-1;
98  }
99  else if (name == "C2v") {
100  order = 4;
101  irs[0] = "a1"; irs[1] = "a2"; irs[2] = "b1"; irs[3] = "b2";
102  ops[0] = "e"; ops[1] = "c2z";ops[2] = "sxz";ops[3] = "syz";
103  c[0][0] = 1; c[0][1] = 1; c[0][2] = 1; c[0][3] = 1;
104  c[1][0] = 1; c[1][1] = 1; c[1][2] =-1; c[1][3] =-1;
105  c[2][0] = 1; c[2][1] =-1; c[2][2] = 1; c[2][3] =-1;
106  c[3][0] = 1; c[3][1] =-1; c[3][2] =-1; c[3][3] = 1;
107  }
108  else if (name == "D2") {
109  order = 4;
110  irs[0] = "a1"; irs[1] = "b1"; irs[2] = "b2"; irs[3] = "b3";
111  ops[0] = "e"; ops[1] = "c2z";ops[2] = "c2y";ops[3] = "c2x";
112  c[0][0] = 1; c[0][1] = 1; c[0][2] = 1; c[0][3] = 1;
113  c[1][0] = 1; c[1][1] = 1; c[1][2] =-1; c[1][3] =-1;
114  c[2][0] = 1; c[2][1] =-1; c[2][2] = 1; c[2][3] =-1;
115  c[3][0] = 1; c[3][1] =-1; c[3][2] =-1; c[3][3] = 1;
116  }
117  else if (name == "D2h") {
118  order = 8;
119  irs[0] = "ag"; irs[1] = "au"; irs[2] = "b1g"; irs[3] = "b1u";
120  irs[4] = "b2g";irs[5] = "b2u";irs[6] = "b3g"; irs[7] = "b3u";
121  ops[0] = "e"; ops[1] = "c2z";ops[2] = "c2y"; ops[3] = "c2x";
122  ops[4] = "i"; ops[5] = "sxy";ops[6] = "sxz"; ops[7] = "syz";
123  c[0][0] = 1; c[0][1] = 1; c[0][2] = 1; c[0][3] = 1; c[0][4] = 1; c[0][5] = 1; c[0][6] = 1; c[0][7] = 1;
124  c[1][0] = 1; c[1][1] = 1; c[1][2] = 1; c[1][3] = 1; c[1][4] =-1; c[1][5] =-1; c[1][6] =-1; c[1][7] =-1;
125  c[2][0] = 1; c[2][1] = 1; c[2][2] =-1; c[2][3] =-1; c[2][4] = 1; c[2][5] = 1; c[2][6] =-1; c[2][7] =-1;
126  c[3][0] = 1; c[3][1] = 1; c[3][2] =-1; c[3][3] =-1; c[3][4] =-1; c[3][5] =-1; c[3][6] = 1; c[3][7] = 1;
127  c[4][0] = 1; c[4][1] =-1; c[4][2] = 1; c[4][3] =-1; c[4][4] = 1; c[4][5] =-1; c[4][6] = 1; c[4][7] =-1;
128  c[5][0] = 1; c[5][1] =-1; c[5][2] = 1; c[5][3] =-1; c[5][4] =-1; c[5][5] = 1; c[5][6] =-1; c[5][7] = 1;
129  c[6][0] = 1; c[6][1] =-1; c[6][2] =-1; c[6][3] = 1; c[6][4] = 1; c[6][5] =-1; c[6][6] =-1; c[6][7] = 1;
130  c[7][0] = 1; c[7][1] =-1; c[7][2] =-1; c[7][3] = 1; c[7][4] =-1; c[7][5] = 1; c[7][6] = 1; c[7][7] =-1;
131  }
132  else {
133  throw "PointGroup: unknown group";
134  }
135  }
136 
137  /// Copy constructor
138  PointGroup(const PointGroup& other) {
139  *this = other;
140  }
141 
142  /// Assignment
143  PointGroup& operator=(const PointGroup& other) {
144  if (this != &other) {
145  name = other.name;
146  order = other.order;
147  for (int ir=0; ir<order; ++ir) {
148  irs[ir] = other.irs[ir];
149  ops[ir] = other.ops[ir];
150  for (int op=0; op<order; ++op) {
151  c[ir][op] = other.c[ir][op];
152  }
153  }
154  }
155  return *this;
156  }
157 
158  /// Destructor
159  virtual ~PointGroup() {}
160 
161  /// Returns irreducible representation corresponding to product
162  int irmul(int ir1, int ir2) const {
163  return ir1 ^ ir2;
164  }
165 
166  /// Applies group operator number op (0,1,...,order-1) to point
167  coordT apply(int op, const coordT& r) const {
168  return apply(ops[op], r);
169  }
170 
171  /// Applies named operator (e, c2z, c2y, c2x, sxy, sxz, syz, i) to point
172  static coordT apply(const std::string& op, const coordT& r) {
173  const double x=r[0], y=r[1], z=r[2];
174  coordT q;
175  if (op == "e") {
176  q[0]=x; q[1]=y; q[2]=z;
177  }
178  else if (op == "c2z") {
179  q[0]=-x; q[1]=-y; q[2]=z;
180  }
181  else if (op == "c2y") {
182  q[0]=-x; q[1]=y; q[2]=-z;
183  }
184  else if (op == "c2x") {
185  q[0]=x; q[1]=-y; q[2]=-z;
186  }
187  else if (op == "sxy") {
188  q[0]=x; q[1]=y; q[2]=-z;
189  }
190  else if (op == "sxz") {
191  q[0]=x; q[1]=-y; q[2]=z;
192  }
193  else if (op == "syz") {
194  q[0]=-x; q[1]=y; q[2]=z;
195  }
196  else if (op == "i") {
197  q[0]=-x; q[1]=-y; q[2]=-z;
198  }
199  else {
200  throw "PointGroup: apply_op_by_name: unknown operator name";
201  }
202  return q;
203  }
204 
205  /// Returns the irrep of the Cartesian axis (0, 1, 2 = x, y, z)
206  int cart_ir(int axis) const {
207  coordT r(0.0);
208  r[axis] = 1.0;
209  for (int ir=0; ir<order; ++ir) {
210  double sum = 0.0;
211  for (int op=0; op<order; ++op) {
212  sum += apply(op,r)[axis] * c[ir][op];
213  }
214  sum /= order;
215  if (sum > 0.9) {
216  return ir;
217  }
218  }
219  throw "PointGroup: cart_ir: problem identifying axis";
220  }
221 
222  /// Returns irreducible cell
223 
224  /// The irreducible cell is defined by considering the cube
225  /// -1<x<1, -1<y<1, -1<z<1. Use symmetry operations to map
226  /// negative coordinates onto positive coordinates if possible.
227  /// This leads to the unique cell for D2h being the cube 0<x<1,
228  /// 0<y<1, 0<z<1.
229  ///
230  /// If a coordinate in the cell is positive, it means that the
231  /// irreducible cell has positive values of that coordinate.
232  coordT ircell() const {
233  double xmin=1.0, ymin=1.0, zmin=1.0;
234  // Loop thru corners of the cube
235  for (int x=-1; x<=1; x+=2) {
236  for (int y=-1; y<=1; y+=2) {
237  for (int z=-1; z<=1; z+=2) {
238  // Find the most positive corner it can be mapped into
239  double rx=x, ry=y, rz=z;
240  for (int op=0; op<order; ++op) {
241  coordT r;
242  r[0] = rx; r[1] = ry; r[2] = rz;
243  coordT q = apply(op,r);
244  double xx = q[0], yy=q[1], zz = q[2];
245  if ((xx>rx) || (xx==rx && yy>ry) || (xx==rx && yy==ry && zz>rz)) {
246  rx=xx; ry=yy; rz=zz;
247  }
248  }
249  xmin = std::min(xmin,rx); ymin = std::min(ymin,ry); zmin = std::min(zmin,rz);
250  }
251  }
252  }
253 
254  coordT r;
255  r[0] = xmin; r[1] = ymin; r[2] = zmin;
256  return r;
257  }
258 
259  static bool test() {
260  std::cout << PointGroup("C1");
261  std::cout << PointGroup("C2");
262  std::cout << PointGroup("Ci");
263  std::cout << PointGroup("Cs");
264  std::cout << PointGroup("C2h");
265  std::cout << PointGroup("C2v");
266  std::cout << PointGroup("D2");
267  std::cout << PointGroup("D2h");
268  return true;
269  }
270 
271  const std::string& get_name() const {
272  return name;
273  }
274 
275  int get_order() const {
276  return order;
277  }
278 
279  const std::string& get_ir_name(int ir) const {
280  return irs[ir];
281  }
282 
283  const std::string& get_op_name(int op) const {
284  return ops[op];
285  }
286 
287  int ctable(int ir, int op) const {
288  return c[ir][op];
289  }
290 };
291 
292 std::ostream& operator<<(std::ostream& s, const PointGroup& g) {
293  int order = g.get_order();
294  s << "\n";
295  s << "Group " << g.get_name() << " - irreducible cell " << g.ircell() << "\n";
296  s << "---------\n";
297  s << " ";
298  for (int op=0; op<order; ++op)
299  s << " " << std::setw(3) << g.get_op_name(op);
300  s << "\n";
301  s << " ";
302  for (int op=0; op<order; ++op)
303  s << " ---";
304 
305  s << "\n";
306 
307  int irx = g.cart_ir(0);
308  int iry = g.cart_ir(1);
309  int irz = g.cart_ir(2);
310 
311  for (int ir=0; ir<order; ++ir) {
312  s << " " << std::left << std::setw(3) << g.get_ir_name(ir) << std::right << " ";
313  for (int op=0; op<order; ++op) {
314  s << " " << std::setw(3) << g.ctable(ir, op);
315  }
316  if (ir == irx) s << " x";
317  if (ir == iry) s << " y";
318  if (ir == irz) s << " z";
319  s << "\n";
320  }
321  s << "\n";
322  return s;
323 }
324 
325 #endif
double q(double t)
Definition: DKops.h:18
Definition: pointgroup.h:50
int order
group order
Definition: pointgroup.h:52
coordT ircell() const
Returns irreducible cell.
Definition: pointgroup.h:232
const std::string & get_name() const
Definition: pointgroup.h:271
PointGroup(const std::string name)
Constructs point group by name (D2h and subgroups only)
Definition: pointgroup.h:61
const std::string & get_ir_name(int ir) const
Definition: pointgroup.h:279
int ctable(int ir, int op) const
Definition: pointgroup.h:287
std::string irs[8]
Names of the irreps.
Definition: pointgroup.h:53
int irmul(int ir1, int ir2) const
Returns irreducible representation corresponding to product.
Definition: pointgroup.h:162
virtual ~PointGroup()
Destructor.
Definition: pointgroup.h:159
static coordT apply(const std::string &op, const coordT &r)
Applies named operator (e, c2z, c2y, c2x, sxy, sxz, syz, i) to point.
Definition: pointgroup.h:172
int get_order() const
Definition: pointgroup.h:275
coordT apply(int op, const coordT &r) const
Applies group operator number op (0,1,...,order-1) to point.
Definition: pointgroup.h:167
std::string ops[8]
Names of the operators.
Definition: pointgroup.h:54
madness::Vector< double, 3 > coordT
Definition: pointgroup.h:58
std::string name
group name
Definition: pointgroup.h:51
int c[8][8]
Character table.
Definition: pointgroup.h:55
int cart_ir(int axis) const
Returns the irrep of the Cartesian axis (0, 1, 2 = x, y, z)
Definition: pointgroup.h:206
const std::string & get_op_name(int op) const
Definition: pointgroup.h:283
PointGroup & operator=(const PointGroup &other)
Assignment.
Definition: pointgroup.h:143
static bool test()
Definition: pointgroup.h:259
PointGroup(const PointGroup &other)
Copy constructor.
Definition: pointgroup.h:138
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
std::ostream & operator<<(std::ostream &s, const PointGroup &g)
Definition: pointgroup.h:292
Defines and implements most of Tensor.
AtomicInt sum
Definition: test_atomicint.cc:46
double g(const coord_1d &r)
Definition: testgconv.cc:49
std::size_t axis
Definition: testpdiff.cc:59
Implement the madness:Vector class, an extension of std::array that supports some mathematical operat...