Index: ThirdParty/mpqc_open/src/lib/math/symmetry/LIBS.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/LIBS.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/LIBS.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,5 @@
+libSCsymmetry.LIBSUF
+#include <math/scmat/LIBS.h>
+#include <util/misc/LIBS.h>
+#include <util/state/LIBS.h>
+#include <util/class/LIBS.h>
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/Makefile
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/Makefile	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/Makefile	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,63 @@
+#
+# Makefile
+#
+# Copyright (C) 1996 Limit Point Systems, Inc.
+#
+# Author: Edward Seidl <seidl@janed.com>
+# Maintainer: LPS
+#
+# This file is part of the SC Toolkit.
+#
+# The SC Toolkit is free software; you can redistribute it and/or modify
+# it under the terms of the GNU Library General Public License as published by
+# the Free Software Foundation; either version 2, or (at your option)
+# any later version.
+#
+# The SC Toolkit is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Library General Public License for more details.
+#
+# You should have received a copy of the GNU Library General Public License
+# along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+# the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+#
+# The U.S. Government is granted a limited license as per AL 91-7.
+#
+
+TOPDIR=../../../..
+ifndef SRCDIR
+  SRCDIR=$(shell pwd)
+endif
+
+include $(SRCDIR)/$(TOPDIR)/lib/GlobalMakefile
+
+CXXSRC = symop.cc rep.cc irrep.cc chartab.cc pointgrp.cc \
+         maketab.cc tetra.cc ico.cc corrtab.cc
+
+LIBOBJ = $(CXXSRC:%.cc=%.$(OBJSUF))
+
+INC =  pointgrp.h
+
+DEPENDINCLUDE = $(INC)
+
+BIN_OR_LIB = LIB
+TARGET_TO_MAKE = libSCsymmetry
+
+TESTPROGS = testpg
+
+DISTFILES = $(CXXSRC) $(INC) Makefile LIBS.h
+
+default:: $(DEPENDINCLUDE)
+
+include $(SRCDIR)/$(TOPDIR)/lib/GlobalRules
+
+LIBS := $(shell $(LISTLIBS) $(INCLUDE) $(SRCDIR)/LIBS.h)
+testpg: testpg.$(OBJSUF) $(LIBS)
+	$(LTLINK) $(CXX) $^ -o $@ $(SYSLIBS)$(LTLINKBINOPTS)
+
+$(LIBOBJ:.$(OBJSUF)=.d): $(DEPENDINCLUDE)
+ifneq ($(DODEPEND),no)
+include $(LIBOBJ:.$(OBJSUF)=.d)
+endif
+
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/chartab.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/chartab.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/chartab.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,350 @@
+//
+// chartab.cc
+//
+// Modifictions are
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+/* chartab.cc -- implementation of the point group classes
+ *
+ *      THIS SOFTWARE FITS THE DESCRIPTION IN THE U.S. COPYRIGHT ACT OF A
+ *      "UNITED STATES GOVERNMENT WORK".  IT WAS WRITTEN AS A PART OF THE
+ *      AUTHOR'S OFFICIAL DUTIES AS A GOVERNMENT EMPLOYEE.  THIS MEANS IT
+ *      CANNOT BE COPYRIGHTED.  THIS SOFTWARE IS FREELY AVAILABLE TO THE
+ *      PUBLIC FOR USE WITHOUT A COPYRIGHT NOTICE, AND THERE ARE NO
+ *      RESTRICTIONS ON ITS USE, NOW OR SUBSEQUENTLY.
+ *
+ *  Author:
+ *      E. T. Seidl
+ *      Bldg. 12A, Rm. 2033
+ *      Computer Systems Laboratory
+ *      Division of Computer Research and Technology
+ *      National Institutes of Health
+ *      Bethesda, Maryland 20892
+ *      Internet: seidl@alw.nih.gov
+ *      June, 1993
+ */
+
+#include <ctype.h>
+
+#include <util/misc/formio.h>
+#include <util/misc/newstring.h>
+#include <math/symmetry/pointgrp.h>
+
+using namespace std;
+using namespace sc;
+
+////////////////////////////////////////////////////////////////////////
+
+CharacterTable::CharacterTable()
+  : g(0), nt(0), pg(C1), nirrep_(0), gamma_(0), symop(0), _inv(0), symb(0)
+{
+}
+
+CharacterTable::CharacterTable(const CharacterTable& ct)
+  : g(0), nt(0), pg(C1), nirrep_(0), gamma_(0), symop(0), _inv(0), symb(0)
+{
+  *this = ct;
+}
+
+CharacterTable::~CharacterTable()
+{
+  if (symb) delete[] symb; symb=0;
+  if (gamma_) delete[] gamma_; gamma_=0;
+  if (symop) delete[] symop; symop=0;
+  if (_inv) delete[] _inv; _inv=0;
+  g=nt=nirrep_=0;
+}
+
+CharacterTable&
+CharacterTable::operator=(const CharacterTable& ct)
+{
+  g=ct.g; nt=ct.nt; pg=ct.pg; nirrep_=ct.nirrep_;
+  
+  if (symb)
+    delete[] symb;
+
+  symb = new_string(ct.symb);
+
+  if (gamma_) delete[] gamma_; gamma_=0;
+  if (ct.gamma_) {
+    gamma_ = new IrreducibleRepresentation[nirrep_];
+    for (int i=0; i < nirrep_; i++) {
+      gamma_[i].init();
+      gamma_[i] = ct.gamma_[i];
+    }
+  }
+
+  if (symop)
+    delete[] symop;
+  symop=0;
+
+  if (ct.symop) {
+    symop = new SymmetryOperation[g];
+    for (int i=0; i < g; i++) {
+      symop[i] = ct.symop[i];
+    }
+  }
+
+  if (_inv)
+    delete[] _inv;
+  _inv=0;
+
+  if (ct._inv) {
+    _inv = new int[g];
+    memcpy(_inv,ct._inv,sizeof(int)*g);
+  }
+
+  return *this;
+}
+
+void
+CharacterTable::print(ostream& os) const
+{
+  if (!g || !nirrep_) return;
+
+  int i;
+
+  os << indent << "point group " << symb << endl << endl;
+
+  for (i=0; i < nirrep_; i++)
+    gamma_[i].print(os);
+
+  os << endl << indent << "symmetry operation matrices:"
+     << endl << endl << incindent;
+  for (i=0; i < g; i++)
+    symop[i].print(os);
+
+  os << decindent << indent << "inverse symmetry operation matrices:"
+     << endl << endl << incindent;
+  for (i=0; i < g; i++)
+    symop[inverse(i)].print(os);
+
+  os << decindent;
+}
+
+CharacterTable::CharacterTable(const char *cpg, const SymmetryOperation& frame)
+  : g(0), nt(0), pg(C1), nirrep_(0), gamma_(0), symop(0), _inv(0), symb(0)
+{
+  // first parse the point group symbol, this will give us the order of the
+  // point group(g), the type of point group (pg), the order of the principle
+  // rotation axis (nt), and the number of irreps (nirrep_)
+
+  if (!cpg) {
+    ExEnv::errn() << "CharacterTable::CharacterTable: null point group" << endl;
+    exit(1);
+  }
+
+  symb = new char[strlen(cpg)+1];
+  size_t i;
+  for (i=0; i < strlen(cpg); i++) symb[i] = tolower(cpg[i]);
+  symb[i] = '\0';
+
+  if (parse_symbol() < 0) {
+    ExEnv::errn() << "CharacterTable::CharacterTable: invalid point group "
+         << cpg << endl;
+    exit(1);
+  }
+
+  if (make_table() < 0) {
+    ExEnv::errn() << "CharacterTable::CharacterTable: could not make table" << endl;
+    exit(1);
+  }
+
+  int ig;
+  for (ig=0; ig < g; ig++)
+    symop[ig] = symop[ig].transform(frame);
+}
+
+CharacterTable::CharacterTable(const char *cpg)
+  : g(0), nt(0), pg(C1), nirrep_(0), gamma_(0), symop(0), _inv(0), symb(0)
+{
+  // first parse the point group symbol, this will give us the order of the
+  // point group(g), the type of point group (pg), the order of the principle
+  // rotation axis (nt), and the number of irreps (nirrep_)
+
+  if (!cpg) {
+    ExEnv::errn() << "CharacterTable::CharacterTable: null point group" << endl;
+    exit(1);
+  }
+
+  symb = new char[strlen(cpg)+1];
+  size_t i;
+  for (i=0; i < strlen(cpg); i++) symb[i] = tolower(cpg[i]);
+  symb[i] = '\0';
+
+  if (parse_symbol() < 0) {
+    ExEnv::errn() << "CharacterTable::CharacterTable: invalid point group "
+         << cpg << endl;
+    exit(1);
+  }
+
+  if (make_table() < 0) {
+    ExEnv::errn() << "CharacterTable::CharacterTable: could not make table" << endl;
+    exit(1);
+  }
+}
+
+int
+CharacterTable::parse_symbol()
+{
+  // default to C1 symmetry
+  g=1; pg=C1; nt=1; nirrep_=1;
+
+  if (!symb) return 0;
+
+  if (!strcmp(symb,"c1")) return 0;
+
+  if (!strcmp(symb,"ci")) {
+    g = 2; pg = CI; nirrep_ = 2; nt = 2;
+    return 0;
+  }
+
+  if(!strcmp(symb,"cs")) {
+    g = 2; pg = CS; nirrep_ = 2; nt = 0;
+    return 0;
+  }
+
+  if (symb[0] == 'c') {
+    int nab,ne;
+
+    if (symb[1] == '\0') return -1;
+
+    nt = atoi(&symb[1]);
+    ne = (nt%2) ? nt/2 : nt/2-1;
+    nab = (nt%2) ? 1 : 2;
+
+    char *vhd = &symb[1];
+    while (*vhd && isdigit(*vhd))
+      vhd++;
+    
+    if (*vhd) {
+      if (*vhd == 'v') {
+        g  = 2*nt; pg = CNV; nirrep_ = 2*nab + ne;
+      } else if (*vhd == 'h') {
+        g  = 2*nt; pg = CNH; nirrep_ = 2*(nab+ne);
+      } else {
+        return -1;
+      }
+    } else {
+      g = nt; pg = CN; nirrep_ = nab+ne;
+    }
+
+    return 0;
+  }
+
+  if (symb[0] == 'd') {
+    int nab,ne;
+
+    if (symb[1] == '\0') return -1;
+
+    nt = atoi(&symb[1]);
+    ne = (nt%2) ? nt/2 : nt/2-1;
+    nab = (nt%2) ? 1 : 2;
+
+    char *vhd = &symb[1];
+    while (*vhd && isdigit(*vhd))
+      vhd++;
+    
+    if (*vhd) {
+      if (*vhd == 'd') {
+        g = 4*nt; pg = DND; nirrep_ = nt+3;
+      } else if (*vhd == 'h') {
+        g = 4*nt; pg = DNH; nirrep_ = 4*nab + 2*ne;
+      } else {
+        return -1;
+      }
+    } else {
+      g = 2*nt; pg = DN; nirrep_ = 2*nab + ne;
+    }
+
+    return 0;
+  }
+
+  if (symb[0] == 's') {
+    if (symb[1] == '\0') return -1;
+
+    nt = atoi(&symb[1]);
+
+    // only S2n groups make sense
+    if (nt%2) return -1;
+
+    g = nt; pg = SN; nirrep_ = nt/2+1;
+
+    return 0;
+  }
+
+  if (symb[0] == 't') {
+    if (symb[1] != '\0') {
+      if (symb[1] == 'd') {
+        g = 24; pg = TD; nirrep_ = 5;
+      } else if(symb[1] == 'h') {
+        g = 24; pg = TH; nirrep_ = 6;
+      } else {
+        return -1;
+      }
+    } else {
+      g = 12; pg = T; nirrep_ = 3;
+    }
+
+    return 0;
+  }
+
+  if (symb[0] == 'o') {
+    if (symb[1] != '\0') {
+      if (symb[1] == 'h') {
+        pg = OH; g = 48; nirrep_ = 10;
+      } else {
+        return -1;
+      }
+    } else {
+      g = 24; pg = O; nirrep_ = 5;
+    }
+
+    return 0;
+  }
+
+  if (symb[0] == 'i') {
+    if (symb[1] != '\0') {
+      if (symb[1] == 'h') {
+        g = 120; pg = IH; nirrep_ = 10;
+      } else {
+        return -1;
+      }
+    } else {
+      g = 60; pg = I; nirrep_ = 5;
+    }
+
+    return 0;
+  }
+
+  return -1;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/corrtab.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/corrtab.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/corrtab.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,248 @@
+//
+// pointgrp.cc
+//
+// Copyright (C) 1997 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <janssen@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <util/misc/formio.h>
+#include <math/symmetry/corrtab.h>
+
+using namespace std;
+using namespace sc;
+
+////////////////////////////////////////////////////////////////////////
+
+CorrelationTable::CorrelationTable():
+  n_(0),
+  ngamma_(0),
+  gamma_(0)
+{
+}
+
+CorrelationTable::CorrelationTable(const Ref<PointGroup>& group,
+                                   const Ref<PointGroup>& subgroup):
+  n_(0),
+  ngamma_(0),
+  gamma_(0)
+{
+  int rc = initialize_table(group,subgroup);
+  if (rc != 0) {
+    ExEnv::err0()
+         << "ERROR: CorrelationTable: " << error(rc) << endl;
+    abort();
+    }
+}
+
+CorrelationTable::~CorrelationTable()
+{
+  clear();
+}
+
+int
+CorrelationTable::initialize_table(const Ref<PointGroup>& group,
+                                   const Ref<PointGroup>& subgroup)
+{
+  clear();
+
+  group_ = group;
+  subgroup_ = subgroup;
+
+  int i, j, k, l;
+
+  CharacterTable ct = group_->char_table();
+  CharacterTable subct = subgroup_->char_table();
+
+  n_ = ct.nirrep();
+  subn_ = subct.nirrep();
+  ngamma_ = new int[n_];
+  gamma_ = new int*[n_];
+
+  // CAN ONLY HANDLE NONDEGENERATE POINT GROUPS
+
+  for (i=0; i<n_; i++) {
+    ngamma_[i] = 0;
+    gamma_[i] = 0;
+    }
+
+  // map the ops in the high order to low order groups
+  int *so_to_subso = new int[ct.order()];
+  int *subso_to_so = new int[subct.order()];
+  for (i=0; i<subct.order(); i++) subso_to_so[i] = -1;
+  for (i=0; i<ct.order(); i++) {
+    SymmetryOperation so = ct.symm_operation(i);
+    int found = 0;
+    so_to_subso[i] = -1;
+    for (j=0; j<subct.order(); j++) {
+      SymmetryOperation subso = subct.symm_operation(j);
+      double sumsquare = 0.0;
+      for (k=0; k<3; k++) {
+        for (l=0; l<3; l++) {
+          double diff = so(k,l)-subso(k,l);
+          sumsquare += diff*diff;
+        }
+      }
+      if (sumsquare < 1.0e-12) {
+        found++;
+        //ExEnv::outn() << scprintf("symmop %d in %s is %d in %s",
+        //                 i,ct.symbol(),j,subct.symbol()) << endl;
+        so_to_subso[i] = j;
+        subso_to_so[j] = i;
+        }
+      }
+    if (found > 1) {
+      delete[] so_to_subso;
+      delete[] subso_to_so;
+      return -1;
+      }
+    }
+  for (i=0; i<subct.order(); i++) {
+    if (subso_to_so[i] == -1) {
+      delete[] so_to_subso;
+      delete[] subso_to_so;
+      return -2;
+      }
+    }
+
+  // compute the correlation table
+  for (i=0; i<ct.nirrep(); i++) {
+    for (j=0; j<subct.nirrep(); j++) {
+      double nmatch = 0.0;
+      for (k=0; k<ct.order(); k++) {
+        double chr = ct.gamma(i).character(k);
+        if (so_to_subso[k] >= 0) {
+          double subchr = subct.gamma(j).character(so_to_subso[k]);
+          nmatch += subchr*chr;
+          }
+        }
+      nmatch /= subct.order();
+      int inmatch = (int)(nmatch+0.5);
+      if (fabs(nmatch-inmatch)>1.0e-6) {
+        delete[] so_to_subso;
+        delete[] subso_to_so;
+        return -4;
+        }
+      if (inmatch > 0) {
+        int *newgamma = new int[ngamma_[i] + inmatch];
+        memcpy(newgamma,gamma_[i],ngamma_[i]*sizeof(int));
+        for (k=0; k<inmatch; k++) newgamma[ngamma_[i]+k] = j;
+        ngamma_[i] += inmatch;
+        delete[] gamma_[i];
+        gamma_[i] = newgamma;
+        }
+      }
+    }
+
+  delete[] so_to_subso;
+  delete[] subso_to_so;
+
+  for (i=0; i<n(); i++) {
+    int degen = ct.gamma(i).degeneracy();
+    int subdegen = 0;
+    for (j=0; j<ngamma(i); j++) {
+      subdegen += subct.gamma(gamma(i,j)).degeneracy();
+      }
+    if (degen != subdegen) {
+      return -3;
+      }
+    }
+
+  return 0;
+}
+
+void
+CorrelationTable::clear()
+{
+  for (int i=0; i<n_; i++) {
+    delete[] gamma_[i];
+    }
+  delete[] ngamma_;
+  delete[] gamma_;
+}
+
+const char *
+CorrelationTable::error(int rc)
+{
+  if (rc == -1) {
+    return "too many symop matches";
+    }
+  else if (rc == -2) {
+    return "not a subgroup or wrong ref frame";
+    }
+  else if (rc == -3) {
+    return "only groups with non-complex characters supported (degen sum)";
+    }
+  else if (rc == -4) {
+    return "only groups with non-complex characters supported (nonint nirr)";
+    }
+  else if (rc != 0) {
+    return "unknown error";
+    }
+  return "no error";
+}
+
+int
+CorrelationTable::degen(int i) const
+{
+  return group_->char_table().gamma(i).degeneracy();
+}
+
+int
+CorrelationTable::subdegen(int i) const
+{
+  return subgroup_->char_table().gamma(i).degeneracy();
+}
+
+void
+CorrelationTable::print(ostream &o) const
+{
+  o << indent
+    << "Correlation Table from "
+    << group_->symbol() << " to " << subgroup_->symbol() << ":" << endl;
+
+  CharacterTable ct = group_->char_table();
+  CharacterTable subct = subgroup_->char_table();
+
+  o << incindent;
+  for (int i=0; i<n(); i++) {
+    o << indent
+      << ct.gamma(i).symbol() << ":";
+    for (int j=0; j<ngamma(i); j++) {
+      o << indent
+        << " " << subct.gamma(gamma(i,j)).symbol();
+      }
+    o << endl;
+    }
+  o << decindent;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ-CONDENSED"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/corrtab.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/corrtab.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/corrtab.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,105 @@
+//
+// corrtab.h
+//
+// Copyright (C) 1997 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#ifndef _math_symmetry_corrtab_h
+#define _math_symmetry_corrtab_h
+
+#include <iostream>
+
+#include <math/symmetry/pointgrp.h>
+
+namespace sc {
+
+// //////////////////////////////////////////////////////////////////
+ 
+/** The CorrelationTable class provides a correlation
+    table between two point groups.
+*/
+class CorrelationTable: public RefCount {
+  private:
+    Ref<PointGroup> group_;
+    Ref<PointGroup> subgroup_;
+
+    int n_;
+    int subn_;
+    int *ngamma_;
+    int **gamma_;
+
+    void clear();
+  public:
+    CorrelationTable();
+
+    /// Create a correlation table for the two groups.
+    CorrelationTable(const Ref<PointGroup>& group,
+                     const Ref<PointGroup>& subgroup);
+
+    ~CorrelationTable();
+
+    /// Returns the higher order point group.
+    Ref<PointGroup> group() const { return group_; }
+    /// Returns the lower order point group.
+    Ref<PointGroup> subgroup() const { return subgroup_; }
+
+    /** Initalize the correlation table.  Returns 0 for success and nonzero
+        for failure.  This will fail if the subgroup is not really a subgroup
+        of group. */
+    int initialize_table(const Ref<PointGroup>& group,
+                         const Ref<PointGroup>& subgroup);
+
+    /// Converts error codes from initialize_table into a text string.
+    const char *error(int errcod);
+
+    /// Returns the number of irreps in the high order group.
+    int n() const { return n_; }
+    /// Returns the number of irreps in the subgroup.
+    int subn() const { return subn_; }
+    /// Returns the degeneracy of the irrep.
+    int degen(int igamma) const;
+    /// Returns the degeneracy of the subgroup irrep.
+    int subdegen(int igamma) const;
+    /// Returns the number of irreps in the low order group that an irrep
+    //from the high order group can be reduced to.
+    int ngamma(int igamma) const { return ngamma_[igamma]; }
+    /** Returns the irreps in the low order group that an irrep from the
+        high order group can be reduced to. */
+    int gamma(int igamma, int i) const { return gamma_[igamma][i]; }
+
+    void print(std::ostream &o=ExEnv::out0()) const;
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ-CONDENSED"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/ico.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/ico.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/ico.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,397 @@
+//
+// ico.cc --- implementation of icosahedral operations
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#include <util/misc/math.h>
+#include <string.h>
+
+#include <math/symmetry/pointgrp.h>
+
+using namespace sc;
+
+// these are the operations which make up T
+static void
+i_ops(SymRep *t1rep, SymRep *t2rep, SymRep *grep, SymRep *hrep)
+{
+  int i;
+  
+  // identity
+  t1rep[0].E();
+  t2rep[0].E();
+  grep[0].E();
+  hrep[0].E();
+    
+  //
+  // 12 C5's
+  //
+  // first the 2 C5's about the z axis
+  t1rep[1].rotation(2.0*(double)M_PI/5.0);
+  t1rep[2].rotation(8.0*(double)M_PI/5.0);
+  
+  t2rep[1] = t1rep[1].operate(t1rep[1]);
+  t2rep[2] = t1rep[2].operate(t1rep[2]);
+
+  grep[1].rotation(2.0*(double)M_PI/5.0);
+  grep[2].rotation(8.0*(double)M_PI/5.0);
+  
+  hrep[1].rotation(2.0*(double)M_PI/5.0);
+  hrep[2].rotation(8.0*(double)M_PI/5.0);
+   
+  // form rotation matrices for the C3 axis about the zx axis (these were
+  // taken from turbomole version 2, which claims they were sort of inherited
+  // from hondo
+  SymRep t1so(3);
+  SymRep gso(4);
+  SymRep hso(5);
+
+  double c2p5 = cos(2.0*(double)M_PI/5.0);
+  double s2p5 = sin(2.0*(double)M_PI/5.0);
+  double cosd = s2p5/((1.0-c2p5)*sqrt(3.0));
+  double cosd2 = cosd*cosd;
+  double sind2 = 1.0 - cosd2;
+  double sind = sqrt(sind2);
+
+  t1so[0][0] =  1.0 - 1.5*cosd2;
+  t1so[1][0] =  0.5*sqrt(3.0)*cosd;
+  t1so[2][0] =  1.5*cosd*sind;
+  t1so[0][1] = -0.5*sqrt(3.0)*cosd;
+  t1so[1][1] = -0.5;
+  t1so[2][1] =  0.5*sqrt(3.0)*sind;
+  t1so[0][2] =  1.5*cosd*sind;
+  t1so[1][2] = -0.5*sqrt(3.0)*sind;
+  t1so[2][2] =  1.0 - 1.5*sind2;
+
+  gso[0][0] = (3.0*sqrt(5.0)+5.0)/20.0;
+  gso[0][1] = cosd*sqrt(3.0)*(sqrt(5.0)-1.0)/4.0;
+  gso[0][2] = 3.0*sqrt(5.0)/10.0;
+  gso[0][3] = -sqrt(5.0-2.0*sqrt(5.0))*sqrt(5.0)/10.0;
+  gso[1][0] = -gso[0][1];
+  gso[1][1] = (1-sqrt(5.0))/4.0;
+  gso[1][2] = cosd*sqrt(3.0)/2.0;
+  gso[1][3] = cosd*sqrt(5-2*sqrt(5.0))*sqrt(3.0)/2.0;
+  gso[2][0] = gso[0][2];
+  gso[2][1] = -gso[1][2];
+  gso[2][2] = (5-3*sqrt(5.0))/20.0;
+  gso[2][3] = sqrt(5.0-2*sqrt(5.0))*(sqrt(5.0)+5)/20;
+  gso[3][0] = -gso[0][3];
+  gso[3][1] = gso[1][3];
+  gso[3][2] = -gso[2][3];
+  gso[3][3] = (sqrt(5.0)+1)/4.0;
+
+  hso[0][0] = -1.0/5.0;
+  hso[0][4] = sqrt(3.0)*(sqrt(5.0)+1)/10.0;
+  hso[0][3] = 3.0*cosd*(3.0*sqrt(5.0)-5.0)/10.0;
+  hso[0][2] = 3.0*cosd*(5.0-sqrt(5.0))/10.0;
+  hso[0][1] = sqrt(3.0)*(sqrt(5.0)-1.0)/10.0;
+  hso[4][0] = hso[0][4];
+  hso[4][4] = (2.0*sqrt(5.0)+1.0)/10.0;
+  hso[4][3] = sqrt(3.0)*cosd*(5.0-2.0*sqrt(5.0))/10.0;
+  hso[4][2] = sqrt(3.0)*cosd*(5.0-3.0*sqrt(5.0))/5.0;
+  hso[4][1] = 2.0/5.0;
+  hso[3][0] = -hso[0][3];
+  hso[3][4] = -hso[4][3];
+  hso[3][3] = -1.0/2.0;
+  hso[3][2] = 0.0;
+  hso[3][1] = sqrt(3.0)*cosd*(5.0-sqrt(5.0))/5.0;
+  hso[2][0] = -hso[0][2];
+  hso[2][4] = -hso[4][2];
+  hso[2][3] = 0.0;
+  hso[2][2] = -1.0/2.0;
+  hso[2][1] = -sqrt(3.0)*sqrt(5.0)*cosd/10.0;
+  hso[1][0] = hso[0][1];
+  hso[1][4] = hso[4][1];
+  hso[1][3] = -hso[3][1];
+  hso[1][2] = -hso[2][1];
+  hso[1][1] = (1.0-2.0*sqrt(5.0))/10.0;
+  
+  // now rotate the first C5's by 2pi/3 degrees about the zx axis (sort of)
+  t1rep[3] = t1rep[1].transform(t1so);
+  t1rep[4] = t1rep[2].transform(t1so);
+
+  grep[3] = grep[1].transform(gso);
+  grep[4] = grep[2].transform(gso);
+
+  hrep[3] = hrep[1].transform(hso);
+  hrep[4] = hrep[2].transform(hso);
+
+  // rotate twice to get the first one aligned along the x axis
+  t1rep[3] = t1rep[3].transform(t1rep[1]).transform(t1rep[1]);
+  t1rep[4] = t1rep[4].transform(t1rep[1]).transform(t1rep[1]);
+
+  grep[3] = grep[3].transform(grep[1]).transform(grep[1]);
+  grep[4] = grep[4].transform(grep[1]).transform(grep[1]);
+
+  hrep[3] = hrep[3].transform(hrep[1]).transform(hrep[1]);
+  hrep[4] = hrep[4].transform(hrep[1]).transform(hrep[1]);
+
+  t2rep[3] = t1rep[4].operate(t1rep[4]);
+  t2rep[4] = t1rep[3].operate(t1rep[3]);
+
+  t2rep[13] = t1rep[2];
+  t2rep[14] = t1rep[1];
+
+  t2rep[15] = t1rep[3];
+  t2rep[16] = t1rep[4];
+  
+  // and then rotate those by 2pi/5 about the z axis 4 times
+  for (i=5; i < 13; i++) {
+    t1rep[i] = t1rep[i-2].transform(t1rep[1]);
+    grep[i] = grep[i-2].transform(grep[1]);
+    hrep[i] = hrep[i-2].transform(hrep[1]);
+
+    t2rep[i] = t2rep[i-2].transform(t2rep[1]);
+    t2rep[i+12] = t2rep[i+10].transform(t2rep[1]);
+  }
+
+  //
+  // 12 C5^2's
+  //
+  // get these from operating on each of the C5's with itself
+  for (i=13; i < 25; i++) {
+    t1rep[i] = t1rep[i-12].operate(t1rep[i-12]);
+    grep[i] = grep[i-12].operate(grep[i-12]);
+    hrep[i] = hrep[i-12].operate(hrep[i-12]);
+  }
+
+  //
+  // 20 C3's
+  //
+  // first we have 2 C3's about the zx axis
+  t1rep[25] = t1so;
+  t1rep[26] = t1so.operate(t1so);
+  
+  grep[25] = gso;
+  grep[26] = gso.operate(gso);
+  
+  hrep[25] = hso;
+  hrep[26] = hso.operate(hso);
+  
+  // and then rotate those by 2pi/5 about the z axis 4 times
+  for (i=27; i < 35; i++) {
+    t1rep[i] = t1rep[i-2].transform(t1rep[1]);
+    grep[i] = grep[i-2].transform(grep[1]);
+    hrep[i] = hrep[i-2].transform(hrep[1]);
+  }
+
+  // now rotate one of the above C3's by 2pi/3 about the zx axis
+  t1rep[35] = t1rep[27].transform(t1so);
+  t1rep[36] = t1rep[28].transform(t1so);
+
+  grep[35] = grep[27].transform(gso);
+  grep[36] = grep[28].transform(gso);
+                      
+  hrep[35] = hrep[27].transform(hso);
+  hrep[36] = hrep[28].transform(hso);
+
+  // and then rotate those by 2pi/5 about the z axis 4 times
+  for (i=37; i < 45; i++) {
+    t1rep[i] = t1rep[i-2].transform(t1rep[1]);
+    grep[i] = grep[i-2].transform(grep[1]);
+    hrep[i] = hrep[i-2].transform(hrep[1]);
+  }
+
+  t2rep[25] = t1rep[35];
+  t2rep[26] = t1rep[36];
+  
+  for (i=27; i < 35; i++)
+    t2rep[i] = t2rep[i-2].transform(t2rep[1]);
+  
+  t2rep[35] = t1rep[26];
+  t2rep[36] = t1rep[25];
+  
+  for (i=37; i < 45; i++)
+    t2rep[i] = t2rep[i-2].transform(t2rep[1]);
+
+  //
+  // 15 C2's
+  //
+  // first we have a C2 about the y axis
+  t1rep[45][0][0] = -1.0;
+  t1rep[45][1][1] =  1.0;
+  t1rep[45][2][2] = -1.0;
+
+  t2rep[45] = t1rep[45];
+  
+  grep[45][0][0] = -1.0;
+  grep[45][1][1] =  1.0;
+  grep[45][2][2] = -1.0;
+  grep[45][3][3] =  1.0;
+  
+  hrep[45][0][0] =  1.0;
+  hrep[45][1][1] =  1.0;
+  hrep[45][2][2] = -1.0;
+  hrep[45][3][3] = -1.0;
+  hrep[45][4][4] =  1.0;
+  
+  // and rotate that by 2pi/5 about the z axis 4 times
+  for (i=46; i < 50; i++) {
+    t1rep[i] = t1rep[i-1].transform(t1rep[1]);
+    t2rep[i] = t2rep[i-1].transform(t2rep[1]);
+    grep[i] = grep[i-1].transform(grep[1]);
+    hrep[i] = hrep[i-1].transform(hrep[1]);
+  }
+
+  // now take the C2 about the y axis and rotate it by 2pi/3 about the zx axis
+  t1rep[50] = t1rep[45].transform(t1so);
+  grep[50] = grep[45].transform(gso);
+  hrep[50] = hrep[45].transform(hso);
+
+  // align this c2 along the x axis
+  t1rep[50] = t1rep[50].transform(t1rep[2]).transform(t1rep[2]);
+  grep[50] = grep[50].transform(grep[2]).transform(grep[2]);
+  hrep[50] = hrep[50].transform(hrep[2]).transform(hrep[2]);
+
+  // and rotate that by 2pi/5 about the z axis 4 times
+  for (i=51; i < 55; i++) {
+    t1rep[i] = t1rep[i-1].transform(t1rep[1]);
+    grep[i] = grep[i-1].transform(grep[1]);
+    hrep[i] = hrep[i-1].transform(hrep[1]);
+  }
+
+  // finally, take a C2 about the y axis, and rotate it by 2pi/3 about the
+  // xz axis, and align it along the x axis
+  t1rep[55] = t1rep[45].transform(t1rep[35]).transform(t1rep[1]);
+  grep[55] = grep[45].transform(grep[35]).transform(grep[1]);
+  hrep[55] = hrep[45].transform(hrep[35]).transform(hrep[1]);
+
+  // and then rotate that by 2pi/5 about the z axis 4 times
+  for (i=56; i < 60; i++) {
+    t1rep[i] = t1rep[i-1].transform(t1rep[1]);
+    grep[i] = grep[i-1].transform(grep[1]);
+    hrep[i] = hrep[i-1].transform(hrep[1]);
+  }
+
+  t2rep[50] = t1rep[55];
+  t2rep[55] = t1rep[50];
+  
+  for (i=51; i < 55; i++) {
+    t2rep[i] = t2rep[i-1].transform(t2rep[1]);
+    t2rep[i+5] = t2rep[i+4].transform(t2rep[1]);
+  }
+}
+
+void
+CharacterTable::i()
+{
+  int i;
+
+  IrreducibleRepresentation& ira = gamma_[0];
+  IrreducibleRepresentation& ir1 = gamma_[1];
+  IrreducibleRepresentation& ir2 = gamma_[2];
+  IrreducibleRepresentation& irg = gamma_[3];
+  IrreducibleRepresentation& irh = gamma_[4];
+
+  ira.init(g,1,"A");
+  ir1.init(g,3,"T1");
+  ir2.init(g,3,"T2");
+  irg.init(g,4,"G");
+  irh.init(g,5,"H");
+
+  // i_ops gives us all the symmetry operations we need
+  i_ops(ir1.rep, ir2.rep, irg.rep, irh.rep);
+    
+  ir1.nrot_ = 1;
+  ir1.ntrans_ = 1;
+
+  for (i=0; i < g; i++) {
+    ira.rep[i][0][0] = 1.0;
+    symop[i] = ir1.rep[i];
+  }
+}
+
+
+void CharacterTable::ih()
+{
+  int i;
+
+  IrreducibleRepresentation& irag = gamma_[0];
+  IrreducibleRepresentation& ir1g = gamma_[1];
+  IrreducibleRepresentation& ir2g = gamma_[2];
+  IrreducibleRepresentation& irgg = gamma_[3];
+  IrreducibleRepresentation& irhg = gamma_[4];
+
+  IrreducibleRepresentation& irau = gamma_[5];
+  IrreducibleRepresentation& ir1u = gamma_[6];
+  IrreducibleRepresentation& ir2u = gamma_[7];
+  IrreducibleRepresentation& irgu = gamma_[8];
+  IrreducibleRepresentation& irhu = gamma_[9];
+
+  irag.init(g,1,"Ag");
+  ir1g.init(g,3,"T1g");
+  ir2g.init(g,3,"T2g");
+  irgg.init(g,4,"Gg");
+  irhg.init(g,5,"Hg");
+
+  irau.init(g,1,"Au");
+  ir1u.init(g,3,"T1u");
+  ir2u.init(g,3,"T2u");
+  irgu.init(g,4,"Gu");
+  irhu.init(g,5,"Hu");
+
+  // i_ops gives us all the symmetry operations we need
+  i_ops(ir1g.rep, ir2g.rep, irgg.rep, irhg.rep);
+    
+  ir1g.nrot_ = 1;
+  ir1u.ntrans_ = 1;
+
+  SymRep ti(3), gi(4), hi(5);
+  ti.i();
+  gi.i();
+  hi.i();
+  
+  for (i=0; i < g/2; i++) {
+    irag.rep[i][0][0] = 1.0;
+    irau.rep[i][0][0] = 1.0;
+
+    irag.rep[i+60][0][0] =  1.0;
+    irau.rep[i+60][0][0] = -1.0;
+
+    ir1g.rep[i+60] = ir1g.rep[i];
+    ir2g.rep[i+60] = ir2g.rep[i];
+    irgg.rep[i+60] = irgg.rep[i];
+    irhg.rep[i+60] = irhg.rep[i];
+    
+    ir1u.rep[i] = ir1g.rep[i];
+    ir2u.rep[i] = ir2g.rep[i];
+    irgu.rep[i] = irgg.rep[i];
+    irhu.rep[i] = irhg.rep[i];
+    
+    ir1u.rep[i+60] = ir1g.rep[i].operate(ti);
+    ir2u.rep[i+60] = ir2g.rep[i].operate(ti);
+    irgu.rep[i+60] = irgg.rep[i].operate(gi);
+    irhu.rep[i+60] = irhg.rep[i].operate(hi);
+    
+    symop[i] = ir1u.rep[i];
+    symop[i+60] = ir1u.rep[i+60];
+  }
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/irrep.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/irrep.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/irrep.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,155 @@
+//
+// irrep.cc
+//
+// Modifications are
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+/* irrep.cc -- implementation of the point group classes
+ *
+ *      THIS SOFTWARE FITS THE DESCRIPTION IN THE U.S. COPYRIGHT ACT OF A
+ *      "UNITED STATES GOVERNMENT WORK".  IT WAS WRITTEN AS A PART OF THE
+ *      AUTHOR'S OFFICIAL DUTIES AS A GOVERNMENT EMPLOYEE.  THIS MEANS IT
+ *      CANNOT BE COPYRIGHTED.  THIS SOFTWARE IS FREELY AVAILABLE TO THE
+ *      PUBLIC FOR USE WITHOUT A COPYRIGHT NOTICE, AND THERE ARE NO
+ *      RESTRICTIONS ON ITS USE, NOW OR SUBSEQUENTLY.
+ *
+ *  Author:
+ *      E. T. Seidl
+ *      Bldg. 12A, Rm. 2033
+ *      Computer Systems Laboratory
+ *      Division of Computer Research and Technology
+ *      National Institutes of Health
+ *      Bethesda, Maryland 20892
+ *      Internet: seidl@alw.nih.gov
+ *      June, 1993
+ */
+
+
+#include <stdlib.h>
+
+#include <util/misc/newstring.h>
+#include <math/symmetry/pointgrp.h>
+#include <util/misc/formio.h>
+
+using namespace std;
+using namespace sc;
+
+/////////////////////////////////////////////////////////////////////////
+
+IrreducibleRepresentation::IrreducibleRepresentation() :
+  g(0), degen(0), nrot_(0), ntrans_(0), complex_(0), symb(0), rep(0), csymb(0)
+{
+}
+
+IrreducibleRepresentation::IrreducibleRepresentation(
+  int order, int d, const char *lab, const char *clab) :
+  g(0), degen(0), nrot_(0), ntrans_(0), complex_(0), symb(0), rep(0), csymb(0)
+{
+  init(order,d,lab,clab);
+}
+
+
+IrreducibleRepresentation::IrreducibleRepresentation(
+  const IrreducibleRepresentation& ir) :
+  g(0), degen(0), nrot_(0), ntrans_(0), complex_(0), symb(0), rep(0), csymb(0)
+{
+  *this = ir;
+}
+
+IrreducibleRepresentation::~IrreducibleRepresentation()
+{
+  init();
+}
+
+IrreducibleRepresentation&
+IrreducibleRepresentation::operator=(const IrreducibleRepresentation& ir)
+{
+  init(ir.g,ir.degen,ir.symb,ir.csymb);
+
+  nrot_ = ir.nrot_;
+  ntrans_ = ir.ntrans_;
+  complex_ = ir.complex_;
+  
+  for (int i=0; i < g; i++)
+    rep[i]=ir.rep[i];
+  
+  return *this;
+}
+
+void
+IrreducibleRepresentation::init(int order, int d, const char *lab,
+                                const char *clab)
+{
+  g=order;
+  degen=d;
+  ntrans_=nrot_=complex_=0;
+
+  delete[] symb;
+  symb = new_string(lab);
+
+  delete[] csymb;
+  if (clab) csymb = new_string(clab);
+  else csymb = 0;
+
+  if (rep) {
+    delete[] rep;
+    rep=0;
+  }
+
+  if (g) {
+    rep = new SymRep[g];
+    for (int i=0; i < g; i++)
+      rep[i].set_dim(d);
+  }
+}
+
+void
+IrreducibleRepresentation::print(ostream& os) const
+{
+  if (!g)
+    return;
+
+  int i,d;
+  
+  os << indent << scprintf("%-5s",symb);
+
+  for (i=0; i < g; i++)
+    os << scprintf(" %6.3f",character(i));
+  os << " | " << ntrans_ << " t, " << nrot_ << " R\n";
+
+  for (d=0; d < nproj(); d++) {
+    os << indent << "     ";
+    for (i=0; i < g; i++)
+      os << scprintf(" %6.3f",p(d,i));
+    os << endl;
+  }
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/maketab.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/maketab.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/maketab.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,1069 @@
+//
+// maketab.cc
+//
+// Modifications are
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+/* maketab.cc
+ *
+ *      THIS SOFTWARE FITS THE DESCRIPTION IN THE U.S. COPYRIGHT ACT OF A
+ *      "UNITED STATES GOVERNMENT WORK".  IT WAS WRITTEN AS A PART OF THE
+ *      AUTHOR'S OFFICIAL DUTIES AS A GOVERNMENT EMPLOYEE.  THIS MEANS IT
+ *      CANNOT BE COPYRIGHTED.  THIS SOFTWARE IS FREELY AVAILABLE TO THE
+ *      PUBLIC FOR USE WITHOUT A COPYRIGHT NOTICE, AND THERE ARE NO
+ *      RESTRICTIONS ON ITS USE, NOW OR SUBSEQUENTLY.
+ *
+ *  Author:
+ *      E. T. Seidl
+ *      Bldg. 12A, Rm. 2033
+ *      Computer Systems Laboratory
+ *      Division of Computer Research and Technology
+ *      National Institutes of Health
+ *      Bethesda, Maryland 20892
+ *      Internet: seidl@alw.nih.gov
+ *      June, 1993
+ */
+
+#include <util/misc/math.h>
+#include <stdio.h>
+#include <string.h>
+
+#include <math/symmetry/pointgrp.h>
+#include <util/misc/formio.h>
+
+using namespace std;
+using namespace sc;
+
+/*
+ * This function will generate a character table for the point group.
+ * This character table is in the order that symmetry operations are
+ * generated, not in Cotton order. If this is a problem, tough.
+ * Also generate the transformation matrices.
+ */
+
+int CharacterTable::make_table()
+{
+  int i,j,ei,gi;
+  char label[4];
+
+  if (!g) return 0;
+
+  gamma_ = new IrreducibleRepresentation[nirrep_];
+
+  symop = new SymmetryOperation[g];
+  SymmetryOperation so;
+
+  _inv = new int[g];
+  
+  // this array forms a reducible representation for rotations about x,y,z
+  double *rot = new double[g];
+  memset(rot,0,sizeof(double)*g);
+
+  // this array forms a reducible representation for translations along x,y,z
+  double *trans = new double[g];
+  memset(trans,0,sizeof(double)*g);
+
+  // the angle to rotate about the principal axis
+  double theta = (nt) ? 2.0*M_PI/nt : 2.0*M_PI;
+
+  switch (pg) {
+
+  case C1:
+    // no symmetry case
+    gamma_[0].init(1,1,"A");
+    gamma_[0].nrot_ = 3;
+    gamma_[0].ntrans_ = 3;
+    gamma_[0].rep[0][0][0] = 1.0;
+
+    symop[0].E();
+
+    break;
+
+  case CI:
+    // equivalent to S2 about the z axis
+    gamma_[0].init(2,1,"Ag");
+    gamma_[0].rep[0][0][0] = 1.0;
+    gamma_[0].rep[1][0][0] = 1.0;
+    gamma_[0].nrot_=3;
+
+    gamma_[1].init(2,1,"Au");
+    gamma_[1].rep[0][0][0] =  1.0;
+    gamma_[1].rep[1][0][0] = -1.0;
+    gamma_[1].ntrans_=3;
+
+    symop[0].E();
+    symop[1].i();
+
+    break;
+
+  case CS: // reflection through the xy plane
+    gamma_[0].init(2,1,"A'","Ap");
+    gamma_[0].rep[0][0][0] = 1.0;
+    gamma_[0].rep[1][0][0] = 1.0;
+    gamma_[0].nrot_=1;
+    gamma_[0].ntrans_=2;
+
+    gamma_[1].init(2,1,"A\"","App");
+    gamma_[1].rep[0][0][0] =  1.0;
+    gamma_[1].rep[1][0][0] = -1.0;
+    gamma_[1].nrot_=2;
+    gamma_[1].ntrans_=1;
+
+    symop[0].E();
+    symop[1].sigma_h();
+
+    break;
+
+  case CN:
+    // clockwise rotation about z axis by theta*i radians
+    //
+    // for odd n, the irreps are A and E1...E(nir-1)
+    // for even n, the irreps are A, B, and E1...E(nir-2)
+    //
+    gamma_[0].init(g,1,"A");
+    for (gi=0; gi < g; gi++)
+      gamma_[0].rep[gi][0][0] = 1.0;
+
+    i=1;
+
+    if (!(nt%2)) {
+      gamma_[1].init(g,1,"B");
+      for (gi=0; gi < g; gi++)
+        gamma_[1].rep[gi][0][0] = (gi%2) ? -1.0 : 1.0;
+
+      i++;
+    }
+
+    ei=1;
+    for (; i < nirrep_; i++, ei++) {
+      IrreducibleRepresentation& ir = gamma_[i];
+
+      if (nt==3 || nt==4)
+        sprintf(label,"E");
+      else
+        sprintf(label,"E%d",ei);
+
+      ir.init(g,2,label);
+      ir.complex_=1;
+
+      // identity
+      ir.rep[0].E();
+
+      // Cn
+      ir.rep[1].rotation(ei*theta);
+
+      // the other n-1 Cn's
+      for (j=2; j < g; j++)
+        ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
+    }
+
+    // identity
+    symop[0].E();
+
+    // Cn
+    symop[1].rotation(theta);
+    
+    // the other n-2 Cn's
+    for (i=2; i < nt; i++)
+      symop[i] = symop[i-1].operate(symop[1]);
+    
+    for (i=0; i < nt ; i++)
+      rot[i] = trans[i] = symop[i].trace();
+
+    break;
+
+  case CNV:
+    // clockwise rotation about z axis by theta*i radians, then
+    // reflect through the xz plane
+    //
+    // for odd n, the irreps are A1, A2, and E1...E(nir-2)
+    // for even n, the irreps are A1, A2, B1, B2, and E1...E(nir-4)
+    //
+
+    gamma_[0].init(g,1,"A1");
+    gamma_[1].init(g,1,"A2");
+
+    for (gi=0; gi < nt; gi++) {
+      // Cn's
+      gamma_[0].rep[gi][0][0] = 1.0;
+      gamma_[1].rep[gi][0][0] = 1.0;
+
+      // sigma's
+      gamma_[0].rep[gi+nt][0][0] =  1.0;
+      gamma_[1].rep[gi+nt][0][0] = -1.0;
+    }
+
+    if (!(nt%2)) {
+      gamma_[2].init(g,1,"B1");
+      gamma_[3].init(g,1,"B2");
+
+      for (gi=0; gi < nt ; gi++) {
+        double ci = (gi%2) ? -1.0 : 1.0;
+        
+        // Cn's
+        gamma_[2].rep[gi][0][0] = ci;
+        gamma_[3].rep[gi][0][0] = ci;
+
+        // sigma's
+        gamma_[2].rep[gi+nt][0][0] =  ci;
+        gamma_[3].rep[gi+nt][0][0] = -ci;
+      }
+    }
+    
+    ei=1;
+    for (i = (nt%2) ? 2 : 4; i < nirrep_; i++, ei++) {
+      IrreducibleRepresentation& ir = gamma_[i];
+
+      char lab[4];
+      if (nt==3 || nt==4)
+        sprintf(lab,"E");
+      else
+        sprintf(lab,"E%d",ei);
+
+      ir.init(g,2,lab);
+
+      // identity
+      ir.rep[0].E();
+
+      // Cn
+      ir.rep[1].rotation(ei*theta);
+
+      // the other n-2 Cn's
+      for (j=2; j < nt; j++)
+        ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
+      
+      // sigma xz
+      ir.rep[nt].sigma_xz();
+
+      SymRep sr(2);
+      sr.rotation(ei*theta/2.0);
+      
+      // the other n-1 sigma's
+      for (j=nt+1; j < g; j++)
+        ir.rep[j] = ir.rep[j-1].transform(sr);
+    }
+
+    // identity
+    symop[0].E();
+    
+    // Cn
+    symop[1].rotation(theta);
+    
+    // the other n-2 Cn's
+    for (i=2; i < nt; i++)
+      symop[i] = symop[i-1].operate(symop[1]);
+    
+    // sigma xz
+    symop[nt].sigma_xz();
+
+    so.rotation(theta/2.0);
+
+    // the other n-1 sigma's
+    for (j=nt+1; j < g; j++)
+      symop[j] = symop[j-1].transform(so);
+
+    for (i=0; i < nt ; i++) {
+      rot[i] = trans[i] = symop[i].trace();
+
+      rot[i+nt] = -symop[i+nt].trace();
+      trans[i+nt] = symop[i+nt].trace();
+    }
+    
+    break;
+
+  case CNH:
+    // lockwise rotation about z axis by theta*i radians,
+    // as well as rotation-reflection about same axis
+
+    //
+    // for odd n, the irreps are A', A", E1'...E(nir/2-1)', E1"...E(nir/2-1)''
+    // for even n, the irreps are Ag, Bg, Au, Bu,
+    //                            E1g...E(nir/2-1)g, E1u...E(nir/2-1)u
+    //
+    gamma_[0].init(g,1, (nt%2) ? "A'" : "Ag", (nt%2) ? "Ap" : 0);
+    gamma_[nirrep_/2].init(g,1, (nt%2) ? "A\"" : "Au", (nt%2) ? "Ap" : 0);
+
+    for (gi=0; gi < nt; gi++) {
+      // Cn's
+      gamma_[0].rep[gi][0][0] = 1.0;
+      gamma_[nirrep_/2].rep[gi][0][0] = 1.0;
+
+      // Sn's
+      gamma_[0].rep[gi+nt][0][0] = 1.0;
+      gamma_[nirrep_/2].rep[gi+nt][0][0] = -1.0;
+    }
+
+    if (!(nt%2)) {
+      gamma_[1].init(g,1,"Bg");
+      gamma_[1+nirrep_/2].init(g,1,"Bu");
+
+      for (gi=0; gi < nt; gi++) {
+        double ci = (gi%2) ? -1.0 : 1.0;
+
+        // Cn's
+        gamma_[1].rep[gi][0][0] = ci;
+        gamma_[1+nirrep_/2].rep[gi][0][0] = ci;
+      
+        // Sn's
+        gamma_[1].rep[gi+nt][0][0] =  ci;
+        gamma_[1+nirrep_/2].rep[gi+nt][0][0] = -ci;
+      }
+    }
+
+    ei=1;
+    for (i = (nt%2) ? 1 : 2; i < nirrep_/2 ; i++, ei++) {
+      IrreducibleRepresentation& ir1 = gamma_[i];
+      IrreducibleRepresentation& ir2 = gamma_[i+nirrep_/2];
+
+      if (nt==3 || nt==4)
+        sprintf(label,(nt%2) ? "E'" : "Eg");
+      else
+        sprintf(label,"E%d%s", ei, (nt%2) ? "'" : "g");
+
+      ir1.init(g,2,label);
+
+      if (nt==3 || nt==4)
+        sprintf(label,(nt%2) ? "E\"" : "Eu");
+      else
+        sprintf(label,"E%d%s", ei, (nt%2) ? "\"" : "u");
+
+      ir2.init(g,2,label);
+
+      ir1.complex_=1;
+      ir2.complex_=1;
+
+      // identity
+      ir1.rep[0].E();
+      ir2.rep[0].E();
+
+      // Cn
+      ir1.rep[1].rotation(ei*theta);
+      ir2.rep[1].rotation(ei*theta);
+
+      for (j=2; j < nt; j++) {
+        ir1.rep[j] = ir1.rep[j-1].operate(ir1.rep[1]);
+        ir2.rep[j] = ir2.rep[j-1].operate(ir2.rep[1]);
+      }
+
+      // Sn's
+      SymRep sr(2);
+      sr.i();
+      
+      for (j=nt; j < g; j++) {
+        ir1.rep[j] = ir1.rep[j-nt];
+        ir2.rep[j] = ir2.rep[j-nt].operate(sr);
+      }
+    }
+
+    // identity
+    symop[0].E();
+
+    // Cn
+    symop[1].rotation(theta);
+    
+    // the other n-2 Cn's
+    for (i=2; i < nt; i++)
+      symop[i] = symop[i-1].operate(symop[1]);
+
+    // Sn's, for odd nt, operate on Cn's with sigma_h, for even nt,
+    // operate Cn's with i
+    if (nt%2)
+      so.sigma_h();
+    else
+      so.i();
+    
+    for (i=0; i < nt ; i++) {
+      symop[i+nt] = symop[i].operate(so);
+      
+      rot[i] = trans[i] = symop[i].trace();
+      trans[i+nt] = symop[i+nt].trace();
+      rot[i+nt] = -trans[i+nt];
+    }
+
+    break;
+
+  case SN:
+    // clockwise rotation-reflection by theta*i radians about z axis
+    //
+    // for odd n/2, the irreps are Ag, Au, E1g...E(nir/2-1)g,E1u...E(nir/2-1)u
+    // for even n/2, the irreps are A, B, E1...E(nir-2)
+    //
+    if ((nt/2)%2) {
+      gamma_[0].init(g, 1, "Ag");
+      gamma_[nirrep_/2].init(g, 1, "Au");
+
+      for (gi=0; gi < nt/2; gi++) {
+        gamma_[0].rep[gi][0][0] = 1.0;
+        gamma_[nirrep_/2].rep[gi][0][0] = 1.0;
+
+        gamma_[0].rep[gi+nt/2][0][0] =  1.0;
+        gamma_[nirrep_/2].rep[gi+nt/2][0][0] = -1.0;
+      }
+
+      ei=1;
+      for (i=1; i < nirrep_/2 ; i++, ei++) {
+        IrreducibleRepresentation& ir1 = gamma_[i];
+        IrreducibleRepresentation& ir2 = gamma_[i+nirrep_/2];
+
+        if (nt==6)
+          sprintf(label,"Eg");
+        else
+          sprintf(label,"E%dg",ei);
+
+        ir1.init(g,2,label);
+        ir1.complex_=1;
+
+        if (nt==6)
+          sprintf(label,"Eu");
+        else
+          sprintf(label,"E%du", ei);
+
+        ir2.init(g,2,label);
+        ir2.complex_=1;
+
+        // identity
+        ir1.rep[0].E();
+        ir2.rep[0].E();
+        
+        // C(n/2)
+        ir1.rep[1].rotation(ei*theta*2.0);
+        ir2.rep[1].rotation(ei*theta*2.0);
+
+        for (j=2; j < nt/2; j++) {
+          ir1.rep[j] = ir1.rep[j-1].operate(ir1.rep[1]);
+          ir2.rep[j] = ir2.rep[j-1].operate(ir2.rep[1]);
+        }
+
+        SymRep sr(2);
+        sr.i();
+      
+        // Sn
+        for (j=nt/2; j < nt; j++) {
+          ir1.rep[j] = ir1.rep[j-nt/2];
+          ir2.rep[j] = ir2.rep[j-nt/2].operate(sr);
+        }
+      }
+
+      // identity
+      symop[0].E();
+
+      // Cn
+      symop[1].rotation(2.0*theta);
+
+      for (i=2; i < nt/2 ; i++)
+        symop[i] = symop[i-1].operate(symop[1]);
+
+      so.i();
+
+      // Sn
+      for (i=nt/2; i < nt; i++)
+        symop[i] = symop[i-nt/2].operate(so);
+      
+      for (i=0; i < nt/2 ; i++) {
+        rot[i] = trans[i] = symop[i].trace();
+
+        trans[i+nt/2] = symop[i+nt/2].trace();
+        rot[i+nt/2] = -trans[i+nt/2];
+      }
+
+    } else {
+      gamma_[0].init(g, 1, "A");
+      gamma_[1].init(g, 1, "B");
+
+      for (gi=0; gi < nt; gi++) {
+        gamma_[0].rep[gi][0][0] = 1.0;
+        gamma_[1].rep[gi][0][0] = (gi%2) ? -1.0 : 1.0;
+      }
+
+      ei=1;
+      for (i=2; i < nirrep_; i++, ei++) {
+        IrreducibleRepresentation& ir = gamma_[i];
+        
+        if (nt==4)
+          sprintf(label,"E");
+        else
+          sprintf(label,"E%d",ei);
+
+        ir.init(g,2,label);
+        ir.complex_ = 1;
+
+        // identity
+        ir.rep[0].E();
+
+        // Sn
+        ir.rep[1].rotation(ei*theta);
+        
+        for (j=2; j < nt; j++)
+          ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
+      }
+
+      // identity
+      symop[0].E();
+
+      // Sn
+      symop[1].rotation(theta);
+      symop[1][2][2] = -1.0;
+
+      for (i=2; i < nt ; i++)
+        symop[i] = symop[i-1].operate(symop[1]);
+
+      for (i=0; i < nt ; i++) {
+        trans[i] = symop[i].trace();
+        rot[i] = (i%2) ? -trans[i] : trans[i];
+      }
+    }
+
+    break;
+
+  case DN:
+    // clockwise rotation about z axis, followed by C2 about x axis
+
+    // D2 is a special case
+    if (nt==2) {
+      gamma_[0].init(g,1,"A");
+      gamma_[1].init(g,1,"B1");
+      gamma_[2].init(g,1,"B2");
+      gamma_[3].init(g,1,"B3");
+
+      for (i=0; i < g; i++) {
+        gamma_[0].rep[i][0][0] = 1.0;
+        gamma_[1].rep[i][0][0] = (i < 2) ? 1.0 : -1.0;
+        gamma_[2].rep[i][0][0] = (i % 2) ? -1.0 : 1.0;
+        gamma_[3].rep[i][0][0] = (i < 2) ?
+          ((i % 2) ? -1.0 : 1.0) : ((i%2) ? 1.0 : -1.0);
+      }
+    } else {
+      // Dn is isomorphic with Cnv
+      //
+      // for odd n, the irreps are A1, A2, and E1...E(nir-2)
+      // for even n, the irreps are A1, A2, B1, B2, and E1...E(nir-4)
+      //
+      gamma_[0].init(g,1,"A1");
+      gamma_[1].init(g,1,"A2");
+
+      for (gi=0; gi < nt; gi++) {
+        // Cn's
+        gamma_[0].rep[gi][0][0] = 1.0;
+        gamma_[1].rep[gi][0][0] = 1.0;
+
+        // C2's
+        gamma_[0].rep[gi+nt][0][0] =  1.0;
+        gamma_[1].rep[gi+nt][0][0] = -1.0;
+      }
+
+      i=2;
+
+      if (!(nt%2)) {
+        gamma_[2].init(g,1,"B1");
+        gamma_[3].init(g,1,"B2");
+
+        for (gi=0; gi < nt ; gi++) {
+          double ci = (gi%2) ? -1.0 : 1.0;
+        
+          // Cn's
+          gamma_[2].rep[gi][0][0] = ci;
+          gamma_[3].rep[gi][0][0] = ci;
+
+          // sigma's
+          gamma_[2].rep[gi+nt][0][0] =  ci;
+          gamma_[3].rep[gi+nt][0][0] = -ci;
+        }
+
+        i = 4;
+      }
+
+      ei=1;
+      for (; i < nirrep_; i++, ei++) {
+        IrreducibleRepresentation& ir = gamma_[i];
+        
+        char lab[4];
+        if (nt==3 || nt==4)
+          sprintf(lab,"E");
+        else
+          sprintf(lab,"E%d",ei);
+
+        ir.init(g,2,lab);
+
+        // identity
+        ir.rep[0].E();
+
+        // Cn
+        ir.rep[1].rotation(ei*theta);
+
+        for (j=2; j < nt; j++)
+          ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
+        
+        // C2(x)
+        ir.rep[nt].c2_y();
+        
+        SymRep sr(2);
+        sr.rotation(ei*theta/2.0);
+      
+        for (j=nt+1; j < 2*nt; j++)
+          ir.rep[j] = ir.rep[j-1].transform(sr);
+      }
+    }
+    
+    // identity
+    symop[0].E();
+
+    // Cn
+    symop[1].rotation(theta);
+
+    for (i=2; i < nt; i++)
+      symop[i] = symop[i-1].operate(symop[1]);
+
+    // C2(x)
+    symop[nt].c2_y();
+    
+    so.rotation(theta/2.0);
+
+    for (i=nt+1; i < 2*nt; i++)
+      symop[i] = symop[i-1].transform(so);
+    
+    for (i=0; i < 2*nt ; i++)
+      rot[i] = trans[i] = symop[i].trace();
+
+    break;
+
+  case DND:
+    // rotation reflection about z axis by theta/2 radians, followed
+    // by c2 about x axis, then reflection through yz plane
+    //
+    // for odd n, the irreps are A1g, A2g, A1u, A2u, E1g...E(nir/2-2)g,
+    //                                               E1u...E(nir/2-2)u
+    // for even n, the irreps are A1, A2, B1, B2, E1...E(nir-4)
+    //
+
+    if (nt%2) {
+      gamma_[0].init(g,1,"A1g");
+      gamma_[1].init(g,1,"A2g");
+
+      for (gi=0; gi < g; gi++) {
+        gamma_[0].rep[gi][0][0] = 1.0;
+        gamma_[1].rep[gi][0][0] = (gi/nt==0 || gi/nt==2) ? 1.0 : -1.0;
+      }
+      
+      i=nirrep_/2;
+      j=i+1;
+      gamma_[i].init(g,1,"A1u");
+      gamma_[j].init(g,1,"A2u");
+
+      for (gi=0; gi < g/2; gi++) {
+        gamma_[i].rep[gi][0][0] = gamma_[0].rep[gi][0][0];
+        gamma_[j].rep[gi][0][0] = gamma_[1].rep[gi][0][0];
+
+        gamma_[i].rep[gi+g/2][0][0] = -gamma_[0].rep[gi][0][0];
+        gamma_[j].rep[gi+g/2][0][0] = -gamma_[1].rep[gi][0][0];
+      }
+
+      ei=1;
+    
+      for (i=2; i < nirrep_/2 ; i++, ei++) {
+        IrreducibleRepresentation& ir1 = gamma_[i];
+        IrreducibleRepresentation& ir2 = gamma_[i+nirrep_/2];
+        
+        if (nt==3) {
+          ir1.init(g,2,"Eg");
+          ir2.init(g,2,"Eu");
+        } else {
+          sprintf(label,"E%dg",ei);
+          ir1.init(g,2,label);
+          sprintf(label,"E%du",ei);
+          ir2.init(g,2,label);
+        }
+
+        // identity
+        ir1.rep[0].E();
+        
+        // Cn
+        ir1.rep[1].rotation(ei*theta);
+        
+        for (j=2; j < nt; j++)
+          ir1.rep[j] = ir1.rep[j-1].operate(ir1.rep[1]);
+        
+        // C2(x)
+        ir1.rep[nt].c2_y();
+        
+        for (j=nt+1; j < 2*nt; j++)
+          ir1.rep[j] = ir1.rep[j-1].transform(ir1.rep[1]);
+        
+        for (j=0; j < 2*nt; j++)
+          ir2.rep[j] = ir1.rep[j];
+        
+        // Sn and sigma d
+        SymRep sr(2);
+        sr.i();
+        
+        for (j=2*nt; j < g; j++) {
+          ir1.rep[j] = ir1.rep[j-2*nt];
+          ir2.rep[j] = ir2.rep[j-2*nt].operate(sr);
+        }
+      }
+
+      // identity
+      symop[0].E();
+
+      // Cn
+      symop[1].rotation(theta);
+
+      for (i=2; i < nt; i++)
+        symop[i] = symop[i-1].operate(symop[1]);
+
+      // C2(x)
+      symop[nt].c2_y();
+
+      for (i=nt+1; i < 2*nt; i++)
+        symop[i] = symop[i-1].transform(symop[1]);
+      
+      // i + n-1 S2n + n sigma
+      so.i();
+      for (i=2*nt; i < g; i++)
+        symop[i] = symop[i-2*nt].operate(so);
+      
+      for (i=0; i < g; i++) {
+        trans[i] = symop[i].trace();
+        rot[i] = (i < g/2) ? trans[i] : -trans[i];
+      }
+
+    } else { // even nt
+
+      gamma_[0].init(g,1,"A1");
+      gamma_[1].init(g,1,"A2");
+      gamma_[2].init(g,1,"B1");
+      gamma_[3].init(g,1,"B2");
+
+      for (gi=0; gi < 2*nt; gi++) {
+        // Sn
+        gamma_[0].rep[gi][0][0] = 1.0;
+        gamma_[1].rep[gi][0][0] = 1.0;
+        gamma_[2].rep[gi][0][0] = (gi%2) ? -1.0 : 1.0;
+        gamma_[3].rep[gi][0][0] = (gi%2) ? -1.0 : 1.0;
+
+        // n C2's and n sigma's
+        gamma_[0].rep[gi+2*nt][0][0] =  1.0;
+        gamma_[1].rep[gi+2*nt][0][0] = -1.0;
+        gamma_[2].rep[gi+2*nt][0][0] = (gi%2) ? -1.0 : 1.0;
+        gamma_[3].rep[gi+2*nt][0][0] = (gi%2) ? 1.0 : -1.0;
+      }
+
+      ei=1;
+      for (i=4; i < nirrep_; i++, ei++) {
+        IrreducibleRepresentation& ir = gamma_[i];
+
+        if (nt==2)
+          sprintf(label,"E");
+        else
+          sprintf(label,"E%d",ei);
+
+        ir.init(g,2,label);
+
+        // identity
+        ir.rep[0].E();
+        
+        // S2n
+        ir.rep[1].rotation(ei*theta/2.0);
+        
+        for (j=2; j < 2*nt; j++)
+          ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
+        
+        // C2(x) + sigma_d
+        ir.rep[2*nt].c2_y();
+        
+        for (j=2*nt+1; j < g; j++)
+          ir.rep[j] = ir.rep[j-1].operate(ir.rep[1]);
+      }
+
+      // identity
+      symop[0].E();
+      
+      // Sn's
+      symop[1].rotation(theta/2.0);
+      symop[1][2][2] = -1.0;
+      
+      for (i=2; i < 2*nt; i++)
+        symop[i] = symop[i-1].operate(symop[1]);
+
+      // C2(x)
+      symop[2*nt].c2_y();
+
+      for (i=2*nt+1; i < g; i++)
+        symop[i] = symop[i-1].operate(symop[1]);
+      
+      for (i=0; i < g; i++) {
+        trans[i] = symop[i].trace();
+        rot[i] = (i%2) ? -trans[i] : trans[i];
+      }
+    }
+
+    break;
+
+  case DNH:
+    // clockwise rotation and rotation-reflection about z axis,
+    // followed by c2 about x axis and then reflection
+    // through xz 
+
+    i=nirrep_/2; j=i+1;
+
+    if (nt%2) {
+      gamma_[0].init(g,1,"A1'");
+      gamma_[1].init(g,1,"A2'");
+      gamma_[i].init(g,1,"A1\"");
+      gamma_[j].init(g,1,"A2\"");
+    } else {
+      if (nt==2) {
+        gamma_[0].init(g,1,"Ag");
+        gamma_[1].init(g,1,"B1g");
+        gamma_[4].init(g,1,"Au");
+        gamma_[5].init(g,1,"B1u");
+      } else {
+        gamma_[0].init(g,1,"A1g");
+        gamma_[1].init(g,1,"A2g");
+        gamma_[i].init(g,1,"A1u");
+        gamma_[j].init(g,1,"A2u");
+      }
+    }
+
+    for (gi=0; gi < nt; gi++) {
+      // E + n-1 Cn's
+      gamma_[0].rep[gi][0][0] = gamma_[1].rep[gi][0][0] =
+        gamma_[i].rep[gi][0][0] = gamma_[j].rep[gi][0][0] = 1.0;
+
+      // n C2's
+      gamma_[0].rep[gi+nt][0][0] = gamma_[i].rep[gi+nt][0][0] =  1.0;
+      gamma_[1].rep[gi+nt][0][0] = gamma_[j].rep[gi+nt][0][0] = -1.0;
+
+      // i + n-1 S2n's
+      gamma_[0].rep[gi+2*nt][0][0] = gamma_[1].rep[gi+2*nt][0][0] =  1.0;
+      gamma_[i].rep[gi+2*nt][0][0] = gamma_[j].rep[gi+2*nt][0][0] = -1.0;
+
+      // n sigma's
+      gamma_[0].rep[gi+3*nt][0][0] = gamma_[j].rep[gi+3*nt][0][0] =  1.0;
+      gamma_[i].rep[gi+3*nt][0][0] = gamma_[1].rep[gi+3*nt][0][0] = -1.0;
+    }
+
+    if (!(nt%2)) {
+      if (nt==2) {
+        gamma_[2].init(g,1,"B2g");
+        gamma_[3].init(g,1,"B3g");
+        gamma_[6].init(g,1,"B2u");
+        gamma_[7].init(g,1,"B3u");
+      } else {
+        gamma_[2].init(g,1,"B1g");
+        gamma_[3].init(g,1,"B2g");
+        gamma_[i+2].init(g,1,"B1u");
+        gamma_[j+2].init(g,1,"B2u");
+      }
+
+      for (gi=0; gi < nt; gi++) {
+        // E + n-1 Cn's
+        gamma_[2].rep[gi][0][0] = gamma_[3].rep[gi][0][0] =
+          gamma_[i+2].rep[gi][0][0] = gamma_[j+2].rep[gi][0][0] =
+          (gi%2) ? -1.0 : 1.0;
+
+        // n C2's
+        gamma_[2].rep[gi+nt][0][0] = gamma_[i+2].rep[gi+nt][0][0] =
+          (gi%2) ? -1.0 : 1.0;
+        gamma_[3].rep[gi+nt][0][0] = gamma_[j+2].rep[gi+nt][0][0] =
+          (gi%2) ? 1.0 : -1.0;
+
+        // i + n-1 S2n's
+        gamma_[2].rep[gi+2*nt][0][0] = gamma_[3].rep[gi+2*nt][0][0] =
+          (gi%2) ? -1.0 : 1.0;
+        gamma_[i+2].rep[gi+2*nt][0][0] = gamma_[j+2].rep[gi+2*nt][0][0] =
+          (gi%2) ? 1.0 : -1.0;
+
+        // n sigma's
+        gamma_[2].rep[gi+3*nt][0][0] = gamma_[j+2].rep[gi+3*nt][0][0] =
+          (gi%2) ? -1.0 : 1.0;
+        gamma_[i+2].rep[gi+3*nt][0][0] = gamma_[3].rep[gi+3*nt][0][0] =
+          (gi%2) ? 1.0 : -1.0;
+      }
+    }
+      
+    ei=1;
+    for (i = (nt%2) ? 2 : 4; i < nirrep_/2 ; i++, ei++) {
+      IrreducibleRepresentation& ir1 = gamma_[i];
+      IrreducibleRepresentation& ir2 = gamma_[i+nirrep_/2];
+      
+      if (nt==3) {
+        ir1.init(g,2,"E'");
+        ir2.init(g,2,"E\"");
+      } else if (nt==4) {
+        ir1.init(g,2,"Eg");
+        ir2.init(g,2,"Eu");
+      } else {
+        sprintf(label,"E%d%s", ei, (nt%2) ? "'" : "g");
+        ir1.init(g,2,label);
+
+        sprintf(label,"E%d%s", ei, (nt%2) ? "\"" : "u");
+        ir2.init(g,2,label);
+      }
+
+      // identity
+      ir1.rep[0].E();
+        
+      // n-1 Cn's
+      ir1.rep[1].rotation(ei*theta);
+
+      for (j=2; j < nt; j++)
+        ir1.rep[j] = ir1.rep[j-1].operate(ir1.rep[1]);
+        
+      // n C2's
+      ir1.rep[nt].c2_y();
+        
+      SymRep sr(2);
+      sr.rotation(ei*theta/2.0);
+      
+      for (j=nt+1; j < 2*nt; j++)
+        ir1.rep[j] = ir1.rep[j-1].transform(sr);
+        
+      sr.i();
+      for (j=0; j < 2*nt; j++) {
+        ir1.rep[j+2*nt] = ir1.rep[j];
+        ir2.rep[j] = ir1.rep[j];
+        ir2.rep[j+2*nt] = ir1.rep[j].operate(sr);
+      }
+    }
+
+    // identity
+    symop[0].E();
+
+    // n-1 Cn's
+    symop[1].rotation(theta);
+    
+    for (i=2; i < nt; i++)
+      symop[i] = symop[i-1].operate(symop[1]);
+    
+    // n C2's
+    symop[nt].c2_y();
+
+    so.rotation(theta/2.0);
+    for (i=nt+1; i < 2*nt; i++)
+      symop[i] = symop[i-1].transform(so);
+
+    if (nt%2)
+      so.sigma_h();
+    else
+      so.i();
+
+    for (i=2*nt; i < g; i++)
+      symop[i] = symop[i-2*nt].operate(so);
+    
+    for (i=0,j=2*nt; i < 2*nt ; i++,j++) {
+      rot[i] = trans[i] = symop[i].trace();
+      trans[j] = symop[j].trace();
+      rot[j] = -trans[j];
+    }
+
+    break;
+
+  case T:
+    t();
+    break;
+
+  case TH:
+    th();
+    break;
+
+  case TD:
+    td();
+    break;
+
+  case O:
+    o();
+    break;
+
+  case OH:
+    oh();
+    break;
+
+  case I:
+    this->i();
+    break;
+
+  case IH:
+    ih();
+    break;
+
+  default:
+    return -1;
+
+  }
+    
+  /* ok, we have the reducible representation of the rotations and
+   * translations, now let's use projection operators to find out how many
+   * rotations and translations there are in each irrep
+   */
+
+  if (pg != C1 && pg != CI && pg != CS && pg != T && pg != TD && pg != TH &&
+      pg != O && pg != OH && pg != I && pg != IH) {
+
+    for (i=0; i < nirrep_; i++) {
+      double nr=0; double nt=0;
+
+      for (j=0; j < gamma_[i].g; j++) {
+        nr += rot[j]*gamma_[i].character(j);
+        nt += trans[j]*gamma_[i].character(j);
+      }
+
+      gamma_[i].nrot_ = (int) ((nr+0.5)/gamma_[i].g);
+      gamma_[i].ntrans_ = (int) ((nt+0.5)/gamma_[i].g);
+    }
+  }
+
+  delete[] rot;
+  delete[] trans;
+  
+  // now find the inverse of each symop
+  for (gi=0; gi < g; gi++) {
+    int gj;
+    for (gj=0; gj < g; gj++) {
+      so = symop[gi].operate(symop[gj]);
+
+      // is so a unit matrix?
+      if (fabs(1.0-so[0][0]) < 1.0e-8 &&
+          fabs(1.0-so[1][1]) < 1.0e-8 &&
+          fabs(1.0-so[2][2]) < 1.0e-8) break;
+    }
+
+    if (gj==g) {
+      ExEnv::err0() << indent
+           << "make_table: uh oh, can't find inverse of " << gi << endl;
+      abort();
+    }
+
+    _inv[gi] = gj;
+  }
+    
+  return 0;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/pointgrp.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/pointgrp.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/pointgrp.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,273 @@
+//
+// pointgrp.cc
+//
+// Modifications are
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+/* pointgrp.cc -- implementation of the point group classes
+ *
+ *      THIS SOFTWARE FITS THE DESCRIPTION IN THE U.S. COPYRIGHT ACT OF A
+ *      "UNITED STATES GOVERNMENT WORK".  IT WAS WRITTEN AS A PART OF THE
+ *      AUTHOR'S OFFICIAL DUTIES AS A GOVERNMENT EMPLOYEE.  THIS MEANS IT
+ *      CANNOT BE COPYRIGHTED.  THIS SOFTWARE IS FREELY AVAILABLE TO THE
+ *      PUBLIC FOR USE WITHOUT A COPYRIGHT NOTICE, AND THERE ARE NO
+ *      RESTRICTIONS ON ITS USE, NOW OR SUBSEQUENTLY.
+ *
+ *  Author:
+ *      E. T. Seidl
+ *      Bldg. 12A, Rm. 2033
+ *      Computer Systems Laboratory
+ *      Division of Computer Research and Technology
+ *      National Institutes of Health
+ *      Bethesda, Maryland 20892
+ *      Internet: seidl@alw.nih.gov
+ *      June, 1993
+ */
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <math.h>
+
+#include <util/misc/formio.h>
+#include <util/state/stateio.h>
+#include <math/symmetry/pointgrp.h>
+
+using namespace std;
+using namespace sc;
+
+////////////////////////////////////////////////////////////////////////
+
+static ClassDesc PointGroup_cd(
+  typeid(PointGroup),"PointGroup",2,"public SavableState",
+  create<PointGroup>, create<PointGroup>, create<PointGroup>);
+
+PointGroup::PointGroup()
+  : symb(0)
+{
+  set_symbol("c1");
+  frame(0,0) = frame(1,1) = frame(2,2) = 1;
+  origin_[0] = origin_[1] = origin_[2] =0;
+}
+
+PointGroup::PointGroup(const char *s)
+  : symb(0)
+{
+  set_symbol(s);
+  frame(0,0) = frame(1,1) = frame(2,2) = 1;
+  origin_[0] = origin_[1] = origin_[2] =0;
+}
+
+PointGroup::PointGroup(const char *s, SymmetryOperation& so)
+  : symb(0)
+{
+  set_symbol(s);
+  frame = so;
+  origin_[0] = origin_[1] = origin_[2] =0;
+}
+
+PointGroup::PointGroup(const char *s, SymmetryOperation& so,
+                       const SCVector3& origin)
+  : symb(0)
+{
+  set_symbol(s);
+  frame = so;
+  origin_ = origin;
+}
+
+PointGroup::PointGroup(const Ref<KeyVal>& kv)
+  : symb(0)
+{
+  if (kv->exists("symmetry")) {
+    char *tmp = kv->pcharvalue("symmetry");
+    set_symbol(tmp);
+    delete[] tmp;
+  }
+  else
+    set_symbol("c1");
+
+  if (kv->exists("symmetry_frame")) {
+    for (int i=0; i < 3; i++)
+      for (int j=0; j < 3; j++) 
+        frame(i,j) = kv->doublevalue("symmetry_frame",i,j);
+  } else {
+    frame(0,0) = frame(1,1) = frame(2,2) = 1;
+  }
+
+  if (kv->exists("origin")) {
+    for (int i=0; i < 3; i++)
+      origin_[i] = kv->doublevalue("origin",i);
+  } else {
+    origin_[0] = origin_[1] = origin_[2] =0;
+  }
+}
+
+PointGroup::PointGroup(StateIn& si) :
+  SavableState(si),
+  symb(0)
+{
+  int i;
+  if (si.version(::class_desc<PointGroup>()) < 2) {
+    ExEnv::errn() << "PointGroup: checkpoint file is too old: cannot read"
+                 << endl;
+    abort();
+  }
+  else {
+    for (i=0; i<3; i++) si.get(origin_[i]);
+  }
+
+  si.getstring(symb);
+  for (i=0; i < 3; i++)
+    for (int j=0; j < 3; j++)
+      si.get(frame(i,j));
+}
+
+PointGroup::PointGroup(const PointGroup& pg)
+  : symb(0)
+{
+  *this = pg;
+}
+
+PointGroup::PointGroup(const Ref<PointGroup>& pg)
+  : symb(0)
+{
+  *this = *pg.pointer();
+}
+
+PointGroup::~PointGroup()
+{
+  if (symb) { delete[] symb; symb=0; }
+}
+
+PointGroup&
+PointGroup::operator=(const PointGroup& pg)
+{
+  set_symbol(pg.symb);
+  frame = pg.frame;
+  origin_ = pg.origin_;
+  return *this;
+}
+
+void
+PointGroup::set_symbol(const char *sym)
+{
+  if (sym) {
+    if (symb) delete[] symb;
+    int len;
+    symb = new char[(len=strlen(sym))+1];
+    for (int i=0; i<len; i++) symb[i] = (char) tolower(sym[i]);
+    symb[len] = '\0';
+  } else {
+    set_symbol("c1");
+  }
+}
+
+void
+PointGroup::save_data_state(StateOut& so)
+{
+  int i;
+  for (i=0; i<3; i++) so.put(origin_[i]);
+
+  so.putstring(symb);
+
+  for (i=0; i < 3; i++)
+    for (int j=0; j < 3; j++)
+      so.put(frame(i,j));
+}
+
+CharacterTable
+PointGroup::char_table() const
+{
+  CharacterTable ret(symb,frame);
+  return ret;
+}
+
+int
+PointGroup::equiv(const Ref<PointGroup> &grp, double tol) const
+{
+  if (strcmp(symb,grp->symb)) return 0;
+
+  for (int i=0; i < 3; i++) {
+    // origin isn't realy used, so don't check
+    //if (fabs(origin_[i] - grp->origin_[i]) > tol) return 0;
+    for (int j=0; j < 3; j++) {
+      if (fabs(frame(i,j) - grp->frame(i,j)) > tol) return 0;
+    }
+  }
+
+  return 1;
+}
+
+void
+PointGroup::print(ostream &o) const
+{
+  int i,j;
+
+  o << indent << "symmetry = " << symb << endl;
+
+  int unit_frame = 1;
+  int zero_origin = 1;
+  for (i=0; i<3; i++) {
+    for (j=0; j<3; j++) {
+      if (i==j && fabs(frame(i,j)-1.0) > 1.0e-10) unit_frame = 0;
+      else if (i != j && fabs(frame(i,j)) > 1.0e-10) unit_frame = 0;
+    }
+    if (fabs(origin_[i]) > 1.0e-10) zero_origin = 0;
+  }
+
+  if (!unit_frame) {
+    o << indent << "symmetry_frame = [";
+    o << incindent;
+    for (i=0; i<3; i++) {
+      o << endl << indent;
+      o << "[";
+      for (j=0; j<3; j++) {
+        o << scprintf(" % 18.16f", frame(i,j));
+      }
+      o << "]";
+    }
+    o << "]" << endl;
+    o << decindent;
+  }
+
+  if (!zero_origin) {
+    o << indent << "origin = [";
+    for (i=0; i<3; i++) {
+      o << scprintf(" % 18.16f", origin_[i]);
+    }
+    o << "]" << endl;
+  }
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/pointgrp.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/pointgrp.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/pointgrp.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,526 @@
+//
+// pointgrp.h
+//
+// Modifications are
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+/* pointgrp.h -- definition of the point group classes
+ *
+ *      THIS SOFTWARE FITS THE DESCRIPTION IN THE U.S. COPYRIGHT ACT OF A
+ *      "UNITED STATES GOVERNMENT WORK".  IT WAS WRITTEN AS A PART OF THE
+ *      AUTHOR'S OFFICIAL DUTIES AS A GOVERNMENT EMPLOYEE.  THIS MEANS IT
+ *      CANNOT BE COPYRIGHTED.  THIS SOFTWARE IS FREELY AVAILABLE TO THE
+ *      PUBLIC FOR USE WITHOUT A COPYRIGHT NOTICE, AND THERE ARE NO
+ *      RESTRICTIONS ON ITS USE, NOW OR SUBSEQUENTLY.
+ *
+ *  Author:
+ *      E. T. Seidl
+ *      Bldg. 12A, Rm. 2033
+ *      Computer Systems Laboratory
+ *      Division of Computer Research and Technology
+ *      National Institutes of Health
+ *      Bethesda, Maryland 20892
+ *      Internet: seidl@alw.nih.gov
+ *      June, 1993
+ */
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#ifndef _math_symmetry_pointgrp_h
+#define _math_symmetry_pointgrp_h
+
+#include <iostream>
+
+#include <util/class/class.h>
+#include <util/state/state.h>
+#include <util/keyval/keyval.h>
+#include <math/scmat/vector3.h>
+
+namespace sc {
+
+// //////////////////////////////////////////////////////////////////
+
+/** The SymmetryOperation class provides a 3 by 3 matrix
+    representation of a symmetry operation, such as a rotation or reflection.
+*/
+class SymmetryOperation {
+  private:
+    double d[3][3];
+
+  public:
+    SymmetryOperation();
+    SymmetryOperation(const SymmetryOperation &);
+    ~SymmetryOperation();
+
+    /// returns the trace of the transformation matrix
+    double trace() const { return d[0][0]+d[1][1]+d[2][2]; }
+
+    /// returns the i'th row of the transformation matrix
+    double* operator[](int i) { return d[i]; }
+
+    /// const version of the above
+    const double* operator[](int i) const { return d[i]; }
+
+    /** returns a reference to the (i,j)th element of the transformation
+        matrix */
+    double& operator()(int i, int j) { return d[i][j]; }
+
+    /// const version of the above
+    double operator()(int i, int j) const { return d[i][j]; }
+
+    /// zero out the symop
+    void zero() { memset(d,0,sizeof(double)*9); }
+
+    /// This operates on this with r (i.e. return r * this).
+    SymmetryOperation operate(const SymmetryOperation& r) const;
+
+    /// This performs the transform r * this * r~
+    SymmetryOperation transform(const SymmetryOperation& r) const;
+    
+    /// Set equal to a unit matrix
+    void unit() { zero(); d[0][0] = d[1][1] = d[2][2] = 1.0; }
+
+    /// Set equal to E
+    void E() { unit(); }
+    
+    /// Set equal to an inversion
+    void i() { zero(); d[0][0] = d[1][1] = d[2][2] = -1.0; }
+
+    /// Set equal to reflection in xy plane
+    void sigma_h() { unit(); d[2][2] = -1.0; }
+
+    /// Set equal to reflection in xz plane
+    void sigma_xz() { unit(); d[1][1] = -1.0; }
+
+    /// Set equal to reflection in yz plane
+    void sigma_yz() { unit(); d[0][0] = -1.0; }
+
+    /// Set equal to a clockwise rotation by 2pi/n
+    void rotation(int n);
+    void rotation(double theta);
+    
+    /// Set equal to C2 about the x axis
+    void c2_x() { i(); d[0][0] = 1.0; }
+
+    /// Set equal to C2 about the x axis
+    void c2_y() { i(); d[1][1] = 1.0; }
+
+    void transpose();
+
+    /// print the matrix 
+    void print(std::ostream& =ExEnv::out0()) const;
+};
+
+// //////////////////////////////////////////////////////////////////
+
+/** The SymRep class provides an n dimensional matrix representation of a
+    symmetry operation, such as a rotation or reflection.  The trace of a
+    SymRep can be used as the character for that symmetry operation.  d is
+    hardwired to 5x5 since the H irrep in Ih is 5 dimensional.
+*/
+class SymRep {
+  private:
+    int n;
+    double d[5][5];
+
+  public:
+    SymRep(int =0);
+    SymRep(const SymmetryOperation&);
+    ~SymRep();
+
+    /// Cast to a SymmetryOperation.
+    operator SymmetryOperation() const;
+    
+    /// returns the trace of the transformation matrix
+    inline double trace() const;
+
+    /// set the dimension of d
+    void set_dim(int i) { n=i; }
+    
+    /// returns the i'th row of the transformation matrix
+    double* operator[](int i) { return d[i]; }
+    /// const version of the above
+    const double* operator[](int i) const { return d[i]; }
+
+    /** returns a reference to the (i,j)th element of the transformation
+        matrix */
+    double& operator()(int i, int j) { return d[i][j]; }
+    /// const version of double& operator()(int i, int j)
+    double operator()(int i, int j) const { return d[i][j]; }
+
+    /// zero out the symop
+    void zero() { memset(d,0,sizeof(double)*25); }
+
+    /// This operates on this with r (i.e. return r * this).
+    SymRep operate(const SymRep& r) const;
+
+    /// This performs the transform r * this * r~
+    SymRep transform(const SymRep& r) const;
+    
+    /// Set equal to a unit matrix
+    void unit() {
+      zero(); d[0][0] = d[1][1] = d[2][2] = d[3][3] = d[4][4] = 1.0;
+    }
+    
+    /// Set equal to the identity
+    void E() { unit(); }
+    
+    /// Set equal to an inversion
+    void i() { zero(); d[0][0] = d[1][1] = d[2][2] = d[3][3] = d[4][4] = -1.0;}
+
+    /// Set equal to reflection in xy plane
+    void sigma_h();
+
+    /// Set equal to reflection in xz plane
+    void sigma_xz();
+
+    /// Set equal to reflection in yz plane
+    void sigma_yz();
+
+    /// Set equal to a clockwise rotation by 2pi/n
+    void rotation(int n);
+    void rotation(double theta);
+    
+    /// Set equal to C2 about the x axis
+    void c2_x();
+
+    /// Set equal to C2 about the x axis
+    void c2_y();
+
+    /// print the matrix 
+    void print(std::ostream& =ExEnv::out0()) const;
+};
+
+inline double
+SymRep::trace() const
+{
+  double r=0;
+  for (int i=0; i < n; i++)
+    r += d[i][i];
+  return r;
+}
+
+// //////////////////////////////////////////////////////////////////
+
+
+class CharacterTable;
+
+/** The IrreducibleRepresentation class provides information associated
+    with a particular irreducible representation of a point group.  This
+    includes the Mulliken symbol for the irrep, the degeneracy of the
+    irrep, the characters which represent the irrep, and the number of
+    translations and rotations in the irrep.  The order of the point group
+    is also provided (this is equal to the number of characters in an
+    irrep).  */
+class IrreducibleRepresentation {
+  friend class CharacterTable;
+
+  private:
+    int g;         // the order of the group
+    int degen;     // the degeneracy of the irrep
+    int nrot_;     // the number of rotations in this irrep
+    int ntrans_;   // the number of translations in this irrep
+    int complex_;  // true if this irrep has a complex representation
+    char *symb;    // mulliken symbol for this irrep
+    char *csymb;    // mulliken symbol for this irrep w/o special characters
+
+    SymRep *rep;   // representation matrices for the symops
+
+  public:
+    IrreducibleRepresentation();
+    IrreducibleRepresentation(const IrreducibleRepresentation&);
+    /** This constructor takes as arguments the order of the point group,
+     the degeneracy of the irrep, and the Mulliken symbol of the irrep.
+     The Mulliken symbol is copied internally. */
+    IrreducibleRepresentation(int,int,const char*,const char* =0);
+
+    ~IrreducibleRepresentation();
+
+    IrreducibleRepresentation& operator=(const IrreducibleRepresentation&);
+
+    /// Initialize the order, degeneracy, and Mulliken symbol of the irrep.
+    void init(int =0, int =0, const char* =0, const char* =0);
+    
+    /// Returns the order of the group.
+    int order() const { return g; }
+
+    /// Returns the degeneracy of the irrep.
+    int degeneracy() const { return degen; }
+
+    /// Returns the value of complex_.
+    int complex() const { return complex_; }
+
+    /// Returns the number of projection operators for the irrep.
+    int nproj() const { return degen*degen; }
+
+    /// Returns the number of rotations associated with the irrep.
+    int nrot() const { return nrot_; }
+
+    /// Returns the number of translations associated with the irrep.
+    int ntrans() const { return ntrans_; }
+
+    /// Returns the Mulliken symbol for the irrep.
+    const char * symbol() const { return symb; }
+
+    /** Returns the Mulliken symbol for the irrep without special
+        characters.
+    */
+    const char * symbol_ns() const { return (csymb?csymb:symb); }
+
+    /** Returns the character for the i'th symmetry operation of the point
+     group. */
+    double character(int i) const {
+      return complex_ ? 0.5*rep[i].trace() : rep[i].trace();
+    }
+
+    /// Returns the element (x1,x2) of the i'th representation matrix.
+    double p(int x1, int x2, int i) const { return rep[i](x1,x2); }
+    
+    /** Returns the character for the d'th contribution to the i'th
+     representation matrix. */
+    double p(int d, int i) const {
+      int dc=d/degen; int dr=d%degen;
+      return rep[i](dr,dc);
+    }
+
+    /** This prints the irrep to the given file, or stdout if none is
+     given.  The second argument is an optional string of spaces to offset
+     by. */
+    void print(std::ostream& =ExEnv::out0()) const;
+};
+
+// ///////////////////////////////////////////////////////////
+/** The CharacterTable class provides a workable character table
+ for all of the non-cubic point groups.  While I have tried to match the
+ ordering in Cotton's book, I don't guarantee that it is always followed.
+ It shouldn't matter anyway.  Also note that I don't lump symmetry
+ operations of the same class together.  For example, in C3v there are two
+ distinct C3 rotations and 3 distinct reflections, each with a separate
+ character.  Thus symop has 6 elements rather than the 3 you'll find in
+ most published character tables. */
+class CharacterTable {
+  public:
+    enum pgroups {C1, CS, CI, CN, CNV, CNH, DN, DND, DNH, SN, T, TH, TD, O,
+                  OH, I, IH};
+
+  private:
+    int g;                               // the order of the point group
+    int nt;                              // order of the princ rot axis
+    pgroups pg;                          // the class of the point group
+    int nirrep_;                         // the number of irreps in this pg
+    IrreducibleRepresentation *gamma_;   // an array of irreps
+    SymmetryOperation *symop;            // the matrices describing sym ops
+    int *_inv;                           // index of the inverse symop
+    char *symb;                          // the Schoenflies symbol for the pg
+
+    /// this determines what type of point group we're dealing with
+    int parse_symbol();
+    /// this fills in the irrep and symop arrays.
+    int make_table();
+
+    // these create the character tables for the cubic groups
+    void t();
+    void th();
+    void td();
+    void o();
+    void oh();
+    void i();
+    void ih();
+
+  public:
+    CharacterTable();
+    /** This constructor takes the Schoenflies symbol of a point group as
+        input. */
+    CharacterTable(const char*);
+    /** This is like the above, but it also takes a reference to a
+        SymmetryOperation which is the frame of reference.  All symmetry
+        operations are transformed to this frame of reference. */
+    CharacterTable(const char*,const SymmetryOperation&);
+
+    CharacterTable(const CharacterTable&);
+    ~CharacterTable();
+
+    CharacterTable& operator=(const CharacterTable&);
+
+    /// Returns the number of irreps.
+    int nirrep() const { return nirrep_; }
+    /// Returns the order of the point group
+    int order() const { return g; }
+    /// Returns the Schoenflies symbol for the point group
+    const char * symbol() const { return symb; }
+    /// Returns the i'th irrep.
+    IrreducibleRepresentation& gamma(int i) { return gamma_[i]; }
+    /// Returns the i'th symmetry operation.
+    SymmetryOperation& symm_operation(int i) { return symop[i]; }
+
+    /** Cn, Cnh, Sn, T, and Th point groups have complex representations.
+        This function returns 1 if the point group has a complex
+        representation, 0 otherwise. */
+    int complex() const {
+      if (pg==CN || pg==SN || pg==CNH || pg==T || pg==TH)
+        return 1;
+      return 0;
+    }
+
+    /// Returns the index of the symop which is the inverse of symop[i].
+    int inverse(int i) const { return _inv[i]; }
+    
+    int ncomp() const {
+      int ret=0;
+      for (int i=0; i < nirrep_; i++) {
+        int nc = (gamma_[i].complex()) ? 1 : gamma_[i].degen;
+        ret += nc;
+      }
+      return ret;
+    }
+
+    /// Returns the irrep component i belongs to.
+    int which_irrep(int i) {
+      for (int ir=0, cn=0; ir < nirrep_; ir++) {
+        int nc = (gamma_[ir].complex()) ? 1 : gamma_[ir].degen;
+        for (int c=0; c < nc; c++,cn++)
+          if (cn==i)
+            return ir;
+      }
+      return -1;
+    }
+
+    /// Returns which component i is.
+    int which_comp(int i) {
+      for (int ir=0, cn=0; ir < nirrep_; ir++) {
+        int nc = (gamma_[ir].complex()) ? 1 : gamma_[ir].degen;
+        for (int c=0; c < nc; c++,cn++)
+          if (cn==i)
+            return c;
+      }
+      return -1;
+    }
+    
+    /// This prints the irrep to the given file, or stdout if none is given.
+    void print(std::ostream& =ExEnv::out0()) const;
+};
+
+// ///////////////////////////////////////////////////////////
+
+/** The PointGroup class is really a place holder for a CharacterTable.  It
+ contains a string representation of the Schoenflies symbol of a point
+ group, a frame of reference for the symmetry operation transformation
+ matrices, and a point of origin.  The origin is not respected by the
+ symmetry operations, so if you want to use a point group with a nonzero
+ origin, first translate all your coordinates to the origin and then set
+ the origin to zero.  */
+class PointGroup: public SavableState {
+  private:
+    char *symb;
+    SymmetryOperation frame;
+    SCVector3 origin_;
+
+  public:
+    PointGroup();
+    /** This constructor takes a string containing the Schoenflies symbol
+        of the point group as its only argument. */
+    PointGroup(const char*);
+    /** Like the above, but this constructor also takes a frame of reference
+        as an argument. */
+    PointGroup(const char*,SymmetryOperation&);
+    /** Like the above, but this constructor also takes a point of origin
+        as an argument. */
+    PointGroup(const char*,SymmetryOperation&,const SCVector3&);
+    /** The PointGroup KeyVal constructor looks for three keywords:
+       symmetry, symmetry_frame, and origin. symmetry is a string
+       containing the Schoenflies symbol of the point group.  origin is an
+       array of doubles which gives the x, y, and z coordinates of the
+       origin of the symmetry frame.  symmetry_frame is a 3 by 3 array of
+       arrays of doubles which specify the principal axes for the
+       transformation matrices as a unitary rotation.
+
+       For example, a simple input which will use the default origin and
+       symmetry_frame ((0,0,0) and the unit matrix, respectively), might
+       look like this:
+
+       <pre>
+       pointgrp<PointGroup>: (
+         symmetry = "c2v"
+       )
+       </pre>
+
+       By default, the principal rotation axis is taken to be the z axis.
+       If you already have a set of coordinates which assume that the
+       rotation axis is the x axis, then you'll have to rotate your frame
+       of reference with symmetry_frame:
+
+       <pre>
+       pointgrp<PointGroup>: (
+         symmetry = "c2v"
+         symmetry_frame = [
+           [ 0 0 1 ]
+           [ 0 1 0 ]
+           [ 1 0 0 ]
+         ]
+       )
+       </pre>
+     */
+    PointGroup(const Ref<KeyVal>&);
+
+    PointGroup(StateIn&);
+    PointGroup(const PointGroup&);
+    PointGroup(const Ref<PointGroup>&);
+    ~PointGroup();
+
+    PointGroup& operator=(const PointGroup&);
+
+    /// Returns 1 if the point groups are equivalent, 0 otherwise.
+    int equiv(const Ref<PointGroup> &, double tol = 1.0e-6) const;
+
+    /// Returns the CharacterTable for this point group.
+    CharacterTable char_table() const;
+    /// Returns the Schoenflies symbol for this point group.
+    const char * symbol() const { return symb; }
+    /// Returns the frame of reference for this point group.
+    SymmetryOperation& symm_frame() { return frame; }
+    /// A const version of the above
+    const SymmetryOperation& symm_frame() const { return frame; }
+    /// Returns the origin of the symmetry frame.
+    SCVector3& origin() { return origin_; }
+    const SCVector3& origin() const { return origin_; }
+
+    /// Sets (or resets) the Schoenflies symbol.
+    void set_symbol(const char*);
+
+    void save_data_state(StateOut& so);
+
+    void print(std::ostream&o=ExEnv::out0()) const;
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/rep.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/rep.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/rep.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,292 @@
+//
+// rep.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#include <util/misc/math.h>
+
+#include <math/symmetry/pointgrp.h>
+#include <util/misc/formio.h>
+
+using namespace std;
+using namespace sc;
+
+/////////////////////////////////////////////////////////////////////////
+
+SymRep::SymRep(int i) :
+  n(i)
+{
+  zero();
+}
+
+SymRep::SymRep(const SymmetryOperation& so) :
+  n(3)
+{
+  memset(d,0,sizeof(double)*25);
+  for (int i=0; i < 3; i++)
+    for (int j=0; j < 3; j++)
+      d[i][j] = so[i][j];
+}
+
+SymRep::~SymRep()
+{
+  n=0;
+}
+
+SymRep::operator SymmetryOperation() const
+{
+  if (n != 3) {
+    ExEnv::err0() << indent << "SymRep::operator SymmetryOperation(): "
+         << "trying to cast to symop when n == " << n << endl;
+    abort();
+  }
+    
+  SymmetryOperation so;
+
+  for (int i=0; i < 3; i++)
+    for (int j=0; j < 3; j++)
+      so[i][j] = d[i][j];
+
+  return so;
+}
+
+SymRep
+SymRep::operate(const SymRep& r) const
+{
+  if (r.n != n) {
+    ExEnv::err0() << indent << "SymRep::operate(): dimensions don't match: "
+         << r.n << " != " << n << endl;
+    abort();
+  }
+  
+  SymRep ret(n);
+
+  for (int i=0; i < n; i++) {
+    for (int j=0; j < n; j++) {
+      double t=0;
+      for (int k=0; k < n; k++)
+        t += r[i][k] * d[k][j];
+      ret[i][j] = t;
+    }
+  }
+
+  return ret;
+}
+
+SymRep
+SymRep::transform(const SymRep& r) const
+{
+  int i,j,k;
+
+  if (r.n != n) {
+    ExEnv::err0() << indent
+         << "SymRep::symm_transform(): dimensions don't match: "
+         << r.n << " != " << n << endl;
+    abort();
+  }
+  
+  SymRep ret(n), foo(n);
+
+  // foo = r * d
+  for (i=0; i < n; i++) {
+    for (j=0; j < n; j++) {
+      double t=0;
+      for (k=0; k < n; k++)
+        t += r[i][k] * d[k][j];
+      foo[i][j] = t;
+    }
+  }
+
+  // ret = (r*d)*r~ = foo*r~
+  for (i=0; i < n; i++) {
+    for (j=0; j < n; j++) {
+      double t=0;
+      for (k=0; k < n; k++)
+        t += foo[i][k]*r[j][k];
+      ret[i][j]=t;
+    }
+  }
+
+  return ret;
+}
+
+void
+SymRep::sigma_h()
+{
+  unit();
+
+  if (n==3) {
+    d[2][2] = -1.0;
+  } else if (n==5) {
+    d[3][3] = d[4][4] = -1.0;
+  }
+}
+  
+void
+SymRep::sigma_xz()
+{
+  unit();
+
+  if (n==2 || n==3 || n==4) {
+    d[1][1] = -1.0;
+    if (n==4)
+      d[2][2] = -1.0;
+  } else if (n==5) {
+    d[2][2] = d[4][4] = -1.0;
+  }
+}
+
+void
+SymRep::sigma_yz()
+{
+  unit();
+
+  if (n==2 || n==3 || n==4) {
+    d[0][0] = -1.0;
+    if (n==4)
+      d[3][3] = -1.0;
+  } else if (n==5) {
+    d[2][2] = d[3][3] = -1.0;
+  }
+}
+  
+void
+SymRep::rotation(int nt)
+{
+  double theta = (nt) ? 2.0*M_PI/nt : 2.0*M_PI;
+  rotation(theta);
+}
+
+void
+SymRep::rotation(double theta)
+{
+  zero();
+
+  double ctheta = cos(theta);
+  double stheta = sin(theta);
+  double c2theta = cos(2*theta);
+  double s2theta = sin(2*theta);
+
+  switch (n) {
+  case 1:
+    d[0][0] = 1.0;
+    break;
+    
+  case 3:
+    d[0][0] = ctheta;
+    d[0][1] = stheta;
+    d[1][0] = -stheta;
+    d[1][1] = ctheta;
+    d[2][2] = 1.0;
+    break;
+
+  case 4:
+  case 2:
+    d[0][0] = ctheta;
+    d[0][1] = stheta;
+    d[1][0] = -stheta;
+    d[1][1] = ctheta;
+    
+    // this is ok since d is hardwired
+    d[2][2] = c2theta;
+    d[2][3] = -s2theta;
+    d[3][2] = s2theta;
+    d[3][3] = c2theta;
+    break;
+    
+  case 5:
+    d[0][0] = 1.0;
+    
+    d[1][1] = c2theta;
+    d[1][2] = s2theta;
+    d[2][1] = -s2theta;
+    d[2][2] = c2theta;
+
+    d[3][3] = ctheta;
+    d[3][4] = -stheta;
+    d[4][3] = stheta;
+    d[4][4] = ctheta;
+    break;
+
+  default:
+    ExEnv::err0() << indent << "SymRep::rotation(): n > 5 (" << n << ")\n";
+    abort();
+  }
+  
+}
+
+void
+SymRep::c2_x()
+{
+  i();
+
+  if (n==2 || n==3 || n==4) {
+    d[0][0] = 1.0;
+    if (n==4)
+      d[3][3] = 1.0;
+  } else if (n==5) {
+    d[0][0] = d[1][1] = d[4][4] = 1.0;
+  }
+}
+  
+void
+SymRep::c2_y()
+{
+  i();
+
+  if (n==2 || n==3 || n==4) {
+    d[1][1] = 1.0;
+    if (n==4)
+      d[2][2] = 1.0;
+  } else if (n==5) {
+    d[0][0] = d[1][1] = d[3][3] = 1.0;
+  }
+}
+  
+
+void
+SymRep::print(ostream& os) const
+{
+  int i;
+
+  os << indent;
+  for (i=0; i < n; i++) os << scprintf("%11d",i+1);
+  os << endl;
+  
+  for (i=0; i < n; i++) {
+    os << indent << scprintf("%3d ",i+1);
+    for (int j=0; j < n; j++)
+      os << scprintf(" %10.7f",d[i][j]);
+    os << endl;
+  }
+  os << endl;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/symmetry.dox
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/symmetry.dox	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/symmetry.dox	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,12 @@
+
+/** \page symmetry The Symmetry Library
+
+The symmetry library implements a few basic ideas in group theory.  Namely,
+it provides a character table for a given point group.  Currently it does
+not support the cubic point groups (e.g. \f$T_d\f$ or \f$I_h\f$).  The only
+classes intended for public consumption are the CharacterTable and
+PointGroup classes.  PointGroup stores the Schoenflies symbol for a point
+group, and can generate a CharacterTable on the fly.  A CharacterTable
+contains IrreducibleRepresentation's and SymmetryOperation's.
+
+*/
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/symop.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/symop.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/symop.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,158 @@
+//
+// symop.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#include <util/misc/math.h>
+
+#include <math/symmetry/pointgrp.h>
+#include <util/misc/formio.h>
+
+using namespace std;
+using namespace sc;
+
+////////////////////////////////////////////////////////////////////////
+
+SymmetryOperation::SymmetryOperation()
+{
+  zero();
+}
+
+SymmetryOperation::SymmetryOperation(const SymmetryOperation &so)
+{
+  for (int i=0; i<3; i++) {
+    for (int j=0; j<3; j++) {
+      d[i][j] = so.d[i][j];
+    }
+  }
+}
+
+SymmetryOperation::~SymmetryOperation()
+{
+}
+
+SymmetryOperation
+SymmetryOperation::operate(const SymmetryOperation& r) const
+{
+  SymmetryOperation ret;
+  for (int i=0; i < 3; i++)
+    for (int j=0; j < 3; j++) {
+      double t=0;
+      for (int k=0; k < 3; k++)
+        t += r.d[i][k]*d[k][j];
+      ret.d[i][j] = t;
+    }
+  return ret;
+}
+
+SymmetryOperation
+SymmetryOperation::transform(const SymmetryOperation& r) const
+{
+  int i,j,k;
+  SymmetryOperation ret,foo;
+  
+  // foo = r * d
+  for (i=0; i < 3; i++) {
+    for (j=0; j < 3; j++) {
+      double t=0;
+      for (k=0; k < 3; k++)
+        t += r.d[i][k] * d[k][j];
+      foo.d[i][j] = t;
+    }
+  }
+
+  // ret = (r*d)*r~ = foo*r~
+  for (i=0; i < 3; i++) {
+    for (j=0; j < 3; j++) {
+      double t=0;
+      for (k=0; k < 3; k++)
+        t += foo.d[i][k]*r.d[j][k];
+      ret.d[i][j]=t;
+    }
+  }
+
+  return ret;
+}
+
+// Clockwise rotation by 2pi/n degrees
+void
+SymmetryOperation::rotation(int n)
+{
+  double theta = (n) ? 2.0*M_PI/n : 2.0*M_PI;
+  rotation(theta);
+}
+
+// Clockwise rotation by theta degrees
+void
+SymmetryOperation::rotation(double theta)
+{
+  zero();
+
+  double ctheta = cos(theta);
+  double stheta = sin(theta);
+
+  d[0][0] = ctheta;
+  d[0][1] = stheta;
+  d[1][0] = -stheta;
+  d[1][1] = ctheta;
+  d[2][2] = 1.0;
+}
+
+void
+SymmetryOperation::transpose()
+{
+  for (int i=1; i<3; i++) {
+    for (int j=0; j<i; j++) {
+      double tmp = d[i][j];
+      d[i][j] = d[j][i];
+      d[j][i] = tmp;
+    }
+  }
+}
+
+void
+SymmetryOperation::print(ostream& os) const
+{
+  os << indent << "        1          2          3\n"
+     << indent << "  1  "
+     << scprintf("%10.7f ", d[0][0])
+     << scprintf("%10.7f ", d[0][1])
+     << scprintf("%10.7f ", d[0][2]) << endl
+     << indent << "  2  "
+     << scprintf("%10.7f ", d[1][0])
+     << scprintf("%10.7f ", d[1][1])
+     << scprintf("%10.7f ", d[1][2]) << endl
+     << indent << "  3  "
+     << scprintf("%10.7f ", d[2][0])
+     << scprintf("%10.7f ", d[2][1])
+     << scprintf("%10.7f ", d[2][2]) << endl << endl;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/testpg.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/testpg.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/testpg.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,40 @@
+
+#include <math/symmetry/pointgrp.h>
+#include <math/symmetry/corrtab.h>
+
+using namespace sc;
+
+main(int argc, char *argv[])
+{
+  Ref<PointGroup> pg = new PointGroup(argv[1]);
+  Ref<PointGroup> pg2;
+  if (argc > 2) pg2 = new PointGroup(argv[2]);
+
+  //pg.char_table().print();
+
+  CharacterTable ct = pg->char_table();
+  CharacterTable ct2;
+  ct2=ct;
+
+  ct.print();
+
+  if (pg2.nonnull()) {
+      pg2->char_table().print();
+      if (argc <= 3) {
+          CorrelationTable corrtab(pg,pg2);
+          corrtab.print();
+        }
+    }
+
+  // test given axis rearrangements
+  for (int i=3; i<argc; i++) {
+      pg2->symm_frame().zero();
+      for (int j=0; j<3; j++) {
+          if (argv[i][j] == 'x') pg2->symm_frame()(j,0) = 1.0;
+          else if (argv[i][j] == 'y') pg2->symm_frame()(j,1) = 1.0;
+          else if (argv[i][j] == 'z') pg2->symm_frame()(j,2) = 1.0;
+        }
+      CorrelationTable corrtab(pg,pg2);
+      corrtab.print();
+    }
+}
Index: ThirdParty/mpqc_open/src/lib/math/symmetry/tetra.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/symmetry/tetra.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/symmetry/tetra.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,573 @@
+//
+// tetra.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Edward Seidl <seidl@janed.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#include <util/misc/math.h>
+#include <string.h>
+
+#include <math/symmetry/pointgrp.h>
+
+using namespace sc;
+
+// these are the operations which make up T
+static void
+t_ops(SymmetryOperation *symop)
+{
+  // identity
+  symop[0].E();
+
+  // C2(x)
+  symop[9].c2_x();
+
+  // C2(y)
+  symop[10].c2_y();
+
+  // C2(z)
+  symop[11].rotation((double)M_PI);
+
+  // a = ( 1, 1, 1)
+  // b = (-1,-1, 1)
+  // c = ( 1,-1,-1)
+  // d = (-1, 1,-1)
+  // C3 (a)
+  symop[1][0][2] =  1.0;
+  symop[1][1][0] =  1.0;
+  symop[1][2][1] =  1.0;
+
+  // C3 (b)
+  symop[2] = symop[1].transform(symop[11]);
+
+  // C3 (c)
+  symop[3] = symop[1].transform(symop[9]);
+
+  // C3 (d)
+  symop[4] = symop[1].transform(symop[10]);
+
+  // C3^2 (a)
+  symop[5][0][1] =  1.0;
+  symop[5][1][2] =  1.0;
+  symop[5][2][0] =  1.0;
+
+  // C3^2 (b)
+  symop[6] = symop[5].transform(symop[11]);
+
+  // C3^2 (c)
+  symop[7] = symop[5].transform(symop[9]);
+
+  // C3^2 (d)
+  symop[8] = symop[5].transform(symop[10]);
+}
+
+// this gives us the operations in Td which aren't in T.
+static void
+td_ops(SymmetryOperation *symop)
+{
+  // S4 (x)
+  symop[0][0][0] = -1.0;
+  symop[0][1][2] = -1.0;
+  symop[0][2][1] =  1.0;
+
+  // S4^3 (x)
+  symop[1][0][0] = -1.0;
+  symop[1][1][2] =  1.0;
+  symop[1][2][1] = -1.0;
+
+  // S4 (y)
+  symop[2][0][2] =  1.0;
+  symop[2][1][1] = -1.0;
+  symop[2][2][0] = -1.0;
+
+  // S4^3 (y)
+  symop[3][0][2] = -1.0;
+  symop[3][1][1] = -1.0;
+  symop[3][2][0] =  1.0;
+
+  // S4 (z)
+  symop[4][0][1] = -1.0;
+  symop[4][1][0] =  1.0;
+  symop[4][2][2] = -1.0;
+
+  // S4^3 (z)
+  symop[5][0][1] =  1.0;
+  symop[5][1][0] = -1.0;
+  symop[5][2][2] = -1.0;
+
+  // a = ( 1, 1, 1)
+  // b = (-1,-1, 1)
+  // c = ( 1,-1,-1)
+  // d = (-1, 1,-1)
+  // sigma (ac)
+  symop[6][0][0] =  1.0;
+  symop[6][1][2] =  1.0;
+  symop[6][2][1] =  1.0;
+
+  // sigma (bd)
+  symop[7][0][0] =  1.0;
+  symop[7][1][2] = -1.0;
+  symop[7][2][1] = -1.0;
+
+  // sigma (ad)
+  symop[8][0][2] =  1.0;
+  symop[8][1][1] =  1.0;
+  symop[8][2][0] =  1.0;
+
+  // sigma (bc)
+  symop[9][0][2] = -1.0;
+  symop[9][1][1] =  1.0;
+  symop[9][2][0] = -1.0;
+
+  // sigma (ab)
+  symop[10][0][1] =  1.0;
+  symop[10][1][0] =  1.0;
+  symop[10][2][2] =  1.0;
+
+  // sigma (dc)
+  symop[11][0][1] = -1.0;
+  symop[11][1][0] = -1.0;
+  symop[11][2][2] =  1.0;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+void
+CharacterTable::t()
+{
+  // t_ops gives us all the symmetry operations we need
+  t_ops(symop);
+
+  int i;
+
+  gamma_[0].init(g,1,"A");
+  for (i=0; i < g; i++)
+    gamma_[0].rep[i][0][0] = 1.0;
+
+  IrreducibleRepresentation& ire = gamma_[1];
+  ire.init(g,2,"E");
+  ire.complex_=1;
+
+  IrreducibleRepresentation& irt = gamma_[2];
+  irt.init(g,3,"T");
+  irt.nrot_ = 1;
+  irt.ntrans_ = 1;
+
+  // the symmetry operation matrices give us a basis for irrep T
+  for (i=0; i < g; i++)
+    irt.rep[i] = symop[i];
+
+  // identity
+  ire.rep[0].E();
+
+  // 4 C3's
+  ire.rep[1].rotation(2.0*(double)M_PI/3.0);
+  ire.rep[2] = ire.rep[1];
+  ire.rep[3] = ire.rep[1];
+  ire.rep[4] = ire.rep[1];
+
+  ire.rep[5] = ire.rep[1].operate(ire.rep[1]);
+  ire.rep[6] = ire.rep[5];
+  ire.rep[7] = ire.rep[5];
+  ire.rep[8] = ire.rep[5];
+  
+  // 3 C2's
+  ire.rep[9].unit();
+  ire.rep[10].unit();
+  ire.rep[11].unit();
+
+}
+
+void
+CharacterTable::th()
+{
+  int i,j;
+
+  SymmetryOperation so;
+  so.i();
+  
+  t_ops(symop);
+  for (i=0; i < 12; i++)
+    symop[i+12] = symop[i].operate(so);
+  
+  gamma_[0].init(g,1,"Ag");
+  gamma_[1].init(g,1,"Au");
+
+  for (i=0; i < 12; i++) {
+    gamma_[0].rep[i][0][0] = 1.0;
+    gamma_[1].rep[i][0][0] = 1.0;
+
+    gamma_[0].rep[i+12][0][0] =  1.0;
+    gamma_[1].rep[i+12][0][0] = -1.0;
+  }
+
+  IrreducibleRepresentation& ireg = gamma_[2];
+  IrreducibleRepresentation& ireu = gamma_[3];
+
+  IrreducibleRepresentation& irtg = gamma_[4];
+  IrreducibleRepresentation& irtu = gamma_[5];
+
+  ireg.init(g,2,"Eg");
+  ireu.init(g,2,"Eu");
+  ireg.complex_=1;
+  ireu.complex_=1;
+
+  irtg.init(g,3,"Tg");
+  irtu.init(g,3,"Tu");
+  irtg.nrot_=1;
+  irtu.ntrans_=1;
+
+  // the symmetry operation matrices form a basis for Tu.  Tg(g)=Tu(g) for
+  // the proper rotations, and = -Tu(g) for the improper ones
+  for (i=0; i < 12; i++) {
+    irtg.rep[i] = symop[i];
+    irtu.rep[i] = symop[i];
+
+    irtg.rep[i+12] = symop[i];
+    irtu.rep[i+12] = symop[i+12];
+  }
+    
+  // identity
+  ireg.rep[0].E();
+  
+  // 4 C3's
+  ireg.rep[1].rotation(2.0*(double)M_PI/3.0);
+  ireg.rep[2] = ireg.rep[1];
+  ireg.rep[3] = ireg.rep[1];
+  ireg.rep[4] = ireg.rep[1];
+
+  // 4 C3^2's
+  ireg.rep[5] = ireg.rep[1].operate(ireg.rep[1]);
+  ireg.rep[6] = ireg.rep[5];
+  ireg.rep[7] = ireg.rep[5];
+  ireg.rep[8] = ireg.rep[5];
+
+  // 3 C2's
+  ireg.rep[9].unit();
+  ireg.rep[10].unit();
+  ireg.rep[11].unit();
+
+  SymRep sr(2);
+  sr.i();
+  
+  for (j=0; j < 12; j++) {
+    ireu.rep[j] = ireg.rep[j];
+    ireg.rep[j+12] = ireg.rep[j];
+    ireu.rep[j+12] = ireg.rep[j].operate(sr);
+  }
+}
+
+void
+CharacterTable::td()
+{
+  // first get the T operations, then the Td operations
+  t_ops(symop);
+  td_ops(&symop[12]);
+  
+  int i;
+  
+  gamma_[0].init(g,1,"A1");
+  gamma_[1].init(g,1,"A2");
+
+  for (i=0; i < 12; i++) {
+    gamma_[0].rep[i][0][0] = 1.0;
+    gamma_[1].rep[i][0][0] = 1.0;
+
+    gamma_[0].rep[i+12][0][0] =  1.0;
+    gamma_[1].rep[i+12][0][0] = -1.0;
+  }
+
+  IrreducibleRepresentation& ire = gamma_[2];
+  ire.init(g,2,"E");
+
+  IrreducibleRepresentation& irt1 = gamma_[3];
+  IrreducibleRepresentation& irt2 = gamma_[4];
+
+  irt1.init(g,3,"T1");
+  irt2.init(g,3,"T2");
+  irt1.nrot_ = 1;
+  irt2.ntrans_ = 1;
+
+  // the symmetry operation matrices form a basis for T2.  T1(g)=T2(g) for
+  // the proper rotations, and = -T2(g) for the improper ones
+  SymmetryOperation so;
+  so.i();
+  
+  for (i=0; i < 12; i++) {
+    irt1.rep[i] = symop[i];
+    irt2.rep[i] = symop[i];
+    irt1.rep[i+12] = symop[i+12].operate(so);
+    irt2.rep[i+12] = symop[i+12];
+  }
+  
+  // identity
+  ire.rep[0].E();
+
+  // 4 C3's
+  ire.rep[1].rotation(2.0*(double)M_PI/3.0);
+  ire.rep[2] = ire.rep[1];
+  ire.rep[3] = ire.rep[1];
+  ire.rep[4] = ire.rep[1];
+
+  // 4 C3^2's
+  ire.rep[5] = ire.rep[1].operate(ire.rep[1]);
+  ire.rep[6] = ire.rep[5];
+  ire.rep[7] = ire.rep[5];
+  ire.rep[8] = ire.rep[5];
+
+  // 3 C2's
+  ire.rep[9].unit();
+  ire.rep[10].unit();
+  ire.rep[11].unit();
+
+  // 6 S4's
+  ire.rep[12].c2_x();
+  ire.rep[13].c2_x();
+
+  ire.rep[14] = ire.rep[12].operate(ire.rep[1]);
+  ire.rep[15] = ire.rep[14];
+  
+  ire.rep[16] = ire.rep[14].operate(ire.rep[1]);
+  ire.rep[17] = ire.rep[16];
+
+  for (i=18; i < 24; i++)
+    ire.rep[i] = ire.rep[i-6];
+}
+
+void
+CharacterTable::o()
+{
+  int i;
+  
+  // first get the T operations, then the O operations
+  t_ops(symop);
+  td_ops(&symop[12]);
+  
+  SymmetryOperation so;
+  so.i();
+
+  for (i=12; i < 24; i++)
+    symop[i] = symop[i].operate(so);
+  
+  gamma_[0].init(g,1,"A1");
+  gamma_[1].init(g,1,"A2");
+
+  for (i=0; i < 12; i++) {
+    gamma_[0].rep[i][0][0] = 1.0;
+    gamma_[1].rep[i][0][0] = 1.0;
+
+    gamma_[0].rep[i+12][0][0] =  1.0;
+    gamma_[1].rep[i+12][0][0] = -1.0;
+  }
+
+  IrreducibleRepresentation& ire = gamma_[2];
+  ire.init(g,2,"E");
+
+  IrreducibleRepresentation& irt1 = gamma_[3];
+  IrreducibleRepresentation& irt2 = gamma_[4];
+
+  irt1.init(g,3,"T1");
+  irt2.init(g,3,"T2");
+  irt1.nrot_ = 1;
+  irt1.ntrans_ = 1;
+
+  // the symmetry operation matrices form a basis for T1.  T2(g)=T1(g) for
+  // the proper rotations, and = -T1(g) for the improper ones
+  
+  for (i=0; i < 12; i++) {
+    irt1.rep[i] = symop[i];
+    irt2.rep[i] = symop[i];
+    irt1.rep[i+12] = symop[i+12];
+    irt2.rep[i+12] = symop[i+12].operate(so);
+  }
+  
+  // identity
+  ire.rep[0].E();
+
+  // 4 C3's
+  ire.rep[1].rotation(2.0*(double)M_PI/3.0);
+  ire.rep[2] = ire.rep[1];
+  ire.rep[3] = ire.rep[1];
+  ire.rep[4] = ire.rep[1];
+
+  // 4 C3^2's
+  ire.rep[5] = ire.rep[1].operate(ire.rep[1]);
+  ire.rep[6] = ire.rep[5];
+  ire.rep[7] = ire.rep[5];
+  ire.rep[8] = ire.rep[5];
+
+  // 3 C2's
+  ire.rep[9].unit();
+  ire.rep[10].unit();
+  ire.rep[11].unit();
+
+  // 6 C4's
+  ire.rep[12].c2_x();
+  ire.rep[13].c2_x();
+
+  ire.rep[14] = ire.rep[12].operate(ire.rep[1]);
+  ire.rep[15] = ire.rep[14];
+  
+  ire.rep[16] = ire.rep[14].operate(ire.rep[1]);
+  ire.rep[17] = ire.rep[16];
+
+  // 6 C2's
+  for (i=18; i < 24; i++)
+    ire.rep[i] = ire.rep[i-6];
+}
+
+void CharacterTable::oh()
+{
+  int i,j;
+  
+  SymmetryOperation so;
+  so.i();
+  
+  // first get the T operations, then the O operations, then the Th
+  // operations, then the Td operations
+  t_ops(symop);
+  td_ops(&symop[36]);
+
+  for (i=0; i < 12; i++) {
+    symop[i+24] = symop[i].operate(so);
+    symop[i+12] = symop[i+36].operate(so);
+  }
+  
+  gamma_[0].init(g,1,"A1g");
+  gamma_[1].init(g,1,"A2g");
+  gamma_[5].init(g,1,"A1u");
+  gamma_[6].init(g,1,"A2u");
+
+  for (i=0; i < 12; i++) {
+    gamma_[0].rep[i][0][0] = 1.0;
+    gamma_[1].rep[i][0][0] = 1.0;
+    gamma_[5].rep[i][0][0] = 1.0;
+    gamma_[6].rep[i][0][0] = 1.0;
+
+    gamma_[0].rep[i+12][0][0] =  1.0;
+    gamma_[1].rep[i+12][0][0] = -1.0;
+    gamma_[5].rep[i+12][0][0] =  1.0;
+    gamma_[6].rep[i+12][0][0] = -1.0;
+
+    gamma_[0].rep[i+24][0][0] =  1.0;
+    gamma_[1].rep[i+24][0][0] =  1.0;
+    gamma_[5].rep[i+24][0][0] = -1.0;
+    gamma_[6].rep[i+24][0][0] = -1.0;
+
+    gamma_[0].rep[i+36][0][0] =  1.0;
+    gamma_[1].rep[i+36][0][0] = -1.0;
+    gamma_[5].rep[i+36][0][0] = -1.0;
+    gamma_[6].rep[i+36][0][0] =  1.0;
+  }
+
+  // the symmetry operation matrices form a basis for T1u.  T2u(g)=T1u(g) for
+  // the proper rotations, and = -T1(g) for the improper ones.
+  // T1g(g)=T1u(g) for the O part, and = -T1u(g) for the ixO part.
+  // T2g(g)=T1g(g) for proper rotations and =-T1g(g) for improper
+
+  gamma_[3].init(g,3,"T1g");
+  gamma_[4].init(g,3,"T2g");
+  gamma_[8].init(g,3,"T1u");
+  gamma_[9].init(g,3,"T2u");
+
+  gamma_[3].nrot_=1;
+  gamma_[8].ntrans_=1;
+  
+  for (i=0; i < 12; i++) {
+    gamma_[3].rep[i] = symop[i];
+    gamma_[4].rep[i] = symop[i];
+    gamma_[8].rep[i] = symop[i];
+    gamma_[9].rep[i] = symop[i];
+    
+    gamma_[3].rep[i+12] = symop[i+12];
+    gamma_[4].rep[i+12] = symop[i+12].operate(so);
+    gamma_[8].rep[i+12] = symop[i+12];
+    gamma_[9].rep[i+12] = symop[i+12].operate(so);
+    
+    gamma_[3].rep[i+24] = symop[i+24].operate(so);
+    gamma_[4].rep[i+24] = symop[i+24].operate(so);
+    gamma_[8].rep[i+24] = symop[i+24];
+    gamma_[9].rep[i+24] = symop[i+24];
+    
+    gamma_[3].rep[i+36] = symop[i+36].operate(so);
+    gamma_[4].rep[i+36] = symop[i+36];
+    gamma_[8].rep[i+36] = symop[i+36];
+    gamma_[9].rep[i+36] = symop[i+36].operate(so);
+  }
+
+  IrreducibleRepresentation& ireg = gamma_[2];
+  IrreducibleRepresentation& ireu = gamma_[7];
+
+  ireg.init(g,2,"Eg");
+  ireu.init(g,2,"Eu");
+    
+  // identity
+  ireg.rep[0].E();
+  
+  // 4 C3's
+  ireg.rep[1].rotation(2.0*(double)M_PI/3.0);
+  ireg.rep[2] = ireg.rep[1];
+  ireg.rep[3] = ireg.rep[1];
+  ireg.rep[4] = ireg.rep[1];
+
+  // 4 C3^2's
+  ireg.rep[5] = ireg.rep[1].operate(ireg.rep[1]);
+  ireg.rep[6] = ireg.rep[5];
+  ireg.rep[7] = ireg.rep[5];
+  ireg.rep[8] = ireg.rep[5];
+
+  // 3 C2's
+  ireg.rep[9].unit();
+  ireg.rep[10].unit();
+  ireg.rep[11].unit();
+
+  // 6 C4's
+  ireg.rep[12].c2_x();
+  ireg.rep[13].c2_x();
+
+  ireg.rep[14] = ireg.rep[12].operate(ireg.rep[1]);
+  ireg.rep[15] = ireg.rep[14];
+  
+  ireg.rep[16] = ireg.rep[14].operate(ireg.rep[1]);
+  ireg.rep[17] = ireg.rep[16];
+
+  // 6 C2's
+  for (i=18; i < 24; i++)
+    ireg.rep[i] = ireg.rep[i-6];
+
+  SymRep sr(2);
+  sr.i();
+  
+  for (j=0; j < 24; j++) {
+    ireu.rep[j] = ireg.rep[j];
+    ireg.rep[j+24] = ireg.rep[j];
+    ireu.rep[j+24] = ireg.rep[j].operate(sr);
+  }
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
